Highlight
Free Access
Issue
A&A
Volume 551, March 2013
Article Number A43
Number of page(s) 23
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201219657
Published online 18 February 2013

© ESO, 2013

1. Introduction

Does the formation of a massive star differ significantly from that of a low-mass star? Arguably, all studies of high mass (≥ 8 M) star formation are centred upon this question. There are several possible reasons to expect differences at higher masses, one being that a massive star is thought to continue to accrete material after reaching the zero-age main sequence (ZAMS, e.g. Hosokawa et al. 2010), which is a consequence of its Kelvin-Helmholtz contraction time-scale being shorter than its accretion time-scale. Therefore, several processes such as radiation pressure (e.g. Yorke 2002) and ionization (e.g. Keto 2002) may halt, decrease or alter accretion on to the star. However, in the earlier stages of protostellar evolution, Hosokawa et al. (2010) also find that the central accreting stars in their simulations become bloated due to accretion, so that the effective temperature and UV luminosity of the protostar remains low. Hence at earlier times these effects may be reduced. Secondly, while the densities and temperatures in pristine molecular clouds can easily explain the formation of low-mass star-forming cores, without the input of some stabilising energy such as an external pressure, micro-turbulence (McKee & Tan 2003), magnetic fields (e.g. Hennebelle et al. 2011; Commerçon et al. 2011) or radiative heating and outflows (e.g. Krumholz et al. 2012), these conditions are not conducive to creating a “monolithic” core that does not fragment, from which the forming massive star can accrete all of its mass. Instead, massive stars may start off embedded in smaller cores that source most of their mass from the surrounding cluster-forming clump (Bonnell et al. 2011; Myers 2011), implying that massive stars can only form in clusters, not in isolation.

Placing these theoretical concerns aside initially, and working under the hypothesis that massive stars form as a scaled-up version of low-mass star formation, in this paper we aim to probe the circumstellar environment of the massive star-forming region AFGL 2591 by combining multi-wavelength observations and modelling to determine whether any of the signatures of its formation differ significantly from those of low-mass protostars.

AFGL 2591 is a well-studied example of a luminous star-forming region (2.1−2.5    ×    105   L at 3.33 kpc, Lada et al. 1984; Henning et al. 1990; Rygl et al. 2012). One of its most prominent features is a one sided conical reflection nebula observed in the near-IR (e.g Minchin et al. 1991; Tamura & Yamashita 1992). Projected within the Cygnus-X star-forming complex, the distance to the source has recently been determined by trigonometric parallax measurements to be 3.33    ±    0.11 kpc, more distant than previously assumed (between 1 and 2 kpc, e.g. Poetzel et al. 1992; Hasegawa & Mitchell 1995; Trinidad et al. 2003; van der Tak & Menten 2005). As will be described below, AFGL 2591 actually consists of several objects. However, as one source, AFGL 2591-VLA 3, dominates the spectral energy distribution (SED) and infrared images and hence the luminosity, the name AFGL 2591 will therefore also be used henceforth to refer to this dominant source.

AFGL 2591 has been studied and modelled by many authors. For example, one-dimensional modelling of the circumstellar geometry via the observed SED has been carried out by Guertler et al. (1991); van der Tak et al. (1999); Mueller et al. (2002b) and de Wit et al. (2009). Improving on these, Preibisch et al. (2003) used two models – one of a disk and the other of an envelope with outflow cavities – to reproduce the 40 AU diameter (assuming d = 1 kpc) bright disk of emission observed in their K band image, and Trinidad et al. (2003) have modelled the millimetre emission from AFGL 2591-VLA 3 as an optically thick disk without an envelope. As part of their comprehensive study, van der Tak et al. (1999) modelled their observed molecular lines (CS, HCN, HCO+) including an outflow cavity in a power-law envelope, finding that a half-opening angle of 30° was able to better reproduce the line profiles. van der Wiel et al. (2011) have modelled a set of six lines detected toward AFGL 2591 as part of the JCMT Spectral Legacy Survey, finding evidence of a cavity or inhomogeneity in the envelope on scales of  ≤ 104AU. In addition, as well as observing a bi-conical outflow structure in 12CO (2–1) on scales of 1−2′′ or several thousand AU, Jiménez-Serra et al. (2012) also uncovered evidence for chemical segregation within the inner 3000 AU of the AFGL 2591-VLA 3 envelope (assuming d = 3 kpc, similar to our assumed distance of 3.33 kpc).

The ionized gas emission in the region surrounding AFGL 2591 has previously been observed by Campbell (1984b); Tofani et al. (1995); Trinidad et al. (2003) and van der Tak & Menten (2005) from 5 to 43 GHz. These observations showed that AFGL 2591 is in fact not an isolated forming star, uncovering four continuum sources in the region. The observed fluxes of two of these, VLA 1 and VLA 2, gave spectral indices consistent with optically thin free-free emission from Hii regions, and a third source VLA 3, was measured to have a steeper spectral index, possibly indicating optically thick emission. VLA 3 is also coincident with the central illuminating source of AFGL 2591 observed at shorter wavelengths (Trinidad et al. 2003). A fourth radio continuum source has also been detected by Campbell (1984a) and Tofani et al. (1995) (their source 4 and n4 respectively), which shall be referred to as VLA 4 in the following sections. In addition, the 3.6 cm images presented in this work uncover a fifth source, VLA 5, to the south west.

Table 1

Summary of VLA observations.

Knots of H2 and [SII] Herbig Haro objects have been detected toward AFGL 2591 (Tamura & Yamashita 1992; Poetzel et al. 1992), suggesting the presence of shocked gas. These are coincident with an east-west bipolar outflow (e.g. Lada et al. 1984; Hasegawa & Mitchell 1995), which extends across 5′ or 4.8 pc at 3.33 kpc, but also contains a more collimated central small-scale component with an extent of 90′′    ×     20′′.

Evidence for the presence of a disk or rotationally flattened material around the central source of AFGL 2591 has been found by several authors. The near-IR imaging polarimetry observations of Minchin et al. (1991) showed that a disk or toroid of material was needed to appropriately scatter the emission. In addition, a large (50′′    ×     80′′ ) flattened “disk” of material has been observed perpendicular to the outflow in observations of CS lines (Yamashita et al. 1987). At smaller scales, van der Tak et al. (2006) and Wang et al. (2012) found evidence for a disk of diameter 800 AU at 1 kpc (corresponding to 2700 AU at 3.33 kpc), which exhibits a systematic velocity gradient in the northeast-southwest direction. This gradient was also found in the SMA observations of Jiménez-Serra et al. (2012), which they found to be consistent with Keplerian-like rotation around a 40 M star. Finally, both OH and water masers have been observed towards VLA 3 (Trinidad et al. 2003; Hutawarakorn & Cohen 2005; Sanna et al. 2012). The Very Large Array (VLA)1 22 GHz water maser observations of Trinidad et al. (2003) uncovered a maser cluster, which included a  ~0.01′′ diameter shell-like structure on the smallest scales. As well as finding largely consistent results, the Very Long Baseline Array 22 GHz water maser observations of Sanna et al. (2012) determined that the maser cluster is arranged in a v-shape, which coincides with the expected location of the outflow walls apparent in the near-IR reflection nebula.

In this paper, we have adopted a multi-wavelength approach to further probe the circumstellar environment of AFGL 2591, and to address the question of whether it forms in a similar manner to its low-mass counterparts. Building on previous work, the modelling presented in this paper includes the first simultaneous radiative transfer model of the near-IR through sub-mm SED, and near-IR images, of AFGL 2591 with a three-dimensional axisymmetric geometry. In addition, we present new multi-configuration VLA 3.6 cm and 7 mm continuum observations that for the first time trace an ionized jet at 3.6 cm, as well as 13CO(1−0), C18O(1−0) and 3 mm continuum Combined Array for Research in Millimeter-wave Astronomy (CARMA2) observations, which trace previously unstudied scales within the molecular outflow and envelope, to derive a self-consistent picture of AFGL 2591.

Section 2 outlines the new observations carried out in this work, describing both the VLA 3.6 cm and 7 mm, as well as the CARMA 13CO, C18O and 3 mm continuum observations. Section 3 presents the archival data used in this paper, including near-IR to sub-mm SED, 2MASS photometry, measured 2MASS brightness profiles and Gemini North near-IR images of AFGL 2591. Section 4 describes the modelling of the near-IR to sub-mm SED and near-IR images. Section 5 presents the results for the SED and near-IR image modelling, as well as for the centimetre and millimetre wavelength datasets. Section 6 presents our discussion, which covers the topics of the jet and outflow of VLA 3, and the properties of AFGL 2591 as a cluster. Our conclusions are given in Sect. 7.

2. Radio interferometric observations

2.1. VLA 3.6 cm and 7 mm continuum

Multi-configuration radio continuum observations at 3.6 cm and 7 mm were conducted between April 2007 and March 2008 with the VLA of the National Radio Astronomy Observatory3. During this time frame, the VLA continuum mode consisted of four 50 MHz bands, two of which were placed at 8.435 GHz, and two at 8.485 GHz for 3.6 cm, and similarly two at 43.315 GHz and two at 43.365 GHz for 7 mm. At 3.6 cm, observations were taken in all four configurations of the VLA (A to D), and at 7 mm, observations were performed in B through D array and combined with A array archive data previously published in van der Tak & Menten (2005). When combined, these observations provided baseline lengths between 35 m and 36.4 km, giving information on angular scales from 0.24′′ to 3′ for 3.6 cm, and 0.05 to 43′′ for 7 mm.

For each observation, Table 1 lists the observed wavelength, configuration, observation date, program code, number of antennas in the array (with the number of antennas with useful data shown in brackets), time on-source, the pointing centre of the target, and the flux density of the gain calibrator determined from bootstrapping the flux from the primary flux calibrator. The number of antennas with useful data was reduced due to antennas being out of the array for upgrading and testing of new VLA capabilities, general hardware issues, and also some residual bad data. For all observations, the gain calibrator was 2015 + 371 and primary flux calibrator was 1331 + 305 (3C286).

Data reduction and imaging were carried out using the Common Astronomy Software Applications (CASA)4 package. As 1331 + 305 was slightly resolved, a model was used for flux calibration, which has a flux of 5.23 Jy at 3.6 cm and 1.45 Jy at 7 mm, allowing data at all uv distances to be used. Antenna gain curves and an opacity correction (zenith opacity = 0.06) were applied for the 7 mm data. The error in absolute flux calibration is approximately 1−2% for the 3.6 cm band, and 3−5% for the 7 mm band. The data were imaged using a CLEAN multi-scale deconvolution routine, using three scales, and Briggs weighting with a robust parameter of 0.5. The rms noise in the final combined images were 30 and 56 μJy beam-1 for the 3.6 cm and 7 mm images respectively, and the synthesised beams were 0.43′′ by 0.40′′, PA = 43° and 0.11′′ by 0.11′′, PA = 43°. These values, as well as the synthesised beam and map rms for each part of the CARMA observations described in Sect. 2.2, are summarised in Table 2.

2.2. CARMA 13CO, C18O and 3 mm continuum

Table 2

Summary of radio interferometric observations.

CARMA observations at λ ~ 3 mm were taken on 10 July 2007 in D configuration (corresponding to baseline lengths between 11−150 m). A total of 14 antennas were used during the observations, five 10-m and nine 6-m in diameter.

The CARMA correlator was comprised of two side bands, placed either side of the chosen LO frequency, which for these observations was 108.607 GHz. Both the upper and lower side bands contained three spectral windows: one wide and two narrow, giving a total of six windows. The two wide spectral windows, centred on 106.7 and 110.5 GHz, had a bandwidth of 500 MHz and a total of 15 channels, and the four narrow spectral windows had a width of approximately 8 MHz and 63 channels, giving a spectral resolution of 122 kHz. The observed lines, 13CO(J = 1−0) and C18O(J = 1−0), lay in the two narrow spectral windows in the upper side band at 110.201 and 109.782 GHz. At these frequencies, the spectral resolution was  ~0.33 km s-1.

The phase centre for AFGL 2591 was 20h29m24.s\hbox{$\fs$}7 + 40°11′18.\hbox{$\farcs$}9 (J2000), and the total time on source was 2.5 h. MWC 349 was used as a gain and bandpass calibrator, and Neptune as a flux calibrator. Data calibration and imaging were carried out in CASA. To do this, the data was exported to CASA from its original MIRIAD format after applying line length corrections and Hanning smoothing. Using a flux for Neptune of  ~4.2 Jy at 108 GHz (Kwon, priv. comm.), fluxes of 1.2 and 1.3 Jy for MWC 349 were derived for the upper and lower side bands respectively. As there was insufficient signal-to-noise on MWC 349 in the narrow spectral windows, band pass calibration was only performed on the two wide spectral windows. However, as the amplitude bandpass shape varies across the 8 MHz band by less than 3%, and similarly for the phase (Friedel, priv. comm.), this should not significantly affect our results. To image the data, Briggs weighting and a robust parameter of 0.5 were used. The beam size and rms noise for each observed line and continuum spectral window are given in Table 2.

3. Archival data

3.1. SED

The observed near-IR to sub-millimetre SED of AFGL 2591 is shown in Fig. 1. The SED data points were collected from the literature; Table 3 lists the wavelengths, flux densities and references for the data displayed in Fig. 1. The 30 data points shown between 3.5 and 45 μm were sampled uniformly in log-space from a highly processed ISO-SWS spectra of AFGL 2591 (Sloan et al. 2003) taken on 7 Nov. 1996, observation ID 35700734. The uncertainties for the averaged ISO-SWS spectrum were assumed to range from 4 to 22% as per the ISO handbook volume V5. As mentioned in Sect. 1, the SED and therefore flux densities of AFGL 2591 given in Table 3 were assumed to be dominated by one source. Evidence to support this assumption will be presented in Sect. 5.2.2, via a comparison of the SEDs of the resolved sources in the region.

thumbnail Fig. 1

SED of AFGL 2591, collated from the literature. The best-fitting models for an envelope with and without a disk (overplotted solid blue and dashed red lines respectively) are discussed in Sect. 5.1. The errors shown are those reset to 10% if the error in the measured flux was  <10%.

Table 3

Observed near-IR to sub-millimetre fluxes for AFGL 2591, collated from the literature.

3.2. Near-IR 2MASS images

A three-colour 2MASS (Skrutskie et al. 2006) J, H, and Ks band image of the near-IR emission observed towards AFGL 2591 is shown in the middle panel of Fig. 2, where red, green and blue are Ks, H and J bands respectively. This Figure shows the conical near-IR reflection nebula (previously observed by e.g. Kleinmann & Lebofsky 1975; Minchin et al. 1991; Tamura & Yamashita 1992), which is consistent with being a blue-shifted outflow lobe from the central source of AFGL 2591.

To obtain integrated J, H and Ks band fluxes for AFGL 2591, irregular aperture photometry was performed on 2MASS Atlas images. Photometric uncertainties on the fluxes are  <1%, however it is likely the fluxes are more uncertain due to the choice of the irregular photometry apertures, and foreground or background stars being included within them. The uncertainty due to contaminating stars is likely dominated by the 2MASS source 20292393 + 4011105 at 20h29m23.s\hbox{$\fs$}93 + 40°11′10.\hbox{$\farcs$}56 (J2000; easily visible as a point source in the Gemini North near-IR images presented in Sect. 5.2), as the apertures were placed to avoid all other point sources visible in the 2MASS images. The fluxes of this contaminating source are 0.011, 0.035 and 0.105 Jy at J, H and Ks bands respectively, with Ks band being an upper limit. Therefore, the uncertainty in the measured 2MASS fluxes of the AFGL 2591 outflow lobe was estimated to be  ~10%.

3.3. 2MASS brightness profiles

Flux profiles of AFGL 2591, summed across strips aligned with and perpendicular to the source outflow axis (PA = 256°, taken from Sanna et al. (2012), with thicknesses of 24.6′′ and 104′′ respectively), were measured for the three 2MASS bands. The background subtracted, normalised profiles are shown in Fig. 3 along with the best-fitting models to both the SED and profiles, for the models consisting of an envelope with and without a disk (blue solid and red dashed lines respectively), which will be discussed further in Sect. 5.1. The background was determined by finding the average value within two strips either side of the main profile. The errors in the profiles shown in each panel of Fig. 3 reflect the uncertainty due to background fluctuations, and are calculated as the standard deviation of the profiles measured in the two background strips, which were assumed to contain minimal source flux. As the strips used to calculate the uncertainty in the perpendicular profiles contained several bright stars, iterative sigma clipping was performed, i.e. pixels lying more than 5 standard deviations away from the mean were then ignored and the mean and standard deviation were recalculated and the process was repeated until no pixels were rejected.

thumbnail Fig. 2

a) Model three-colour J, H and Ks band image for the envelope with disk model (RGB: K, H, and J bands respectively). b) Observed 2MASS three-colour image of AFGL 2591. c) Model J, H, and Ks band image for the envelope without disk model. The model images have been normalised to the total integrated fluxes given in Table 3, in order that the morphology of the emission can be easily compared. Stretch: red: Ks band, 130−300 MJy sr-1; green: H band, 95−150 MJy sr-1, blue: J band, 30−60 MJy sr-1.

3.4. Near-IR gemini north NIRI images

High resolution near-IR images of AFGL 2591 were obtained from the Gemini website, which were taken under program GN-2001A-SV-20 as part of the commissioning of Gemini North’s NIRI instrument, and were available under public release6. The total integration time was 2 min for J band and 1 min at H and K′ bands. We aligned the images using Chandra X-ray sources in the field (Evans et al. 2010), giving a positional accuracy of 0.6′′, and calibrated the flux level of the images in MJy sr-1 using aperture photometry of bright sources in both the Gemini North and 2MASS images. The FWHMs of the PSFs in the images range between 0.3−0.4′′. The peak position of the central source of AFGL 2591 in the J band image was found to be 20h29m24.s\hbox{$\fs$}86 + 40°11′19.\hbox{$\farcs$}5 (J2000), which was taken to be the position of the powering object. As these images were saturated in K′ band, 2MASS images were instead used for the modelling described in Sect. 4. However, in Sects. 5.2 and 5.3 these high-resolution images are compared to the radio continuum and molecular line data observed for this work.

4. SED and near-IR image modelling

4.1. The dust continuum radiative transfer code

The SED of AFGL 2591 was modelled between near-IR and sub-mm wavelengths using the 3D Monte Carlo dust radiative transfer code Hyperion described in Robitaille (2011) and the genetic search algorithm used to model IRAS 20126 + 4104 in Johnston et al. (2011).

For a given source and surrounding dust geometry sampled in a spherical-polar grid, Hyperion models the nonisotropic scattering and thermal emission by/from dust, calculating radiative equilibrium dust temperatures (using the technique of Lucy 1999), and producing spectra and multi-wavelength images. As dust and gas are assumed to be coupled in the code we are able to probe the bulk of the material which surrounds AFGL 2591. In this section, we describe the density structure of the disk, envelope and outflow cavity in the model.

We model the circumstellar geometry of AFGL 2591 with the same 3D axisymmetric envelope with disk geometry successfully employed to model the SEDs and scattered light images of low-mass protostars (e.g. Robitaille et al. 2007; Wood et al. 2002; Whitney et al. 2003). The axisymmetric three dimensional flared disk density is described between inner and outer radii Rmin and Rdisk by ρdisk(ϖ,z)=ρ0(R0ϖ)αexp12[zh(ϖ)]2,\begin{equation} \rho_{\rm disk} (\varpi,z)=\rho_0 \left (\frac{R_0}{\varpi}\right )^{\alpha} \exp{ \left\{ -{1\over 2} \left[\frac{z}{h( \varpi )}\right]^2 \right\} } , \label{discdensity} \end{equation}(1)where R0 = 100 AU, ϖ is the cylindrical radius, z is the height above the disk midplane and ρ0 is set by the total disk mass Mdisk, by integrating the disk density over z, ϖ the cylindrical radius, and φ the azimuthal angle. The scaleheight h(ϖ) increases with radius: h(ϖ)=h0ϖ/R0)(β\hbox{$h(\varpi)=h_0\left ( {\varpi /{R_0}} \right )^\beta$} where h0 is the scaleheight at R0. We have assumed the parameter β = 1.25, derived for irradiated disks around low mass stars (D’Alessio et al. 1998), and have taken α = 2.25, to provide a surface density of Σ ~ ϖ-1. This equation is essentially the same as Eq. (1) given in Johnston et al. (2011), however their second term is removed as this pertains to accreting α-disks, whereas our disk is not accreting. Further, Eq. (1) is expressed in terms of R0 instead of R so that the scaleheight is defined at a radius which is more conceptually accessible.

The density of the circumstellar envelope is taken to be that of a rotationally flattened collapsing spherical cloud (Ulrich 1976; Terebey et al. 1984) with radius Renvmax\hbox{$R^{\rm max}_{\rm env}$}, ρenv(r,μ)=ρenv0(rRc)3/2(1+μμ0)1/2(μμ0+2μ02Rcr)-1\begin{equation} \rho_{\rm env} (r,\mu) = \rho_{{\rm env}0} \left( {{r}\over{R_{\rm c}}}\right)^{-3/2}\left(1+{{\mu}\over{\mu_0}} \right)^{-1/2} \left( {{\mu}\over{\mu_0}}+{{2\mu_0^2R_{\rm c}}\over{r}}\right)^{-1} \end{equation}(2)where ρenv0 is the density scaling of the envelope, r is the spherical radius, Rc is the centrifugal radius, μ is the cosine of the polar angle (μ = cosθ), and μ0 is the cosine of the polar angle of a streamline of infalling particles in the envelope as r → ∞. The equation for the streamline is given by μ03+μ0(r/Rc1)μ(r/Rc)=0\begin{equation} \mu_0^3 + \mu_0(r/R_{\rm c}-1)-\mu(r/R_{\rm c})=0\; \end{equation}(3)which can be solved for μ0. In our model, we assume the centrifugal radius is also the radius at which the disk forms Rc    =    Rdisk.

We chose to describe the envelope density via a density scaling factor instead of the envelope accretion rate and stellar mass used in Johnston et al. (2011), as this does not require the assumption of evolutionary tracks – which are currently very uncertain for massive stars – to provide a physical stellar radius and temperature for a given stellar mass. Instead, we have varied the stellar radius and temperature as parameters (described in Sect. 4.2).

To reproduce the morphology of the emission in the near-IR images, we also include a bipolar cavity in the model geometry with density ρcav, and a shape described by z(ϖ)=z0ϖbcav\begin{equation} z(\varpi) = z_0 \varpi^{{\rm b}_{\rm cav}} \label{cavity} \end{equation}(4)where the shape of the cavity is determined by the parameter bcav which we set to be bcav = 1.5, and z0 is chosen so that the cavity half-opening angle at z = 10   000 AU is θcav. The cavity was included by resetting the envelope density to ρcav inside the region defined by Eq. (4) only where the envelope density was initially larger.

The inner radius of the dust disk and envelope Rmin is expressed in terms of the dust destruction radius Rsub, which was found empirically by Whitney et al. (2004) to be Rsub=R(Tsub/T)-2.1\begin{equation} R_{\rm sub} = R_{\star}(T_{\rm sub}/T_{\star})^{-2.1} \end{equation}(5)where the temperature at which dust sublimates is assumed to be Tsub = 1600 K, and T is the temperature of the star. We adopted the dust opacity and scattering properties from Kim et al. (1994) for the circumstellar dust. This dust model has a ratio of total-to-selective extinction RV = 3.6 and a grain size distribution with an average particle size slightly larger than the diffuse interstellar medium.

For the genetic algorithm, which is described in detail in Sect. 4.3 of Johnston et al. (2011), the size of the first generation N was set to 1000, and the size of subsequent generations M was set to be 200. The code was taken to be converged when the χ2 value of the best fitting model decreased less than 5% in 20 generations.

thumbnail Fig. 3

Black error bars: normalised flux profiles for the three 2MASS bands, aligned with (top) and perpendicular to (bottom) the outflow axis. Blue and red lines: the profiles of the best-fitting models to both the observed SED and profiles, for envelope with disk and without disk models respectively. Vertical grey dashed lines mark the position of the central source.

Table 4

Assumed ranges for model parameters as input to the genetic search algorithm.

4.2. Input assumptions: model parameters

A set of plausible ranges for parameters describing the model, in which the genetic algorithm searched for the best fit, is given in Table 4. In addition to this, the stellar radius and temperature were also required to lie above the zero-age main sequence defined as log TZAMS vs. log (0.9    ×    RZAMS) for the solar metallicity models of Schaller et al. (1992), where the factor of 0.9 is to ensure that enough models are sampled around the ZAMS, rather than strictly above the ZAMS. The stellar temperature T and a stellar surface gravity of log g = 4 were then used to select a model from the stellar atmosphere grid of Castelli & Kurucz (2004). To sample the cavity half-opening angle θcav, we first determined the relation between the half-opening angle and the inclination such that the equation of the cavity, Eq. (4), projected onto the plane of the sky produced the right opening angle observed in the Gemini near-IR images, i=arcsin[0.826 (sinθcav)bcosθcav]·\begin{equation} i = \arcsin{\left[ \frac{ 0.826~(\sin{\theta_{\rm cav}})^b}{ \cos{\theta_{\rm cav}}}\right]}\cdot \end{equation}(6)When sampling the half-opening angle, the inclination was first sampled and Eq. (6) was then used to determine an upper bound for θcav, as radiative transfer effects only work to make the cavity look larger than it actually is. For instance, light can be scattered into the envelope making it look larger, whereas because the density of the cavity is constant, there is no definite edge of emission seen within the cavity. The largest half-opening angle possible, for an inclination of 90°, was therefore 54° at a radius of 10 000 AU.

4.3. Model fitting

As part of the genetic algorithm used here and described in Johnston et al. (2011), the models were fit to the data using the SED fitting tool of Robitaille et al. (2007), where Av, the external extinction, was a free parameter in the fit. The distance was set to be 3.33 kpc, and the visual extinction Av was allowed to vary between 0 and 40 mag.

As described in Robitaille et al. (2007), the fitting automatically takes into account the different circular aperture sizes when fitting the models. However, in the case of the ISO fluxes, the ISO apertures are not circular, so instead the fluxes were determined by computing images at the ISO wavelengths, and summing up the flux weighted by the slit transmission maps (found within the Observer’s SWS Interactive Analysis software, OSIA7) which give the transmission as a function of location on the detector. The ISO slits were rotated to the correct orientation for the observation and applied to the model. During the SED fitting, each data point was weighted by the distance in wavelength between its two adjacent data points. This was done to allow the algorithm to search for a model which reproduced the shape of the entire SED, not only the sections with a high density of data.

The models were also simultaneously fit to the six flux profiles shown in Fig. 3, measured from the 2MASS images. The model flux profiles were found in the same way as the observed profiles, by summing the flux in the model images within 24.6′′ and 104′′ wide strips aligned with and perpendicular to the outflow axis, and normalising them to the total flux. To produce the convolved model J, H and Ks band images, the model images were convolved with the FWHM of the 2MASS point-spread functions (PSF) for each band: 3.15′′, 3.23′′ and 3.42′′ for J, H and Ks bands respectively, which were determined by fitting Gaussian PSFs to stars in the observed images. The model images shown in Fig. 2 are high signal-to-noise versions of the images used for fitting.

For both the SED and profiles, the errors were reset to 10% before fitting if the uncertainty in the measured flux was  <10%. This was done to take into account non-measurement errors such as variability. The combined reduced χ2 for each model given as input to the genetic algorithm was calculated as χcombined2=χSED2+χprofiles2\begin{equation} \chi^2_{\rm combined}= \chi_{\rm SED}^2 + \chi_{\rm profiles}^2 \end{equation}(7)where χSED2\hbox{$\chi_{\rm SED}^2$} is the best-fitting reduced χ2 value for the SED alone, and χprofiles2\hbox{$\chi_{\rm profiles}^2$} is the overall reduced χ2 for the profile fits.

5. Results

5.1. Results of SED and image modelling

In this paper, we model AFGL 2591 as a single source; as we will show in Sect. 5.2.2, the second brightest object in the region VLA 1 does not significantly contribute to its infrared flux, therefore this is a reasonable approximation.

The genetic code was run three times – firstly with the model parameter ranges given in Table 4 as input for all parameters, which will be referred to as the “envelope with disk” model, and secondly with the parameter ranges given in Table 4 for all parameters except the disk mass and disk scale height at 100 AU, Mdisk and h0, which were set to zero and were therefore not treated as model parameters. This second run will be referred to as the “envelope without disk” model. The envelope without disk model was run to ascertain whether the SED and images could be adequately produced without a disk, i.e. a simpler model which had two fewer parameters. For the envelope without disk model, Rdisk will instead be referred to as the centrifugal radius of the envelope, Rc, as these two parameters are interchangeable in the models. The third run had exactly the same setup as the envelope with disk model, but was used to determine the repeatability of the experiment and any parameter degeneracies, and will be referred to as the control model.

The genetic search algorithm codes were stopped after 28, 29 and 49 generations for the envelope with disk, without disk and control models respectively, when the convergence criterion was reached. This corresponds to 6545 models or 2940 CPU hours for the envelope with disk model, 6278 models or 3774 CPU hours for the envelope without disk model, and 10536 models or 3960 CPU hours for the control model.

Table 5

Parameters of the genetic algorithm best-fitting models.

The resulting SEDs and profiles of the best-fitting models for the envelope with and without disk runs are shown against the data in Figs. 1 and 3, and the model images are compared to the observed images in Fig. 2. The parameters of the best-fitting models (i.e. the model in each run with the lowest overall χ2 value after convergence), are also given in Table 5. The interstellar extinction AV for the best-fitting models of the envelope with disk, without disk and control runs was 0.0, 0.0 and 0.8 respectively.

The minimum reduced χ2 for the best-fitting envelope with disk model is 14.9, and the corresponding value for the best-fitting without disk model is 18.7, therefore the model with a disk provides a marginally better fit to the SED and profiles. Comparing the model SEDs to the data shown in Fig. 1, it can be seen that they are qualitatively very similar, although the model with a disk fits the SED slightly better. This difference is partially caused by their disagreement in the mid-IR regime, where the without disk model does not follow the data as closely. Figure 3 also shows very similar fits to the 2MASS image profiles for each model, however the without disk model fits some of the profiles marginally better, being more centrally peaked. We note that although the reduced χ2 for the best-fitting model with a disk is lower, and thus provides evidence that an envelope with an embedded disk describes the source better than a model without a disk, the difference is not large enough to prove this conclusively. We note that neither model provides a good fit to the J and H band profiles along the outflow axis. This is likely to be caused by the fact that the emission at these wavelengths is not dominated by the central source, but instead by scattered light in the surrounding medium, which in reality is highly non-uniform, as seen in Figs. 2 and 7. Thus to reproduce the images more accurately, future modelling will have to take such inhomogeneities into account, for example the “loops” in the outflow cavity seen in the near-IR images (Preibisch et al. 2003).

To understand how well the parameters of AFGL 2591 were determined, we took the best-fitting models for each run and varied in turn each parameter across the original parameter ranges, while holding the other parameters constant. We then fit the resultant SED and images for each new model, and calculated their overall reduced χ2-values. From these we could therefore determine the χ2 surfaces in the vicinity of the the best-fitting model parameters and hence understand their uncertainties. Figure 4 shows such a plot for each varied parameter, for all three runs. Firstly, we can see the best-fitting values and shape of the χ2 surfaces for stellar radius and temperature are different for each run, which is due to a degeneracy between these parameters. The χ2 surfaces for the stellar luminosity, which is a combination of these two parameters, is shown as a panel in Fig. 4: this surface is very well defined and very similar for all runs, showing that fitting the SED allows this parameter to be accurately and non-ambiguously determined. In the next panel we show the χ2 surface for the parameter orthogonal to the stellar luminosity, R2T\hbox{$R_{\star}^{2}T_{\star}$}, is flat and thus poorly determined. Therefore the SED or image profiles do not provide a handle on either the stellar radius and temperature, only the stellar luminosity. Similarly, the envelope density scaling factor ρenv0 and Rdisk or Rc are degenerate, as they both determine the envelope density, which can be seen from Eq. (2). The final two panels of Fig. 4, show the orthogonal parameter combinations ρenv Rc3/2\hbox{$\rho_{\rm env}~R_{\rm c}^{3/2}$} and ρenv Rc2/3\hbox{$\rho_{\rm env}~R_{\rm c}^{-2/3}$}. Here ρenv Rc3/2\hbox{$\rho_{\rm env}~R_{\rm c}^{3/2}$} is similar between the three runs, and has a sharp χ2 surface, while the best-fitting ρenv Rc2/3\hbox{$\rho_{\rm env}~R_{\rm c}^{-2/3}$} value changes by an order of magnitude. Thus we cannot uniquely determine these two parameters from our modelling, but we can determine a parameter that combines these which governs the overall envelope density and mass, ρenv Rc3/2\hbox{$\rho_{\rm env}~R_{\rm c}^{3/2}$}.

Looking at the remaining parameters, it can be seen that although the minima in their χ2 surfaces are not sharp, the shape of their χ2 surfaces is similar for the three runs, with Rmin tending to smaller values, Renvmax\hbox{$R^{\rm max}_{\rm env}$} tending to intermediate or large values, and the inclination preferring intermediate values. The χ2 surface for the cavity half-opening angle θcav changes between each run, however is generally poorly constrained. It is also important to note that the half-opening angle is somewhat degenerate with the inclination, as the inclination of a model determines the possible upper limit for the half-opening angle as described in Sect. 4.2. The shape of the χ2 surface for the cavity density ρcav is variable at higher values, but always has its minimum at low values, between 10-21 and  ~2    ×    10-20 g cm-3\hbox{${}^{-3}$}. For the two runs which include a disk, the χ2 surfaces of the remaining parameters Mdisk and h0 are almost flat across the sampled ranges, although they slightly prefer values close to the upper limit. However, for h0 the χ2 rises steeply toward the upper limit of the parameter range after its minimum at 11 AU.

thumbnail Fig. 4

Plots of the χ2 surface along the axis of each model parameter for each run, determined by varying each parameter while holding the other best-fitting parameters constant. Red dashed line: envelope without disk run, light blue line: envelope with disk run, and black line: control run.

In the remainder of this section, we compare the results of our modelling to the recent results from other studies of AFGL 2591. Firstly, our modelling finds values of roughly 2−5 Rsub for the envelope and disk inner radius Rmin. Using Eq. (5), we find Rsub = 34−39 AU for the three best-fitting models, and hence we expect values of Rmin ~ 70–200 AU. This value is in agreement with the radius of the inner cavity adopted by Jiménez-Serra et al. (2012), who in turn took this value from the disk of emission observed by Preibisch et al. (2003), likely to be tracing the inner rim of the envelope and disk. At 1 kpc, Preibisch et al. (2003) found the radius of this emission was 40 AU, which is thus 130 AU at 3.33 kpc.

We compare the maximum envelope radius Renvmax\hbox{$R^{\rm max}_{\rm env}$} to the results of the JCMT line observations of van der Wiel et al. (2011), which were able to trace the largest envelope scales, with the maximum radius in their model set to 200 000 AU (scaled to 3.33 kpc). This is in agreement with our SED modelling results, which prefer large envelope radii, on the scale of several hundreds of thousand AU.

We find the cavity half-opening angle θcav to be poorly constrained by our modelling, although we determined an upper limit for the cavity half-opening angle of 54° in Sect. 4.2 (defined at 104 AU). Wang et al. (2012) adopted an opening angle of 60°, and hence a half-opening angle of 30°, in their modelling. The models used in Jiménez-Serra et al. (2012) and van der Wiel et al. (2011) do not include an outflow cavity.

We find the inclination to prefer intermediate values, between roughly 30−65°. This is in agreement with the adopted inclination of 30° by Wang et al. (2012), who in turn used the value determined by van der Tak et al. (2006). However, Jiménez-Serra et al. (2012) found an inclination of 20° by modelling the position–velocity diagram of combined emission from several lines. Thus, previous results appear to favour an inclination which is slightly closer to face-on. As the inclinations for the envelope with and without disk models were 59 and 46° respectively, we therefore adopted an inclination of 50° to determine several of the physical properties of the jet and outflow in Sect. 6.1, as a compromise between these two values.

We found low cavity densities ρcav were preferred for all runs, with best-fitting values between 10-21 and  ~2    ×    10-20 g cm-3\hbox{${}^{-3}$}. In Wang et al. (2012) and van der Tak et al. (1999), their model cavity densities were set to zero.

Our modelling found the mass of the disk to be on the order of 10−20 M, whereas using 0.5′′ 203.4 GHz dust continuum observations Wang et al. (2012) determined the disk mass to be  ~1−3 M. However, as the disk mass was one of the most poorly determined parameters by our modelling, and as the interferometric dust continuum observations of Wang et al. (2012) directly measure the column density and thus mass of the disk, without the intervening emission or scattering from the envelope, we suggest that the result of Wang et al. (2012) is more robust.

It was not possible to compare our best-fitting values for the disk scale-height h0 to the model of Wang et al. (2012), as they instead used a sin5θ function to describe the vertical structure of the disk, where θ is the angle from the polar axis, as opposed to our Gaussian profile with z.

The best-fitting luminosity of AFGL 2591 is 2.3    ×    105   L for both the envelope with and without disk models respectively, which is in agreement with previous results (2.1−2.5    ×    105   L at 3.33 kpc, Lada et al. 1984; Henning et al. 1990; Rygl et al. 2012).

The envelope density parameter ρenv0 Rc3/2\hbox{$\rho_{\rm{env}0}~R_{\rm c}^{3/2}$} is a proxy for the mass of the envelope, which dominates the total mass of the enclosed circumstellar material at large radii, as it describes the scaling of Eq. (2). We find the total mass of the envelope plus disk is 1200 M for the envelope plus disk model, and 2700 M for the envelope without disk model. The total mass of gas associated with the region observed by Minh & Yang (2008) is 2    ×    104   M, which is scaled to 3.33 kpc and within a radius of 5.8 pc. Thus although our model traces a sizeable fraction of the envelope material within a radius of several pc, it may miss more diffuse material on the scales observed by Minh & Yang (2008), and thus could underestimate the true total mass associated with the region. To determine whether the mass distribution on smaller scales of our best-fitting models was consistent with other results, we compared the total mass within a radius of 2700 AU to the model of Jiménez-Serra et al. (2012). Within this radius, which is half the determined size of their envelope 5400 AU, the total mass of the model of Jiménez-Serra et al. (2012) is 1 M. Comparing this to our model we find an enclosed mass within 2700 AU of 0.1 M for the model without a disk, and 1 M for the model with disk. Therefore in the case where there is a disk, the enclosed masses are in agreement.

5.2. 3.6 cm and 7 mm continuum

thumbnail Fig. 5

a) Map of the 3.6 cm continuum emission surrounding AFGL 2591. Contours are -3, 3, 4, 5, 7, 10, 15, 20, 30, 40, 50...100  ×  rms noise = 30 μJy beam-1. Greyscale: –0.03 to 3.77 mJy beam-1 (1.2  ×  peak value). The synthesised beam is shown in the bottom left corner: 0.43′′    ×    0.40′′, PA = 43°. b) Close-up of the 3.6 cm continuum emission towards VLA 3. Contours, greyscale and beam as in panel a).

Figure 5 presents the observed multi-configuration 3.6 cm image of AFGL 2591. Panel a) shows the entire region surrounding the source, including the four sources first observed by Campbell (1984a) at 6 cm (VLA 1 through 4). In addition, a new, low surface brightness source was detected: VLA 5, which varies in peak flux from 3−5 σ across the extent of the source. With the addition of shorter baselines, more extended emission is recovered compared to the 3.6 cm A array images of Tofani et al. (1995) and Trinidad et al. (2003). Most noticeably, the emission from source VLA 3, which is coincident with the position of the central source in the near-IR images, is better-determined; a close-up of the emission from this source is shown in panel b) of Fig. 5. The 3.6 cm emission from VLA 3 is consistent with a compact core plus jet morphology. The dominant side of the jet extends to the east, with a deconvolved width and length of  <0.2′′ and 1.2′′ (<670 AU and 4000 AU at 3.33 kpc), position angle of  ~100°, and an opening angle (at a radius of 4000 AU) of  <10° derived from its deconvolved width and length. The east jet is resolved, ending in a “knot” at 20h29m24.s\hbox{$\fs$}98 + 40°11′19.\hbox{$\farcs$}3 (J2000). There is also a smaller corresponding jet to the west, seen as a slight extension of the source in this direction. The morphology of VLA 3 is consistent with a jet-system which is approximately parallel to the larger-scale flow (observed by e.g. Hasegawa & Mitchell 1995), and therefore these observations confirm the proposal of Trinidad et al. (2003) that VLA 3 is its powering source.

Figure 6 shows the 7 mm multi-configuration image of the region surrounding AFGL 2591, in which only VLA 1 and VLA 3 were detected at this resolution. VLA 3 is compact, with a radius of  <500 AU, and is only slightly resolved in this image, showing no evidence of the eastern jet seen in Fig. 5 at 3.6 cm, which may partially be due to the 7 mm images having almost twice the rms noise. However, the A array-only image of van der Tak & Menten (2005), which has a synthesised beam size of 43 by 37 mas, shows that the source is slightly extended to the south-west, in a similar direction to the counter-jet seen in the 3.6 cm image. This extension is also seen in our multi-configuration 7 mm image. Using VLBA water maser observations, Sanna et al. (2012) found that this emission most likely arises from the outflow cavity walls, which the masers they observed were also found to trace.

Table 6 provides the measured positions, peak and integrated fluxes, as well as deconvolved sizes for the sources observed in the 3.6 cm (8.4 GHz) and 7 mm (43.3 GHz) images. These were measured using a custom-made irregular aperture photometry program, with which the integrated flux density was measured above 1 sigma. Uncertainties in the aperture fluxes were calculated to be a combination of the uncertainty due to the image noise over the aperture, and the maximum VLA absolute flux error, which is 2% and 5% for 3.6 cm and 7 mm respectively. Similarly, the peak flux uncertainty was found by combining the 1 sigma flux density with the VLA absolute flux error. The deconvolved source size at the 3 − σ level was measured by taking the major axis of the source to be along the direction which it is most extended, to within a position angle of 10°.

Table 6

Measured properties of the observed 3.6 cm (8.4 GHz) and 7 mm (43 GHz) continuum sources.

thumbnail Fig. 6

Map of the 7 mm continuum emission towards AFGL 2591. Contours are –4, 4, 7, 10, 15, 20, 25, 30  ×  rms noise = 56 μJy beam-1. Greyscale: –0.06 to 2.14 mJy beam-1 (1.2  ×  peak value). The inset panel shows a close-up of the 7 mm emission from VLA 3. The synthesised beam is shown in the bottom left corner of both images: 0.11′′    ×    0.11′′, PA = 43°.

thumbnail Fig. 7

Three-colour JHK′  Germini-North image of AFGL 2591 overlaid with the 3.6 cm contours from Fig. 5. Stretch of Gemini image: red, K′ band: 150−2500 MJy sr-1; green, H band: 150−500 MJy sr-1; blue, J band: 40−80 MJy sr-1.

Figure 7 compares the 3.6 cm emission from the region to the near-IR emission recorded in the Gemini North images. The five detected Hii regions line-up with several features in the near-IR image. Firstly, the peak of VLA 3 is coincident with that of the central source of AFGL 2591, at the apex of the one sided reflection nebula. VLA 1 is instead anti-correlated with the diffuse near-IR emission. This dent or cavity in the cloud can be seen more clearly in the K band bispectrum speckle interferometry image of Preibisch et al. (2003). VLA 2 also appears to be coincident with the close binary discovered by Preibisch et al. (2003), and VLA 4 with a source in the Gemini North image at 20h29m24.s\hbox{$\fs$}3 + 40°11′28′′ (J2000; however this source is too faint to be seen in Fig. 7). In addition, VLA 5 may be powered by the bright source at 20h29m23.s\hbox{$\fs$}96 + 40°11′09.\hbox{$\farcs$}25 (J2000; which is the same source as 2MASS 20292393 + 4011105, previously mentioned in Sect. 3.2). The sources VLA 1 and VLA 3 are discussed in detail below.

5.2.1. The central source of AFGL 2591, VLA 3

Figure 8 shows the SED in frequency space of the central source of AFGL 2591, VLA 3, in black error bars. The fluxes shown for VLA 3, which are taken from the literature and this work, are listed in Tables 3 (for λ < 1 mm) and 7 (for λ > 1 mm). At λ > 1 mm, fluxes have been listed from the observations with the largest available uv coverage.

A greybody was fitted to the fluxes between 100 and 850 μm, giving a temperature of 25  ± 7 K and a dust emissivity exponent β of 2.3    ±    0.6. For these four fitted fluxes, the image photometry apertures contained most if not all of the source flux. It can be seen in Fig. 8 that the fluxes at millimetre wavelengths (87−226 GHz) fall short of the flux expected from the fitted greybody. However, as these fluxes are taken from interferometric observations, it is likely that they only contribute a fraction of the actual flux. For instance, van der Tak & Menten (2005) note that their quoted 226 GHz continuum flux represents only 5% of the total flux, due to lack of shorter baselines. An arrow is shown on Fig. 8 representing this correction, which agrees well with the flux expected at 226 GHz from the fitted greybody.

Extrapolating this greybody to 7 mm, a flux of 2.3 mJy is expected, but the measured flux of VLA 3 at 7 mm is 2.9 mJy. However, it is likely that a portion of the 7 mm flux is due to ionized gas emission as well as dust emission, and a significant fraction of this combined emission is resolved out by the interferometer. To estimate the contribution from dust emission at 7 mm, a model image was made of the best-fitting envelope with disk model found from the SED and image modelling presented in Sect. 5.1. The 7 mm VLA A to D array observations were simulated from the model image using the CASA task simdata, then combined and imaged. The resultant 7 mm model image did not have any significant emission above 3σ = 0.17 mJy beam-1 over the 0.1 arcsec2 photometry aperture for VLA 3, giving an upper limit of 0.57 mJy integrated over the aperture. Hence, we estimate that  >2.3 mJy of the observed 7 mm emission is from ionized gas, and  <0.57 mJy is due to dust emission. This is in agreement with the results of Sanna et al. (2012) who found that the morphology of the 7 mm emission seen in previous VLA A array observations most likely arises from photoionization of the outflow cavity walls.

thumbnail Fig. 8

Radio-SEDs of AFGL 2591-VLA 3 and VLA 1. The SED of AFGL 2591-VLA 3 is show in black error bars, and the resolved fluxes of AFGL 2591-VLA 3 from Marengo et al. (2000) are shown as red squares. Green squares show the fluxes of VLA 1 – errors are not shown as they are smaller than the markers. The black and green solid lines show greybody fits to the SEDs of AFGL 2591-VLA 3 between 100 and 850 μm, and VLA 1 at mid-IR fluxes, respectively (see Sects. 5.2.1 and 5.2.2). The green dashed line shows the fit to the long-wavelength fluxes for VLA 1.

To measure the spectral index α (where Sν ∝ να) of VLA 3, which may give further insight into the nature of the emission, the 3.6 cm and 7 mm data were re-imaged using uv distances which were common to both datasets, from 4.5 to 1037 Kλ, and the fluxes remeasured. The integrated fluxes were found to be 1.43 ± 0.03 mJy for 3.6 cm and 3.3    ±    0.17 mJy for 7 mm (which increased compared to the original 7 mm flux due to a reduction in resolution, so that fainter emission was brought above the noise level). Therefore the spectral index between 3.6 cm and 7 mm was found to be 0.5 ± 0.02, where the quoted error in the spectral index is due solely to the photometric errors, similar to the spectral index expected for an ionized wind, 0.6. However, as this value is measured from only two fluxes, it is likely to be more uncertain than the photometric errors given. This spectral index is also an upper limit, as we estimate a fraction of the 7 mm flux is due to dust emission. By applying the common uv distance range above, we derive  <0.51 mJy for the contribution from dust emission, and therefore  >2.8 mJy from ionized gas emission. Using this lower limit for the ionized gas emission at 7 mm, we estimate the spectral index to be between 0.4 and 0.5.

Table 7

Overview of observed fluxes for VLA 1 to VLA 5.

thumbnail Fig. 9

Continuum emission towards AFGL 2591 at 106.7 and 110.5 GHz (~2.8 and 2.7 mm respectively), shown in both contours and greyscale. The rms noise in the images is σ = 2.2 and 2.1 mJy beam-1 respectively. The greyscales extend from −3    ×    σ to 65.5 and 54.0 mJy beam-1. Contours are at −3,3,4,5,7,10,15   and   20    ×    σ. The crosses show the positions of VLA 1 to 3. The synthesised beams are 4.9′′    ×    4.1′′ PA = −4.8° and 4.3′′    ×    3.5′′ PA = 95° respectively.

5.2.2. VLA 1

Figure 8 shows the SED of VLA 1 as green squares. This source was detected in the mid-IR by Marengo et al. (2000), as well as at several radio frequencies, for which the fluxes are listed in Table 7. The spectrum of VLA 1 appears to be flat at radio wavelengths, with a fitted spectral index α of 0.0 ± 0.03. The fitted power-law is shown as a dashed line. However, we first need to verify whether only a fraction of the flux of VLA 1 is being recovered at each wavelength. The flux of VLA 1 measured in our 3.6 cm A-D array image is not larger than the fluxes measured in the A array only images in previous works (A-D array: 80 mJy, this work; A array only: 82 mJy, 94 mJy, Tofani et al. (1995) and Trinidad et al. (2003) respectively). In fact, there is a small decrease in the flux of VLA 1 for the larger uv coverage data. Therefore it is likely that most of the flux from VLA 1 has been recovered by observations which are sensitive to scales up to 7′′, the largest observable angular scale of VLA A array observations at 3.6 cm. This is the case for all of the interferometric fluxes shown for VLA 1, hence the calculated flat spectral index of VLA 1 is not likely to be due to instrumental effects. The observed spectral index is close to that expected from optically thin free-free ionized gas emission (α = −0.1).

The mid-IR emission from VLA 1 is likely to arise from thermal emission from dust within the Hii region. To estimate how much VLA 1 contributes to the total flux of the region, a greybody spectrum with a dust opacity index β = 2 was fit to the fluxes measured by Marengo et al. (2000). The best-fitting greybody has a temperature of 111    ±    14 K and peaks at 26    ±    3.3 μm with a flux of 112 Jy, which is uncertain by a factor of two. Therefore the thermal dust emission from VLA 1 never contributes more than  ~30% to the total flux from the region at wavelengths between 8 and 100 μm and a negligible fraction of the observed flux (<0.5%) at other wavelengths, and therefore should not significantly affect the SED and image modelling results given in Sect. 5.1. The mid-IR fluxes for only the resolved central source of AFGL 2591 from Marengo et al. (2000) are plotted on Fig. 8 (red squares) and agree well with the ISO-SWS data, which were used to cover the mid-IR to far-IR section of the SED.

5.3. 13CO, C18O and 3 mm continuum

This section details the properties of the dense molecular gas surrounding AFGL 2591 inferred from the C18O molecular line and millimetre continuum observations.

Figure 9 shows the continuum emission from the region at 106.7 and 110.5 GHz, or 2.81 and 2.7 mm, in both of the wide  ~3 mm continuum bands observed. The morphology of the emission is similar to that seen in the millimetre continuum maps at various wavelengths presented by van der Tak et al. (1999), van der Tak et al. (2006), and resolved further by Wang et al. (2012). We found that the sources in the 3 mm images were systematically offset from the 3.6 cm peak positions, which we determined to be caused by inaccurate coordinates used for the phase calibrator MWC 349 in the mm observations. We therefore shifted the continuum and line maps by 0.89′′ in RA and 0.18′′ in declination to agree with the 7 mm position of MWC 349 A reported in Rodríguez et al. (2007). In Fig. 9, the peak positions of VLA 1 to 3 are shown as crosses. The main two sources that are detected are VLA 1 and VLA 3, but there is also extended emission towards the position of VLA 2. Although the large uncertainties in the two measured 3 mm fluxes do not allow calculation of accurate spectral indices for VLA 1 and 2 using these values alone, the difference in the spectral index of the two sources can be seen from the difference in their relative brightnesses between the two images, indicating VLA 1 has a flat spectral index, and VLA 3 has a rising spectral index at shorter wavelengths. In addition, when we compared both images convolved with the beam of the 106.7 GHz image, the morphology of the emission was very similar, apart from the clear difference in the flux of VLA 3. The sources were fit with 2D Gaussians using the CASA task imfit to find the fluxes, which are given in Table 7.

thumbnail Fig. 10

C18O (black line) and 13CO (grey line) spectral profiles measured at the position of AFGL 2591-VLA 3: 20h29m24.s\hbox{$\fs$}88 + 40°11′19.\hbox{$\farcs$}5 (J2000).

Figure 10 presents C18O and 13CO line profiles at the position of AFGL 2591-VLA 3. The 13CO observations were missing significant flux on extended scales which, in combination with self-absorption effects, is likely to be the cause of the dip in the 13CO line profile seen in Fig. 10. As there was a significant amount of missing flux in the 13CO line, and the emission was distributed incoherently across the maps, it was not possible to interpret this emission. Therefore the 13CO data will not be discussed further, however channel maps of the 13CO emission are included in Figs. 11 and 12.

Figure 13 presents channel maps of the C18O emission towards AFGL 2591 between –9.0 and –2.7 km s-1 (where –5.7 km s-1 is the rest velocity of the cloud), and Fig. 14 shows the C18O intensity-weighted first moment map. In both figures the positions of VLA 1 through 3 are marked by crosses, however in Fig. 14, the position and size of VLA 1 at 3.6 cm is shown instead as a circle. A blue-shifted velocity feature can be seen in Fig. 14 to the west of VLA 3, and coincident with VLA 2, which decreases smoothly in velocity to the west VLA 3, ranging from –5.3 to –8.0 km s-1. This can be seen in the C18O channel maps (Fig. 13), which also show that the shape of the emission becomes narrower away from the line centre and towards the west. This is most obvious in the –7.7 and –8.0 km s-1 channels. In addition, Fig. 14 shows that there is a smaller intensity peak  ~6′′ to the south of VLA 3. Here, the velocity instead decreases from northwest to southeast from approximately –5 to –6 km s-1.

In Fig. 15, integrated C18O emission in several velocity ranges is compared with the Near-IR Gemini image, and the positions of bow shocks detected in the K′ band Gemini North image (first noted by Preibisch et al. 2003). The bow shocks are shown approximately as white semi-ellipses, as they are too faint to be seen directly from Fig. 15. The white dashed line also shows the direction of the ionized jet seen in the 3.6 cm images (Fig. 5), with PA 100°. The contours show the high-velocity red and blue-shifted gas (red and blue contours respectively). Although the “high”-velocity channels referred to in this work are high velocity relative to the observed line, it should be noted that more extended, higher velocity gas has been detected by Hasegawa & Mitchell (1995) using single dish 12CO observations with 14.3′′ resolution, which extends in velocity from –45 to 35 km s-1. Therefore, the observations presented here trace the higher density but comparatively lower-velocity gas, within a few km s-1 of the line centre. To minimise confusion due to overlap in velocity of various components, the contours in Fig. 15 only show integrated C18O emission from –4.3 to –3.3 km s-1 (red) and –8.0 to –7.0 km s-1 (blue).

6. Discussion

6.1. The jet and outflow from VLA 3

thumbnail Fig. 13

C18O channel map at 0.3 km s-1 resolution between –9.0 and –2.7 km s-1. The rms in the map is σ = 0.1 Jy beam-1. The map peak flux is 1.1 Jy beam-1. Contours are at –3, 3, 4, 5, 6, 7, 8, 9, 10, 11    ×    σ. The synthesised beam is shown in the bottom left-hand corner (4.5′′    ×     3.6′′, PA 93°). A scale size of 50 000 AU is represented by a bar in the bottom right-hand panel. The cross shows the position of the central source of AFGL 2591 from the Gemini near-IR J band image. The rest velocity of the cloud is –5.7 km s-1.

In Fig. 15, the extrapolated direction of the ionized jet observed at 3.6 cm points directly through the red-shifted C18O contours to the bow-shocks seen in the K′ band image. Hence this outlines a coherent picture in which the red-shifted outflow lobe of AFGL 2591 consists of a small-scale 4000 AU ionized jet that is part of a jet or wind extending out to  ~0.4 pc where the flow terminates against the surrounding cloud as bow-shocks. However it is interesting to note that the position angle of the jet and bow shocks do not exactly align with that of the blue-shifted reflection nebula or the elongation of the blue-shifted outflow observed by Hasegawa & Mitchell (1995). Therefore other stars in the vicinity, such as the powering star(s) of VLA 1, may be causing precession of the jet. If, as well as the emission lying along the position angle of the ionized jet, the clump of emission at 20h29m26.s\hbox{$\fs$}7 + 40°11′23′′ (J2000) is also included, these three clumps may be tracing an arc of emission describing the densest parts of the red-shifted outflow lobe, created as the jet precesses. If so, the extent of these clumps corresponds to an observed opening angle of  ~40° at their distance from the source,  ~0.4 pc. At the same distance, the observed blue-shifted outflow lobe is larger with an opening angle of  ~60°. The opening angles of the biconical outflow observed by Jiménez-Serra et al. (2012) and the blue shifted outflow cone observed by Sanna et al. (2012), on scales  <2′′ or  <6700 AU, are  ~90–110°. The larger observed opening angles on smaller scales suggests that the outflow cavity can not be described by a cone with a constant opening angle, but is instead better described by a power-law cavity similar to that used in our radiative transfer models.

Several suggestions regarding the nature of the mm and cm continuum emission from VLA 3 have been made by previous studies, including a core-halo Hii region, an ionized wind with a dust disk (Trinidad et al. 2003), or emission from a spherical gravitationally confined Hii region (van der Tak & Menten 2005). At 3.6 cm, the deeper observations presented here clearly show that as well as a central compact core, the source also exhibits a non-spherical jet-like morphology. Therefore, we calculate several properties of the emission below, assuming it originates from a jet.

Without assuming a specific ionization mechanism, the mass loss rate of the jet can be estimated using the model of Reynolds (1986), which describes the emission from a partially optically thick ionized jet (his Eq. (19)): 10-6M yr-1=9.38×10-2(υ100kms-1)(1x0)(μmp)×  [(SνmJy)(ν10GHz)α]0.75(dkpc)1.5(νm10GHz)0.45+0.75α×  θ0.75(T104K)-0.075(sini)-0.25 F-0.75\begin{eqnarray} &&\frac{\dot{M}}{10^{-6}\,{M}_{\odot}~\rm{yr}^{-1}} = 9.38\times10^{-2} \left( \frac{\upsilon}{100\,\rm{km\,s}^{-1}} \right) \left( \frac{1}{x_0} \right) \left( \frac{\mu}{m_p} \right)\notag \\ &&\qquad\times ~~ \left[ \left( \frac{S_{\nu}}{\rm{mJy}} \right) \left( \frac{\nu}{10\,\rm{GHz}} \right)^{-\alpha} \right]^{0.75} \left( \frac{d}{\rm kpc} \right)^{1.5} \left( \frac{\nu_{m}}{10\,\rm{GHz}} \right)^{-0.45+0.75\alpha} \notag \\ &&\qquad\times ~~ \theta^{0.75} \left(\frac{T}{10^4\,\rm{K}}\right)^{-0.075} (\sin{i})^{-0.25}~F^{-0.75} \end{eqnarray}(8)\arraycolsep1.75ptwhere υ is the velocity of the jet, x0 is the ionization fraction; μ/mp is the mean particle mass per hydrogen atom of the ionized material, given by 1/(1 + x0); Sν is the observed flux at the frequency ν; α is the spectral index; d is the distance to the source; νm is the turn-over frequency below which the emission becomes optically thick; θ is the opening angle of the flow, defined as the ratio of the projected width to the radius at the base of the jet; T is the temperature of the ionized gas; i is the inclination measured from the line-of-sight and F is a function of the spectral index and the dependance of optical depth on radius (see Reynolds 1986, Eq. (17)).

Assuming an isothermal, uniformly ionized jet with a density gradient of ρ ~ r-2; υ = 500 km s-1 for the velocity of the jet, found for the blue-shifted HH objects towards AFGL 2591 (Poetzel et al. 1992); x0 = 0.1 (commonly found for low-mass sources, e.g. Bacciotti & Eislöffel 1999); vm = 50 GHz; θ ~  0.2′′/1.2′′, the ratio of the maximum deconvolved width to the maximum radius; T = 104 K ; a spectral index between 0.4 and 0.5; and i = 50°, taken from the results of the SED and near-IR image profile modelling presented in Sect. 5.1; F was found to be between 2.3 and 1.8, and the mass loss rate of the jet observed at 3.6 cm was therefore determined to be in the range 0.77−1.0    ×    10-5   M   yr-1.

This value is approximately a thousand times larger than the jet mass loss rates seen for low-mass protostars, which are commonly found to be on the order of 10-8   M   yr-1 (e.g. Podio et al. 2006). The mass loss rate in this jet can also be compared to the mass loss rate of the small-scale red-shifted molecular flow observed by Hasegawa & Mitchell (1995) in 12CO(J = 3−2), which corresponds well to the position and direction of the outflow lobe suggested by the ionized jet, red-shifted C18O emission and bow shocks discussed in Sect. 5.3. For their optically thick case, Hasegawa & Mitchell find a mass loss rate of 6.3    ×    10-5   M yr-1 at 1 kpc corresponding to 7.0    ×    10-4   M yr-1 at 3.33 kpc.

Multiplying the ionized jet mass loss rate by the velocity υ, the momentum rate in the jet is 3.9−5.2 × 10-3   M yr-1 km s-1, which is very similar to the momentum rates determined by Hasegawa & Mitchell (1995) for both the large and small scale red-shifted 12CO outflows (7.7 and 7.4    ×    10-3   M yr-1 km s-1 at 3.33 kpc with i = 50°, scaled from 7.5 and 7.2    ×    10-4   M yr-1 km s-1 at 1 kpc with i = 45°). Thus the jet itself would have very close to the required momentum to drive the observed larger scale red-shifted outflow.

The ionization mechanism of the jet can be modelled as a plane-parallel shock in a homogeneous neutral “stellar” wind (e.g. Curiel et al. 1989; Anglada 1996). The momentum rate of the jet can be expressed as: (Myr-1kms-1)=3.13×10-4η(Sνd2mJykpc2)(υ200kms-1)0.32×(T104K)-0.45(ν5GHz)0.1(τ1eτ)\begin{eqnarray} &&\left( \frac{\dot{P}}{{M}_{\odot}\,\rm{yr}^{-1}\,\rm{km\,s}^{-1}} \right) = \frac{3.13\times10^{-4}}{\eta} \left( \frac{S_{\nu}d^{2}}{\rm{mJy}\,\rm{\rm kpc}^2} \right) \left( \frac{\upsilon_{\star}}{200\,\rm{km\,s}^{-1}} \right)^{0.32} \notag\\ &&\qquad\times \left( \frac{T}{10^4\,\rm{K}} \right)^{-0.45} \left( \frac{\nu}{5\,{\rm GHz}} \right)^{0.1} \left( \frac{\tau}{1-{\rm e}^{-\tau}}\right) \arraycolsep0pt \end{eqnarray}(9)where η is the shock efficiency or ionization fraction, found to be  ~0.1 for low-mass sources (e.g. Bacciotti & Eislöffel 1999), υ is the initial velocity of the stellar wind or jet, taken to be 500 km s-1, and τ is the optical depth of the emitting gas. Here, the flux Sν is measured at 3.6 cm. The optical depth can be determined using Eq. (6) of Anglada et al. (1998): α=2+ln[(1eτ1)/(1eτ1(ν1/ν2)2.1)]ln(ν1/ν2)·\begin{equation} \alpha = 2 + \frac{\ln{\left[(1-{\rm e}^{-\tau_1})/\left(1-{\rm e}^{-\tau_1(\nu_1/\nu_2)^{2.1}}\right)\right]}}{\ln{(\nu_1/\nu_2)}}\cdot \end{equation}(10)Using the estimated spectral index of 0.4−0.5, the optical depth of the emission at 3.6 cm was estimated to be between 2.1 and 2.6, giving a value for the required momentum rate in the jet of 2.8−3.4    ×    10-2   M yr-1 km s-1. Therefore the required momentum rate in the jet is larger than the value of 3.9−5.2    ×    10-3   M yr-1 km s-1 calculated above by a factor of  ~5 to 9, implying that the 3.6 cm emission cannot be solely caused by shocks in the jet. However, if the shock efficiency is increased to near unity, the momentum rate of the jet would be sufficient to ionize it. As this is unlikely, photoionization by the central star therefore must provide a significant fraction of the emission at 3.6 cm.

thumbnail Fig. 14

Shown in colour-scale: intensity-weighted first moment map of the C18O emission towards AFGL 2591. The corresponding velocities are shown on the colour-bar on the right. The black open circle represents the position and size of VLA 1 at 3.6 cm. The crosses indicate the position of VLA 2 and VLA 3. A scale size of 20 000 AU is represented by a bar in the upper left-hand corner. Shown in contours: the C18O integrated intensity map from –2.3 to –8.7 km s-1. The peak flux of the integrated intensity map is 2.5 Jy beam-1 km s-1. Contours are at –3, 3, 4, 5, 6, 7    ×    σ = 0.33 Jy beam-1 km s-1. The synthesised beam is shown in the lower right-hand corner (4.5′′ ×  3.6′′, PA 93°). The rest velocity of the cloud is –5.7 km s-1.

Turning to the C18O observations, the blue-shifted emission seen in Fig.14 and the blue contours in Fig. 15 may be tracing a blue-shifted outflow from VLA 2 or VLA 3. If this is the case, the outflow can be seen to be narrower and more collimated at more negative velocities, seen in Fig. 13. In addition, a “Hubble-like” velocity trend (i.e. an systematic increase flow velocity with distance from the source, e.g. Lada & Fich 1996) can be seen across the blue-shifted emission in Fig. 14. Therefore there is a high-velocity collimated component at the centre and front of this blue-shifted flow, which smoothly transitions into a surrounding lower-velocity wide-angle wind.

Although the integrated C18O emission peaks on VLA 2, thus making it a possibility that VLA 2 is driving the outflow, there are several lines of evidence which suggest this is not the case. Firstly, the Hii region surrounding VLA 2 appears to be more evolved than VLA 3, and hence VLA 2 is less likely to be driving an outflow. Secondly, the emission at the position of VLA 3 is at the rest velocity of the cloud (5.7 km s-1) instead of VLA 2, where it is approximately 6.3 km s-1. Thirdly, the blue-shifted emission seen in C18O agrees very well with the high-velocity large-scale emission seen in Fig. 3 of Hasegawa & Mitchell (1995). Therefore, if it were the case that the blue-shifted C18O emission was created by VLA 2, this source would also be responsible for the large-scale blue-shifted outflow seen by Hasegawa & Mitchell (1995). However, because this source has a much lower total luminosity compared to VLA 3, which dominates the emission at all wavelengths, and as this would leave VLA 3 without a corresponding blue-shifted outflow to the red-shifted one observed, we consider it more likely the blue-shifted emission is the inner part of the large-scale blue-shifted outflow of AFGL 2591-VLA 3.

Assuming that the observed blue-shifted C18O emission is the inner part of the entrained large-scale blue-shifted outflow of AFGL 2591-VLA 3, the outflow properties were calculated, and are given in Table 8. The blue-shifted outflow mass was found assuming the gas was optically thin (derived in Appendix A),

thumbnail Fig. 15

Gemini North three-colour JHK′ image of AFGL 2591 overlaid with contours of C18O emission integrated over –4.3 to −3.3 km s-1 (red) and –8.0 to –7.0 km s-1 (blue). The three arcs in the east of the image show the positions of bow shocks visible in the Gemini North image, which are too faint to be seen directly in the figure. The dashed line shows the direction of the ionized jet observed at 3.6 cm (PA 100°), and yellow diamonds show roughly the peak positions of the H2 knots observed in Tamura & Yamashita (1992). Red-shifted contours are –3, 3, 4, 5, 6  ×  0.1 Jy beam-1 km s-1 and blue-shifted contours are –3, 3, 4, 5, 6, 7  ×  0.1 Jy beam-1 km s-1. The synthesised beam is shown in the lower left-hand corner (4.5′′ ×  3.6′′, PA 93°). Near-IR image stretch as in Fig. 7.

M outflow = 2.75 × 10 -5   [ 12 CO C 18 O ] ( T + 0.882 ) e 5.27 / T   d 2   S ν   d υ   ( M ) \begin{equation} M_{\rm{outflow}} = 2.75 \times 10^{-5} ~ \left[ \frac{^{12}\rm{CO}}{\rm{C}^{18}\rm{O}}\right] \frac{(T+0.882)}{{\rm e}^{-5.27/T}}~d^2~ \int S_{\nu} ~{\rm d}\upsilon~({M}_{\odot}) \label{massinflow} \end{equation}(11)

where 12COC18O][\hbox{$\left[ \frac{^{12}\rm{CO}}{\rm{C}^{18}\rm{O}}\right]$} is the abundance ratio of 12CO to C18O, T is the temperature of the gas assuming local thermodynamic equilibrium, dkpc is the distance to the source in kpc, and Sν dυ\hbox{$\int S_{\nu} ~{\rm d}\upsilon$} is the integral of the flux over the velocities of the blue-shifted emission. The flux in each channel was found by integrating the emission within an irregular aperture around only the blue-shifted outflow emission (i.e. not including emission from the resolved peak to the south of the central position of AFGL 2591). The temperature was assumed to be 25 K from the greybody fit to the far-IR and sub-mm fluxes in Sect. 5.2.1, and the abundance ratio [12CO/C18O]\hbox{$\left[{^{12}\rm{CO}}/{\rm{C}^{18}\rm{O}}\right]$} was calculated using the results given in Wilson & Rood (1994), where DGC\hbox{$D_{_{\rm GC}}$} is the galactocentric distance in kpc,

[ 12 CO C 18 O ] = ( 58.8 ± 11.8 ) D GC + ( 37.1 ± 82.6 ) , \begin{equation} \left[ \frac{^{12}\rm{CO}}{\rm{C}^{18}\rm{O}}\right]= (58.8 \pm 11.8) D_{_{\rm GC}} +(37.1 \pm 82.6), \end{equation}(12)

giving a value of 12COC18O][=538\hbox{$\left[ \frac{^{12}\rm{CO}}{\rm{C}^{18}\rm{O}}\right]=538$} for AFGL 2591 for a galactocentric distance of 8.5 kpc.

It should be noted that the shape of the C18O line profile shown in Fig. 10 may indicate that this line is partially optically thick, so that the estimates of the outflow properties depending on the derived mass in Table 8 will be lower limits. In addition, as the C18O emission only traces the densest inner parts of the flow, it only represents a small fraction of the mass and momentum from the entire large-scale outflow.

To find the momentum MoutflowV and kinetic energy 12MoutflowV2\hbox{$\frac{1}{2}M_{\rm{outflow}}V^2$} in the flow, the final term of Eq. (11) was replaced by Sν(υυLSR)/sini dυ\hbox{$\int S_{\nu} (\upsilon - \upsilon_{\rm LSR})/\sin{i} ~{\rm d}\upsilon$} or 12Sν((υυLSR)/sini)2 dυ\hbox{$\frac{1}{2}\int S_{\nu} ((\upsilon - \upsilon_{_{\rm LSR}})/\sin{i})^2 ~{\rm d}\upsilon$} respectively, where υ is the velocity of each channel, υLSR = −5.7 km s-1 is the rest velocity of the cloud, and i is the inclination of the flow relative to the line of sight, taken to be 50° from the results of the SED modelling presented in Sect. 5.1.

Table 8

Derived properties of the inner part of the blue-shifted outflow traced by the C18O emission.

The mass weighted average velocity of the flow was calculated as the momentum over the mass in the flow:  ⟨ V ⟩  = MoutflowV/Moutflow, and the projected length of the flow was measured as the distance from the central source of AFGL 2591 to the furthest 3σ contour from the source in the blue-shifted emission shown in Fig. 15, which was found to be  ~14.2′′ or 47 300 AU at 3.33 kpc. Therefore the actual length of the flow l was calculated as l = lproj/sini. The dynamical timescale of the flow was determined using the equation tdyn = l/ ⟨ V ⟩ , however it should be noted that, as with any outflow, this value may not be an exact measure of the time taken to displace the gas, as the gas may have been set in motion at some distance from the central source. In addition, this assumes a constant outflow velocity, but the gas is observed to be accelerating away from the source. Nevertheless, this provides a characteristic timescale, which can be compared to other studies. Comparing to dynamical timescales found for large scale flows toward forming massive stars (of order 104−105 yr for L > 104   L, Beuther et al. 2002), the timescale for the inner part of the flow seen in C18O is at the upper end of these values. Scaling the blue-shifted dynamical timescale of 6.3 × 104 yr for the 12CO outflow observed by Hasegawa & Mitchell (1995) to the new distance of 3.33 kpc, we obtain a dynamical timescale of 2.1 × 105 yr, similar to the value we find for the smaller-scale C18O outflow.

6.2. Cluster properties

The section below deals with the overall properties of the forming cluster surrounding AFGL 2591-VLA 3. The maximum envelope radius for AFGL 2591-VLA 3, found by our modelling in Sect. 5.1, is on parsec scales (180 000 AU and 420 000 AU or 0.87 and 2.0 pc for the envelope with and without disk models respectively). This is also corroborated by the large-scale envelope observed in the 12CO and 13CO maps of van der Wiel et al. (2011), as well as in the 80′′ or 270 000 AU diameter disk-like structure seen in the CS(1−0) maps of Yamashita et al. (1987). This coherent parsec-scale structure, containing at least five forming stars (VLA 1 through 5), therefore constitutes a cluster-forming clump.

The spectral types of the stars required to power VLA 1, 2, 4 and 5 were calculated using the 3.6 cm fluxes, by the same method as outlined in Johnston et al. (2009), assuming they are optically thin Hii regions with one powering star, giving spectral types of B0.5, B1, B2 and B1 at 3.33 kpc for VLA 1, 2, 4 and 5 respectively. As VLA 3 has a spectral index of between 0.4 and 0.5, its emission is not optically thin, and thus such an estimate would not be accurate. As the morphology is consistent with a compact core plus jet, it may be that some of its emission is not from an HC Hii region, but is instead created by shock-ionized gas. However, assuming photoionization dominates the emission from VLA 3, the spectral type of the powering star of VLA 3 was estimated to be O8–O9 under the assumption of a photoionized wind, using the mass loss rate of the jet found in Sect. 6.1: 0.77−1.0    ×    10-5   M   yr-1, the assumed wind velocity: 500 km s-1, and Table 4 of Rodriguez & Canto (1983). For the above spectral type estimates, we note that the fluxes in Table 6 may only represent a fraction of the total flux at these wavelengths, due to the finite range in spatial scales which the observations probe, and therefore they provide lower limits.

Following Sanna et al. (2012), we take the mass of the most luminous star VLA 3 to be 38 M. Assuming sources VLA 1 to 5 are cluster members, the lowest luminosity powering star (VLA 4 which is B2) corresponds to a ZAMS mass of  ~9 M, found using the ZAMS luminosities of Panagia (1973) and the theoretical HR diagram of Maeder & Meynet (1989). Thus if this cluster has a fully populated IMF, extrapolating a Krupa IMF (Kroupa 2001) down to 1 M suggests there should be  ~100 stars in this cluster with mass above 1 M.

The roughly circular shape of VLA 1 (with a radius of  ~2′′ or  ~6700 AU at 3.33 kpc) and its optically thin spectral index, supports the view that it is a UC Hii region excited by one or more massive stars. However, an exciting star for VLA 1 is not seen in either the 2MASS or Gemini North near-IR images. Assuming a single star, the ZAMS B0.5 absolute magnitude given in Panagia (1973) of −2.8, and a distance of 3.33 kpc, the spectral type of the calculated exciting star of VLA 1 corresponds to a 2MASS Ks band magnitude of 11.1. The detection limit of the 2MASS point source catalog in Ks band is 14.3 mag (Skrutskie et al. 2006), which indicates that there must be more than 3.2 mag of extinction at Ks band for this source to be obscured, corresponding to 26 mag of extinction in the V band (Cardelli et al. 1989). These high extinctions, in addition to the indentation that VLA 1 has made in the ambient cloud traced by near-IR emission, suggest that VLA 1 lies within the same cloud core as AFGL 2591-VLA 3.

The dent in the cloud at the position of VLA 1 provides support to the idea that this UC Hii region is at the same distance, and is expanding into the material illuminated by the central source of AFGL 2591. Therefore it seems likely that the object at approximately 20h29m24.s\hbox{$\fs$}6 + 40°11′16.\hbox{$\farcs$}4 (J2000), suggested to be an edge-on disk by Preibisch et al. (2003), is instead the brightened rim of the expanding UC Hii region, where it is interacting with the surrounding cloud. Hence, through its expansion, VLA 1 may have been able to affect the formation of other members of the cluster.

The secondary peak in Fig. 14 shows a velocity gradient across it in roughly the NW-SE direction. There are no sources at this position in any of the images, limiting the possibility that this emission is from an outflow of a nearby object. Instead the velocity of the emission appears to decrease radially away from VLA 1, which is marked as a circle in Fig. 14. Indeed, Fig. 7 indicates that VLA 1 is preferentially expanding towards the southeast, into possibly lower density gas. Hence the origin of this velocity gradient in Fig. 14 may be due to the expansion of the gas around VLA 1 towards the observer, with a projected velocity of  ~1 km s-1 relative to the ambient cloud. Further, it is interesting to note that the lack of C18O emission at the position of VLA 1 again suggests that VLA 1 is at the same distance as AFGL 2591-VLA 3. The lack of C18O emission towards VLA 1 could be due to the fact it has swept-up the molecular cloud material around it, or the molecular gas has been dissociated by the high temperatures (~104 K) in the Hii region.

7. Conclusions

We have carried out a study of the emission from the luminous star-forming region AFGL 2591 from near-IR through to cm wavelengths, to determine whether the properties of the source can be explained by the same processes thought to govern the formation of low-mass stars.

In Sect. 4, we described a Monte Carlo dust continuum radiative transfer dust code to model the continuum absorption, emission and scattering through two scaled-up azimuthally symmetric dust geometries, the first consisting of a rotationally flattened envelope with outflow cavities, and the second also including a flared disk. In Sect. 5.1 our results show that in both cases, these models were able to reproduce the main features of the SED and 2MASS images of AFGL 2591.

The observed 3.6 cm images presented in Sect. 5.2 detect and more finely resolve the central powering source AFGL 2591-VLA 3, which has a compact core plus collimated jet morphology. The east component of the jet extends 4000 AU from the central source, with an opening angle of <10° at this radius. This jet is also roughly parallel to the outflow seen at larger scales by Hasegawa & Mitchell (1995), confirming that VLA 3 is the powering source of the outflow from AFGL 2591. The multi-configuration 7 mm image detects both VLA 1 and VLA 3. At this wavelength, VLA 3 does not show a jet morphology, but instead presents compact (<500 AU) emission, some of which (<0.57 mJy of 2.9 mJy) is estimated to be due to dust emission. Comparison with ancillary 0.3−0.4′′ Gemini near-IR images show that VLA 3 is coincident with the illuminating source of AFGL 2591, at the apex of the observed near-IR reflection nebula. Further, VLA 2 is coincident with the close binary first seen by Preibisch et al. (2003), and VLA 4 and a newly detected source VLA 5, are coincident with near-IR sources.

By fitting greybodies to the SEDs of VLA 1 and VLA 3 at wavelengths shorter than 1 mm, the contribution of VLA 1 to the total SED from AFGL 2591 was determined to be  <30% at wavelengths between 8 and 100 μm, and to be negligible at other wavelengths. Therefore contamination by VLA 1, or fainter sources in the region (e.g. VLA 2), should not significantly change the emergent SED of the central powering object which dominates the luminosity from the region, AFGL 2591-VLA 3.

The nature of the emission from VLA3 was also investigated. The spectral index of VLA 3 was calculated (between 3.6 cm and 7 mm) to be between 0.4 and 0.5, similar to that of an ionized wind. In Sect. 6.1, the 3.6 cm emission was modelled as an ionized jet, from which the mass loss rate in the jet was determined to be in the range 0.77−1.0    ×    10-5   M yr-1. In addition, it was found that the jet may have enough momentum to drive the larger-scale flow. However, assuming a shock efficiency of 10%, the momentum rate of the jet was found to not be sufficient to ionize it via only shocks. Thus a significant portion of the 3.6 cm emission is due to photoionization.

In Sect. 5.3, the observed C18O emission towards the source was found to be tracing the densest accelerated or entrained material. The blue-shifted emission, to the west of VLA 3 and coincident with VLA 2, shows an increasing velocity gradient to the west, which is more collimated at higher velocities, and is coincident with the inner part of the blue-shifted outflow observed by Hasegawa & Mitchell (1995). Thus although we cannot rule out that the that this emission is an outflow launched by VLA 2, there is substantial evidence that the C18O emission is tracing the densest parts of the entrained gas in the large-scale outflow from VLA 3. In addition, the dense red-shifted gas traced by C18O lies along the position angle of the ionized jet seen at 3.6 cm, directly between the jet and bow shocks seen in the near-IR Gemini image. Therefore it is likely this emission also traces a small, dense, segment of the red-shifted flow from VLA 3.

The above results show that many of the properties of AFGL 2591-VLA 3 resemble those of low-mass protostars. For example, the geometry of AFGL 2591 determined from modelling the SED and near-IR images of the source is consistent with that of a rotationally flattened envelope with and without a flared disk. Furthermore, the outflow of AFGL 2591-VLA 3 is comprised of both a collimated ionized jet, and a wide angle wind, a property seen in the outflow phenomena of many low-mass sources (e.g. Arce et al. 2007). However, although the emission at 3.6 cm from AFGL 2591-VLA 3 has the morphology of a jet, it is unlikely that it can be explained solely by shocks in a neutral wind or outflow. Therefore some part of the compact emission from the star may instead be provided by a hypercompact Hii region. Yet if this is the case, its presence has not disrupted the accretion or outflow processes of AFGL 2591-VLA 3. Thus, in this manner, the formation of the central dominant object VLA 3 of AFGL 2591 does not appear to be significantly different to that of low-mass protostars. However, it is important to note that this star is not forming in isolation, evidenced by the four other cm sources observed within the bounds of the envelope determined for AFGL 2591-VLA 3, with a radius found from our modelling to be on parsec scales. Therefore another way to view this picture is that AFGL 2591-VLA 3 is able to source its accreting material from the shared gas reservoir of a small cluster while still exhibiting the phenomena expected during the formation of low-mass stars.

Online material

thumbnail Fig. 11

13CO channel map at 0.3 km s-1 resolution between –12.0 and –5.7 km s-1. The rms in the map is σ = 0.1 Jy beam-1. The map peak flux is 1.3 Jy beam-1. Contours are at –3, 3, 5, 10, 15     ×    σ. The synthesised beam is shown in the bottom left-hand corner (4.4′′    ×     3.7′′, PA 96°). A scale size of 0.5 pc is represented by a bar in the bottom right-hand panel. The cross shows the position of the central source of AFGL 2591 from the Gemini near-IR J band image. The rest velocity of the cloud is –5.7 km s-1.

thumbnail Fig. 12

13CO channel map at 0.3 km s-1 resolution between –5.7 and 0.6 km s-1. The rms in the map is σ = 0.1 Jy beam-1. The map peak flux is 1.3 Jy beam-1. Contours are at –3, 3, 5, 10, 15  ×    σ. The synthesised beam is shown in the bottom left-hand corner (4.4′′    ×    3.7′′, PA 96°). A scale size of 0.5 pc is represented by a bar in the bottom right-hand panel. The cross shows the position of the central source of AFGL 2591 from the Gemini near-IR J band image. The rest velocity of the cloud is –5.7 km s-1.


1

Now the Jansky Very Large Array.

2

Support for CARMA construction was derived from the states of California, Illinois, and Maryland, the James S. McDonnell Foundation, the Gordon and Betty Moore Foundation, the Kenneth T. and Eileen L. Norris Foundation, the University of Chicago, the Associates of the California Institute of Technology, and the National Science Foundation. Ongoing CARMA development and operations are supported by the National Science Foundation under a cooperative agreement, and by the CARMA partner universities.

3

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Acknowledgments

We would like to thank the referee for providing insightful comments which improved this work. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work makes use of observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the Science and Technology Facilities Council (United Kingdom), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), Ministério da Ciência e Tecnologia (Brazil) and Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina). We also used APLpy, an open-source plotting package for Python hosted at http://aplpy.github.com.

References

  1. Anglada, G. 1996, in Radio Emission from the Stars and the Sun, eds. A. R. Taylor, & J. M. Paredes, ASP Conf. Ser., 93, 3 [Google Scholar]
  2. Anglada, G., Villuendas, E., Estalella, R., et al. 1998, AJ, 116, 2953 [NASA ADS] [CrossRef] [Google Scholar]
  3. Arce, H. G., Shepherd, D., Gueth, F., et al. 2007, Protostars and Planets V, 245 [Google Scholar]
  4. Bacciotti, F., & Eislöffel, J. 1999, A&A, 342, 717 [NASA ADS] [Google Scholar]
  5. Beuther, H., Schilke, P., Sridharan, T. K., et al. 2002, A&A, 383, 892 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Bonnell, I. A., Smith, R. J., Clark, P. C., & Bate, M. R. 2011, MNRAS, 410, 2339 [NASA ADS] [CrossRef] [Google Scholar]
  7. Campbell, B. 1984a, ApJ, 287, 334 [NASA ADS] [CrossRef] [Google Scholar]
  8. Campbell, B. 1984b, ApJ, 282, L27 [NASA ADS] [CrossRef] [Google Scholar]
  9. Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245 [NASA ADS] [CrossRef] [Google Scholar]
  10. Castelli, F., & Kurucz, R. L. 2004, [arXiv:astro-ph/0405087] [Google Scholar]
  11. Commerçon, B., Hennebelle, P., & Henning, T. 2011, ApJ, 742, L9 [NASA ADS] [CrossRef] [Google Scholar]
  12. Curiel, S., Rodriguez, L. F., Bohigas, J., et al. 1989, Astrophys. Lett. Commun., 27, 299 [NASA ADS] [Google Scholar]
  13. D’Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411 [NASA ADS] [CrossRef] [Google Scholar]
  14. de Wit, W. J., Hoare, M. G., Fujiyoshi, T., et al. 2009, A&A, 494, 157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Di Francesco, J., Johnstone, D., Kirk, H., MacKenzie, T., & Ledwosinska, E. 2008, ApJS, 175, 277 [NASA ADS] [CrossRef] [Google Scholar]
  16. Evans, I. N., Primini, F. A., Glotfelty, K. J., et al. 2010, ApJS, 189, 37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Guertler, J., Henning, T., Kruegel, E., & Chini, R. 1991, A&A, 252, 801 [NASA ADS] [Google Scholar]
  18. Hasegawa, T. I., & Mitchell, G. F. 1995, ApJ, 451, 225 [NASA ADS] [CrossRef] [Google Scholar]
  19. Hennebelle, P., Commerçon, B., Joos, M., et al. 2011, A&A, 528, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Henning, T., Pfau, W., & Altenhoff, W. J. 1990, A&A, 227, 542 [NASA ADS] [Google Scholar]
  21. Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478 [NASA ADS] [CrossRef] [Google Scholar]
  22. Hutawarakorn, B., & Cohen, R. J. 2005, MNRAS, 357, 338 [NASA ADS] [CrossRef] [Google Scholar]
  23. Jiménez-Serra, I., Zhang, Q., Viti, S., Martín-Pintado, J., & de Wit, W.-J. 2012, ApJ, 753, 34 [NASA ADS] [CrossRef] [Google Scholar]
  24. Johnston, K. G., Shepherd, D. S., Aguirre, J. E., et al. 2009, ApJ, 707, 283 [NASA ADS] [CrossRef] [Google Scholar]
  25. Johnston, K. G., Keto, E., Robitaille, T. P., & Wood, K. 2011, MNRAS, 836 [Google Scholar]
  26. Keto, E. 2002, ApJ, 580, 980 [NASA ADS] [CrossRef] [Google Scholar]
  27. Kim, S., Martin, P. G., & Hendry, P. D. 1994, ApJ, 422, 164 [NASA ADS] [CrossRef] [Google Scholar]
  28. Kleinmann, S. G., & Lebofsky, M. J. 1975, ApJ, 201, L91 [NASA ADS] [CrossRef] [Google Scholar]
  29. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  30. Krumholz, M. R., Klein, R. I., & McKee, C. F. 2012, ApJ, 754, 71 [NASA ADS] [CrossRef] [Google Scholar]
  31. Lada, C. J., & Fich, M. 1996, ApJ, 459, 638 [NASA ADS] [CrossRef] [Google Scholar]
  32. Lada, C. J., Thronson, Jr., H. A., Smith, H. A., Schwartz, P. R., & Glaccum, W. 1984, ApJ, 286, 302 [NASA ADS] [CrossRef] [Google Scholar]
  33. Lucy, L. B. 1999, A&A, 344, 282 [NASA ADS] [Google Scholar]
  34. Maeder, A., & Meynet, G. 1989, A&A, 210, 155 [NASA ADS] [Google Scholar]
  35. Marengo, M., Jayawardhana, R., Fazio, G. G., et al. 2000, ApJ, 541, L63 [NASA ADS] [CrossRef] [Google Scholar]
  36. McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  37. Minchin, N. R., Hough, J. H., McCall, A., et al. 1991, MNRAS, 251, 508 [NASA ADS] [Google Scholar]
  38. Minh, Y. C., & Yang, J. 2008, J. Korean Astron. Soc., 41, 139 [NASA ADS] [CrossRef] [Google Scholar]
  39. Miville-Deschênes, M.-A., & Lagache, G. 2005, ApJS, 157, 302 [NASA ADS] [CrossRef] [Google Scholar]
  40. Mueller, K. E., Shirley, Y. L., Evans, II, N. J., & Jacobson, H. R. 2002a, in Hot Star Workshop III: The Earliest Phases of Massive Star Birth, eds. P. Crowther, ASP Conf. Ser., 267, 395 [Google Scholar]
  41. Mueller, K. E., Shirley, Y. L., Evans, II, N. J., & Jacobson, H. R. 2002b, ApJS, 143, 469 [NASA ADS] [CrossRef] [Google Scholar]
  42. Myers, P. C. 2011, ApJ, 743, 98 [NASA ADS] [CrossRef] [Google Scholar]
  43. Panagia, N. 1973, AJ, 78, 929 [NASA ADS] [CrossRef] [Google Scholar]
  44. Podio, L., Bacciotti, F., Nisini, B., et al. 2006, A&A, 456, 189 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Poetzel, R., Mundt, R., & Ray, T. P. 1992, A&A, 262, 229 [NASA ADS] [Google Scholar]
  46. Preibisch, T., Balega, Y. Y., Schertl, D., & Weigelt, G. 2003, A&A, 412, 735 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Reynolds, S. P. 1986, ApJ, 304, 713 [NASA ADS] [CrossRef] [Google Scholar]
  48. Robitaille, T. P. 2011, A&A, 536, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Robitaille, T. P., Whitney, B. A., Indebetouw, R., & Wood, K. 2007, ApJS, 169, 328 [NASA ADS] [CrossRef] [Google Scholar]
  50. Rodriguez, L. F., & Canto, J. 1983, Rev. Mex. Astron. Astrofis., 8, 1631 [Google Scholar]
  51. Rodríguez, L. F., Gómez, Y., & Tafoya, D. 2007, ApJ, 663, 1083 [NASA ADS] [CrossRef] [Google Scholar]
  52. Rosenblum, B., Nethercot, A. H., & Townes, C. H. 1958, Phys. Rev., 109, 400 [NASA ADS] [CrossRef] [Google Scholar]
  53. Rygl, K. L. J., Brunthaler, A., Sanna, A., et al. 2012, A&A, 539, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Sanna, A., Reid, M. J., Carrasco-González, C., et al. 2012, ApJ, 745, 19 [NASA ADS] [CrossRef] [Google Scholar]
  55. Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269 [Google Scholar]
  56. Scoville, N. Z., Sargent, A. I., Sanders, D. B., et al. 1986, ApJ, 303, 416 [NASA ADS] [CrossRef] [Google Scholar]
  57. Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
  58. Sloan, G. C., Kraemer, K. E., Price, S. D., & Shipman, R. F. 2003, ApJS, 147, 379 [NASA ADS] [CrossRef] [Google Scholar]
  59. Tamura, M., & Yamashita, T. 1992, ApJ, 391, 710 [NASA ADS] [CrossRef] [Google Scholar]
  60. Terebey, S., Shu, F. H., & Cassen, P. 1984, ApJ, 286, 529 [NASA ADS] [CrossRef] [Google Scholar]
  61. Tofani, G., Felli, M., Taylor, G. B., & Hunter, T. R. 1995, A&AS, 112, 299 [NASA ADS] [Google Scholar]
  62. Trinidad, M. A., Curiel, S., Cantó, J., et al. 2003, ApJ, 589, 386 [NASA ADS] [CrossRef] [Google Scholar]
  63. Ulrich, R. K. 1976, ApJ, 210, 377 [NASA ADS] [CrossRef] [Google Scholar]
  64. van der Tak, F. F. S., & Menten, K. M. 2005, A&A, 437, 947 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. van der Tak, F. F. S., van Dishoeck, E. F., Evans,II, N. J., Bakker, E. J., & Blake, G. A. 1999, ApJ, 522, 991 [NASA ADS] [CrossRef] [Google Scholar]
  66. van der Tak, F. F. S., Walmsley, C. M., Herpin, F., & Ceccarelli, C. 2006, A&A, 447, 1011 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. van der Wiel, M. H. D., van der Tak, F. F. S., Spaans, M., et al. 2011, A&A, 532, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Wang, K.-S., van der Tak, F. F. S., & Hogerheijde, M. R. 2012, A&A, 543, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Whitney, B. A., Wood, K., Bjorkman, J. E., & Wolff, M. J. 2003, ApJ, 591, 1049 [Google Scholar]
  70. Whitney, B. A., Indebetouw, R., Bjorkman, J. E., & Wood, K. 2004, ApJ, 617, 1177 [NASA ADS] [CrossRef] [Google Scholar]
  71. Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191 [NASA ADS] [CrossRef] [Google Scholar]
  72. Wood, K., Wolff, M. J., Bjorkman, J. E., & Whitney, B. 2002, ApJ, 564, 887 [NASA ADS] [CrossRef] [Google Scholar]
  73. Yamashita, T., Sato, S., Tamura, M., et al. 1987, PASJ, 39, 809 [NASA ADS] [Google Scholar]
  74. Yorke, H. W. 2002, in Hot Star Workshop III, The Earliest Phases of Massive Star Birth, ed. P. Crowther, ASP Conf. Ser., 267, 165 [Google Scholar]

Appendix A: Mass determination from C18O emission

Starting from Eq. (A.1) from the appendix in Scoville et al. (1986):

N total = 3 k 8 π 3 B e μ 2 e h B e J l ( J l + 1 ) / kT ( J l + 1 ) T + h B e / 3 k ( 1 exp ( / k T ) )   0 τ υ   d υ \appendix \setcounter{section}{1} \begin{equation} N_{\rm total}= \frac{3k}{8 \pi^3 B_{\rm e} {\mu}^2} \frac{{\rm e}^{hB_{\rm e} J_{\rm l}(J_{\rm l}+1)/kT}}{(J_{\rm l}+1)} \frac{T+hB_{\rm e}/3k}{\left(1- \exp{\left(-h\nu/kT\right)} \right)} ~ \int_0^{\infty} \tau_{\upsilon} ~{\rm d}\upsilon \label{scovilleA1} \end{equation}(A.1)

where k is the Boltzmann constant, Be is the rotational constant for the molecule, μ is the dipole moment of the molecule in e.s.u. (1 Debye = 10-18 e.s.u.), h is the planck constant, Jl is the rotational quantum number for the lower energy level, T is the excitation temperature assumed to characterise the populations in all of the energy levels of the gas, ν is the frequency of the emission and τυ is the optical depth over the Doppler line profile as a function of velocity υ.

The solution to the equation of radiative transfer (assuming background terms are negligible) is

Tb=/ke/kT1(1eτυ).\appendix \setcounter{section}{1} \begin{equation} T_{\rm b} = \frac{h\nu/k}{{\rm e}^{h\nu/kT} -1} \left(1 - {\rm e}^{-\tau_{\upsilon}}\right). \label{solution} \end{equation}(A.2)Multiplying Eq. (A.1) by the ratio of the LHS to RHS of Eq. (A.2) and simplifying the terms containing e/kT gives Ntotal=Ntotal Tb(e/kT1)/k(1eτυ)=3k28π3Beμ2ehBeJl(Jl+1)/kT(Jl+1)T+hBe/3ke/kT0Tb τυ(1eτυ) dυ.\appendix \setcounter{section}{1} \begin{eqnarray} \label{nbar} {N}_{\rm total} = &&\, \frac{ {N}_{\rm total} ~ T_{\rm b} ({\rm e}^{h\nu/kT} -1)} {h\nu/k (1 - {\rm e}^{-\tau_{\upsilon}}) }\notag\\ = && \frac{3k^2}{8 \pi^3 B_{\rm e} {\mu}^2 h \nu} \frac{{\rm e}^{hB_{\rm e} J_{\rm l}(J_{\rm l}+1)/kT}}{(J_{\rm l}+1)} \frac{T+hB_{\rm e}/3k}{{\rm e}^{-h\nu/kT}} \int_0^{\infty} \frac{T_{\rm b} ~ \tau_{\upsilon}}{(1-{\rm e}^{-\tau_{\upsilon}})} ~{\rm d}\upsilon. \end{eqnarray}(A.3)

The total mass of gas in the source is given by:

M H 2 = N total [ H 2 12 CO ] [ 12 CO C 18 O ] μ G m H 2 π θ 2 4 d 2 \appendix \setcounter{section}{1} \begin{equation} M_{\rm H_2}= {N}_{\rm total} \left[ \frac{\rm H_2}{^{12}\rm{CO}}\right] \left[ \frac{^{12}\rm{CO}}{\rm{C}^{18}\rm{O}}\right] \mu_{\rm G} m_{\rm H_2} \frac{\pi \theta^2}{4} d^2 \label{massh2} \end{equation}(A.4)

where H212CO][\hbox{$\left[ \frac{\rm H_2}{^{12}\rm{CO}}\right]$} is the H2 to 12CO abundance ratio, 12COC18O][\hbox{$\left[ \frac{^{12}\rm{CO}}{\rm{C}^{18}\rm{O}}\right]$} is the 12CO to C18O abundance ratio, μG is the mean atomic weight of the gas, mH2 is the molecular mass of H2, θ is the angular diameter of a uniform disk source, and d is the distance to the source.

Substituting Ntotal (Eq. (A.3)) into Eq. (A.4) we find:

M H 2 = 3 k 2 8 π 3 B e μ 2 e h B e J l ( J l + 1 ) / kT ( J l + 1 ) T + h B e / 3 k e / kT 0 T b   τ υ ( 1 e τ υ )   d υ × [ H 2 12 CO ] [ 12 CO C 18 O ] μ G m H 2 π θ 2 4 d 2 . \appendix \setcounter{section}{1} \begin{equation} \begin{split} M_{\rm H_2}= \frac{3k^2}{8 \pi^3 B_{\rm e} {\mu}^2 h \nu} & \frac{{\rm e}^{hB_{\rm e} J_{\rm l}(J_{\rm l}+1)/kT}}{(J_{\rm l}+1)} \frac{T+hB_{\rm e}/3k}{{\rm e}^{-h\nu/kT}} \int_0^{\infty} \frac{T_{\rm b} ~ \tau_{\upsilon}}{(1-{\rm e}^{-\tau_{\upsilon}})} ~{\rm d}\upsilon \\ & \times \left[ \frac{\rm H_2}{^{12}\rm{CO}}\right] \left[ \frac{^{12}\rm{CO}}{\rm{C}^{18}\rm{O}}\right]\mu_{\rm G} m_{\rm H_2} \frac{\pi \theta^2}{4} d^2. \label{massh2_2} \end{split} \end{equation}(A.5)Constants used as well as several values specific to the (1−0) transition of C18O are listed below, expressed in c.g.s. units. Values for Be are taken from Rosenblum et al. (1958), and values for μ are taken from the JPL line catalog8. Jl=0ν=1.09782182×1011 HzBe=5.5135449×1010 cycles/sμ=0.11079 Debye=0.11079×10-18 e.s.u[H212CO]=104μg=1.36mH2=3.34524316×10-24g.\appendix \setcounter{section}{1} \begin{eqnarray*} & J_{\rm l} = 0 \\ & \nu = 1.09782182 \times 10^{11}~ \rm{Hz} \\ & B_{\rm e} = 5.5135449 \times 10^{10}~ \rm{cycles/s}\\ & \mu = 0.11079 ~\rm{Debye} = 0.11079 \times 10^{-18}~\rm{e.s.u} \\ & \left[ \frac{\rm H_2}{^{12}\rm{CO}}\right] = 10^4 \\ & \mu_{\rm g} = 1.36 \\ & m_{\rm H_2} = 3.34524316 \times 10^{-24}\,\rm{g}. \\ \end{eqnarray*}Entering these into Eq. (A.5) we obtain: MH2=1.88×10-7 [12COC18O](T+0.882)e5.27/T×θarcsec2dkpc2 0Tb τυ(1eτυ) dυ  M.\appendix \setcounter{section}{1} \begin{eqnarray} M_{\rm H_2} &&= \, 1.88 \times 10^{-7} ~ \left[ \frac{^{12}\rm{CO}}{\rm{C}^{18}\rm{O}}\right] \frac{(T+0.882)}{{\rm e}^{-5.27/T}} \notag\\ &&\times \, \theta_{\rm arcsec}^2 d_{\rm kpc}^2~ \int_0^{\infty} \frac{T_{\rm b} ~ \tau_{\upsilon}}{(1-{\rm e}^{-\tau_{\upsilon}})} ~{\rm d}\upsilon {~~M}_{\odot}. \end{eqnarray}(A.6)Similarly, the mass of gas may be obtained from the line flux density Sν, integrated over a map, Tb=c22k1ν2Sν,\appendix \setcounter{section}{1} \begin{eqnarray} T_{\rm b} &&= \frac{c^2}{2k}\frac{1}{\nu^2} S_{\nu}, \end{eqnarray}(A.7)

we find an alternative expression for the mass of the gas:

M H 2 = 2.75 × 10 -5   [ 12 CO C 18 O ] ( T + 0.882 ) e 5.27 / T × d kpc 2   0 S ν   τ υ ( 1 e τ υ )   d υ     M \appendix \setcounter{section}{1} \begin{eqnarray} \label{finalwithtau} M_{\rm H_2} &&= \,2.75 \times 10^{-5} ~ \left[ \frac{^{12}\rm{CO}}{\rm{C}^{18}\rm{O}}\right] \frac{(T+0.882)}{{\rm e}^{-5.27/T}} \notag\\ && \times \,d_{\rm kpc}^2~ \int_0^{\infty} \frac{S_{\nu} ~ \tau_{\upsilon}}{(1-{\rm e}^{-\tau_{\upsilon}})} ~{\rm d}\upsilon {~~M}_{\odot} \end{eqnarray}(A.8)

where Sν is measured in Jy.

If the assumption is made that the emitting gas is optically thin (τ ≪ 1) over the line, Eq. (A.8) simplifies to MH2=2.75×10-5 [12COC18O](T+0.882)e5.27/Tdkpc2 0Sν dυ  M.\appendix \setcounter{section}{1} \begin{eqnarray} M_{\rm H_2} = 2.75 \times 10^{-5} ~ \left[ \frac{^{12}\rm{CO}}{\rm{C}^{18}\rm{O}}\right] \frac{(T+0.882)}{{\rm e}^{-5.27/T}} d_{\rm kpc}^2~ \int_0^{\infty} S_{\nu} ~{\rm d}\upsilon {~~M}_{\odot}. \end{eqnarray}(A.9)The abundance ratio [12CO/C18O]\hbox{$\left[{^{12}\rm{CO}}/{\rm{C}^{18}\rm{O}}\right]$} can be calculated using the results given in Wilson & Rood (1994), where DGC\hbox{$D_{_{\rm GC}}$} is the galactocentric distance in kpc:

[ 12 CO C 18 O ] = ( 58.8 ± 11.8 ) D GC + ( 37.1 ± 82.6 ) . \appendix \setcounter{section}{1} \begin{equation} \left[ \frac{^{12}\rm{CO}}{\rm{C}^{18}\rm{O}}\right]= (58.8 \pm 11.8) D_{_{\rm GC}} +(37.1 \pm 82.6). \end{equation}(A.10)

All Tables

Table 1

Summary of VLA observations.

Table 2

Summary of radio interferometric observations.

Table 3

Observed near-IR to sub-millimetre fluxes for AFGL 2591, collated from the literature.

Table 4

Assumed ranges for model parameters as input to the genetic search algorithm.

Table 5

Parameters of the genetic algorithm best-fitting models.

Table 6

Measured properties of the observed 3.6 cm (8.4 GHz) and 7 mm (43 GHz) continuum sources.

Table 7

Overview of observed fluxes for VLA 1 to VLA 5.

Table 8

Derived properties of the inner part of the blue-shifted outflow traced by the C18O emission.

All Figures

thumbnail Fig. 1

SED of AFGL 2591, collated from the literature. The best-fitting models for an envelope with and without a disk (overplotted solid blue and dashed red lines respectively) are discussed in Sect. 5.1. The errors shown are those reset to 10% if the error in the measured flux was  <10%.

In the text
thumbnail Fig. 2

a) Model three-colour J, H and Ks band image for the envelope with disk model (RGB: K, H, and J bands respectively). b) Observed 2MASS three-colour image of AFGL 2591. c) Model J, H, and Ks band image for the envelope without disk model. The model images have been normalised to the total integrated fluxes given in Table 3, in order that the morphology of the emission can be easily compared. Stretch: red: Ks band, 130−300 MJy sr-1; green: H band, 95−150 MJy sr-1, blue: J band, 30−60 MJy sr-1.

In the text
thumbnail Fig. 3

Black error bars: normalised flux profiles for the three 2MASS bands, aligned with (top) and perpendicular to (bottom) the outflow axis. Blue and red lines: the profiles of the best-fitting models to both the observed SED and profiles, for envelope with disk and without disk models respectively. Vertical grey dashed lines mark the position of the central source.

In the text
thumbnail Fig. 4

Plots of the χ2 surface along the axis of each model parameter for each run, determined by varying each parameter while holding the other best-fitting parameters constant. Red dashed line: envelope without disk run, light blue line: envelope with disk run, and black line: control run.

In the text
thumbnail Fig. 5

a) Map of the 3.6 cm continuum emission surrounding AFGL 2591. Contours are -3, 3, 4, 5, 7, 10, 15, 20, 30, 40, 50...100  ×  rms noise = 30 μJy beam-1. Greyscale: –0.03 to 3.77 mJy beam-1 (1.2  ×  peak value). The synthesised beam is shown in the bottom left corner: 0.43′′    ×    0.40′′, PA = 43°. b) Close-up of the 3.6 cm continuum emission towards VLA 3. Contours, greyscale and beam as in panel a).

In the text
thumbnail Fig. 6

Map of the 7 mm continuum emission towards AFGL 2591. Contours are –4, 4, 7, 10, 15, 20, 25, 30  ×  rms noise = 56 μJy beam-1. Greyscale: –0.06 to 2.14 mJy beam-1 (1.2  ×  peak value). The inset panel shows a close-up of the 7 mm emission from VLA 3. The synthesised beam is shown in the bottom left corner of both images: 0.11′′    ×    0.11′′, PA = 43°.

In the text
thumbnail Fig. 7

Three-colour JHK′  Germini-North image of AFGL 2591 overlaid with the 3.6 cm contours from Fig. 5. Stretch of Gemini image: red, K′ band: 150−2500 MJy sr-1; green, H band: 150−500 MJy sr-1; blue, J band: 40−80 MJy sr-1.

In the text
thumbnail Fig. 8

Radio-SEDs of AFGL 2591-VLA 3 and VLA 1. The SED of AFGL 2591-VLA 3 is show in black error bars, and the resolved fluxes of AFGL 2591-VLA 3 from Marengo et al. (2000) are shown as red squares. Green squares show the fluxes of VLA 1 – errors are not shown as they are smaller than the markers. The black and green solid lines show greybody fits to the SEDs of AFGL 2591-VLA 3 between 100 and 850 μm, and VLA 1 at mid-IR fluxes, respectively (see Sects. 5.2.1 and 5.2.2). The green dashed line shows the fit to the long-wavelength fluxes for VLA 1.

In the text
thumbnail Fig. 9

Continuum emission towards AFGL 2591 at 106.7 and 110.5 GHz (~2.8 and 2.7 mm respectively), shown in both contours and greyscale. The rms noise in the images is σ = 2.2 and 2.1 mJy beam-1 respectively. The greyscales extend from −3    ×    σ to 65.5 and 54.0 mJy beam-1. Contours are at −3,3,4,5,7,10,15   and   20    ×    σ. The crosses show the positions of VLA 1 to 3. The synthesised beams are 4.9′′    ×    4.1′′ PA = −4.8° and 4.3′′    ×    3.5′′ PA = 95° respectively.

In the text
thumbnail Fig. 10

C18O (black line) and 13CO (grey line) spectral profiles measured at the position of AFGL 2591-VLA 3: 20h29m24.s\hbox{$\fs$}88 + 40°11′19.\hbox{$\farcs$}5 (J2000).

In the text
thumbnail Fig. 13

C18O channel map at 0.3 km s-1 resolution between –9.0 and –2.7 km s-1. The rms in the map is σ = 0.1 Jy beam-1. The map peak flux is 1.1 Jy beam-1. Contours are at –3, 3, 4, 5, 6, 7, 8, 9, 10, 11    ×    σ. The synthesised beam is shown in the bottom left-hand corner (4.5′′    ×     3.6′′, PA 93°). A scale size of 50 000 AU is represented by a bar in the bottom right-hand panel. The cross shows the position of the central source of AFGL 2591 from the Gemini near-IR J band image. The rest velocity of the cloud is –5.7 km s-1.

In the text
thumbnail Fig. 14

Shown in colour-scale: intensity-weighted first moment map of the C18O emission towards AFGL 2591. The corresponding velocities are shown on the colour-bar on the right. The black open circle represents the position and size of VLA 1 at 3.6 cm. The crosses indicate the position of VLA 2 and VLA 3. A scale size of 20 000 AU is represented by a bar in the upper left-hand corner. Shown in contours: the C18O integrated intensity map from –2.3 to –8.7 km s-1. The peak flux of the integrated intensity map is 2.5 Jy beam-1 km s-1. Contours are at –3, 3, 4, 5, 6, 7    ×    σ = 0.33 Jy beam-1 km s-1. The synthesised beam is shown in the lower right-hand corner (4.5′′ ×  3.6′′, PA 93°). The rest velocity of the cloud is –5.7 km s-1.

In the text
thumbnail Fig. 15

Gemini North three-colour JHK′ image of AFGL 2591 overlaid with contours of C18O emission integrated over –4.3 to −3.3 km s-1 (red) and –8.0 to –7.0 km s-1 (blue). The three arcs in the east of the image show the positions of bow shocks visible in the Gemini North image, which are too faint to be seen directly in the figure. The dashed line shows the direction of the ionized jet observed at 3.6 cm (PA 100°), and yellow diamonds show roughly the peak positions of the H2 knots observed in Tamura & Yamashita (1992). Red-shifted contours are –3, 3, 4, 5, 6  ×  0.1 Jy beam-1 km s-1 and blue-shifted contours are –3, 3, 4, 5, 6, 7  ×  0.1 Jy beam-1 km s-1. The synthesised beam is shown in the lower left-hand corner (4.5′′ ×  3.6′′, PA 93°). Near-IR image stretch as in Fig. 7.

In the text
thumbnail Fig. 11

13CO channel map at 0.3 km s-1 resolution between –12.0 and –5.7 km s-1. The rms in the map is σ = 0.1 Jy beam-1. The map peak flux is 1.3 Jy beam-1. Contours are at –3, 3, 5, 10, 15     ×    σ. The synthesised beam is shown in the bottom left-hand corner (4.4′′    ×     3.7′′, PA 96°). A scale size of 0.5 pc is represented by a bar in the bottom right-hand panel. The cross shows the position of the central source of AFGL 2591 from the Gemini near-IR J band image. The rest velocity of the cloud is –5.7 km s-1.

In the text
thumbnail Fig. 12

13CO channel map at 0.3 km s-1 resolution between –5.7 and 0.6 km s-1. The rms in the map is σ = 0.1 Jy beam-1. The map peak flux is 1.3 Jy beam-1. Contours are at –3, 3, 5, 10, 15  ×    σ. The synthesised beam is shown in the bottom left-hand corner (4.4′′    ×    3.7′′, PA 96°). A scale size of 0.5 pc is represented by a bar in the bottom right-hand panel. The cross shows the position of the central source of AFGL 2591 from the Gemini near-IR J band image. The rest velocity of the cloud is –5.7 km s-1.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.