© 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    ×    10⁵   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  ≤ 10⁴AU. In addition, as well as observing a bi-conical outflow structure in ¹²CO (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.

Knots of H₂ 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)¹ 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 ¹³CO(1−0), C¹⁸O(1−0) and 3 mm continuum Combined Array for Research in Millimeter-wave Astronomy (CARMA²) 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 ¹³CO, C¹⁸O 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 Observatory³. 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)⁴ 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 ¹³CO, C¹⁸O and 3 mm continuum

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, ¹³CO(J = 1−0) and C¹⁸O(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 20^h29^m247 + 40°11′189 (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 V⁵. 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.

3.2. Near-IR 2MASS images

A three-colour 2MASS (Skrutskie et al. 2006) J, H, and K_s 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 K_s, 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 K_s 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 20^h29^m2393 + 40°11′1056 (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 K_s bands respectively, with K_s 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.

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 release⁶. 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 20^h29^m2486 + 40°11′195 (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 R^min and R_disk by (1)where R₀ = 100 AU, ϖ is the cylindrical radius, z is the height above the disk midplane and ρ₀ is set by the total disk mass M_disk, by integrating the disk density over z, ϖ the cylindrical radius, and φ the azimuthal angle. The scaleheight h(ϖ) increases with radius: where h₀ is the scaleheight at R₀. 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 R₀ 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 , (2)where ρ_env0 is the density scaling of the envelope, r is the spherical radius, R_c is the centrifugal radius, μ is the cosine of the polar angle (μ = cosθ), and μ₀ 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 (3)which can be solved for μ₀. In our model, we assume the centrifugal radius is also the radius at which the disk forms R_c    =    R_disk.

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 (4)where the shape of the cavity is determined by the parameter b_cav which we set to be b_cav = 1.5, and z₀ 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 R^min is expressed in terms of the dust destruction radius R_sub, which was found empirically by Whitney et al. (2004) to be (5)where the temperature at which dust sublimates is assumed to be T_sub = 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 R_V = 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 χ² value of the best fitting model decreased less than 5% in 20 generations.

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 T_ZAMS vs. log (0.9    ×    R_ZAMS) 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, (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 A_v, the external extinction, was a free parameter in the fit. The distance was set to be 3.33 kpc, and the visual extinction A_v 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, OSIA⁷) 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 K_s 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 K_s 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 χ² for each model given as input to the genetic algorithm was calculated as (7)where is the best-fitting reduced χ² value for the SED alone, and is the overall reduced χ² 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, M_disk and h₀, 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, R_disk will instead be referred to as the centrifugal radius of the envelope, R_c, 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.

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 χ² value after convergence), are also given in Table 5. The interstellar extinction A_V 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 χ² 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 χ² 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 χ²-values. From these we could therefore determine the χ² 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 χ² surfaces for stellar radius and temperature are different for each run, which is due to a degeneracy between these parameters. The χ² 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 χ² surface for the parameter orthogonal to the stellar luminosity, , 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 R_disk or R_c 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 and . Here is similar between the three runs, and has a sharp χ² surface, while the best-fitting 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, .

Looking at the remaining parameters, it can be seen that although the minima in their χ² surfaces are not sharp, the shape of their χ² surfaces is similar for the three runs, with R^min tending to smaller values, tending to intermediate or large values, and the inclination preferring intermediate values. The χ² 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 χ² 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. For the two runs which include a disk, the χ² surfaces of the remaining parameters M_disk and h₀ are almost flat across the sampled ranges, although they slightly prefer values close to the upper limit. However, for h₀ the χ² rises steeply toward the upper limit of the parameter range after its minimum at 11 AU.

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 R_sub for the envelope and disk inner radius R^min. Using Eq. (5), we find R_sub = 34−39 AU for the three best-fitting models, and hence we expect values of R^min ~ 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 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 10⁴ 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. 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 h₀ to the model of Wang et al. (2012), as they instead used a sin⁵θ 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    ×    10⁵   L_⊙ for both the envelope with and without disk models respectively, which is in agreement with previous results (2.1−2.5    ×    10⁵   L_⊙ at 3.33 kpc, Lada et al. 1984; Henning et al. 1990; Rygl et al. 2012).

The envelope density parameter 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    ×    10⁴   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

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 20^h29^m2498 + 40°11′193 (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°.

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 20^h29^m243 + 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 20^h29^m2396 + 40°11′0925 (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 arcsec² 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.

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.

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. ¹³CO, C¹⁸O and 3 mm continuum

This section details the properties of the dense molecular gas surrounding AFGL 2591 inferred from the C¹⁸O 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.

Figure 10 presents C¹⁸O and ¹³CO line profiles at the position of AFGL 2591-VLA 3. The ¹³CO 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 ¹³CO line profile seen in Fig. 10. As there was a significant amount of missing flux in the ¹³CO line, and the emission was distributed incoherently across the maps, it was not possible to interpret this emission. Therefore the ¹³CO data will not be discussed further, however channel maps of the ¹³CO emission are included in Figs. 11 and 12.

Figure 13 presents channel maps of the C¹⁸O 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 C¹⁸O 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 C¹⁸O 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 C¹⁸O 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 ¹²CO 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 C¹⁸O 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

In Fig. 15, the extrapolated direction of the ionized jet observed at 3.6 cm points directly through the red-shifted C¹⁸O 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 20^h29^m267 + 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)): (8)\arraycolsep1.75ptwhere υ is the velocity of the jet, x₀ is the ionization fraction; μ/m_p is the mean particle mass per hydrogen atom of the ionized material, given by 1/(1 + x₀); 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); x₀ = 0.1 (commonly found for low-mass sources, e.g. Bacciotti & Eislöffel 1999); v_m = 50 GHz; θ ~  0.2′′/1.2′′, the ratio of the maximum deconvolved width to the maximum radius; T = 10⁴ 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 ¹²CO(J = 3−2), which corresponds well to the position and direction of the outflow lobe suggested by the ionized jet, red-shifted C¹⁸O 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 ¹²CO 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: (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): (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.

Turning to the C¹⁸O 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 C¹⁸O 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 C¹⁸O 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 C¹⁸O 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 C¹⁸O 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),

(11)

where is the abundance ratio of ¹²CO to C¹⁸O, T is the temperature of the gas assuming local thermodynamic equilibrium, d_kpc is the distance to the source in kpc, and 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 was calculated using the results given in Wilson & Rood (1994), where is the galactocentric distance in kpc,

(12)

giving a value of for AFGL 2591 for a galactocentric distance of 8.5 kpc.

It should be noted that the shape of the C¹⁸O 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 C¹⁸O 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 M_outflowV and kinetic energy in the flow, the final term of Eq. (11) was replaced by or 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.

The mass weighted average velocity of the flow was calculated as the momentum over the mass in the flow:  ⟨ V ⟩  = M_outflowV/M_outflow, 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 = l_proj/sini. The dynamical timescale of the flow was determined using the equation t_dyn = 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 10⁴−10⁵ yr for L > 10⁴   L_⊙, Beuther et al. 2002), the timescale for the inner part of the flow seen in C¹⁸O is at the upper end of these values. Scaling the blue-shifted dynamical timescale of 6.3 × 10⁴ yr for the ¹²CO outflow observed by Hasegawa & Mitchell (1995) to the new distance of 3.33 kpc, we obtain a dynamical timescale of 2.1 × 10⁵ yr, similar to the value we find for the smaller-scale C¹⁸O 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 ¹²CO and ¹³CO 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 K_s band magnitude of 11.1. The detection limit of the 2MASS point source catalog in K_s band is 14.3 mag (Skrutskie et al. 2006), which indicates that there must be more than 3.2 mag of extinction at K_s 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 20^h29^m246 + 40°11′164 (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 C¹⁸O 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 C¹⁸O 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 (~10⁴ 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 C¹⁸O 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 C¹⁸O 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 C¹⁸O 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.

---

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

Appendix A: Mass determination from C¹⁸O emission

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

(A.1)

where k is the Boltzmann constant, B_e 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, J_l 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

(A.2)Multiplying Eq. (A.1) by the ratio of the LHS to RHS of Eq. (A.2) and simplifying the terms containing e^−hν/kT gives (A.3)

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

(A.4)

where is the H₂ to ¹²CO abundance ratio, is the ¹²CO to C¹⁸O abundance ratio, μ_G is the mean atomic weight of the gas, m_H2 is the molecular mass of H₂, θ is the angular diameter of a uniform disk source, and d is the distance to the source.

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

(A.5)Constants used as well as several values specific to the (1−0) transition of C¹⁸O are listed below, expressed in c.g.s. units. Values for B_e are taken from Rosenblum et al. (1958), and values for μ are taken from the JPL line catalog⁸. Entering these into Eq. (A.5) we obtain: (A.6)Similarly, the mass of gas may be obtained from the line flux density S_ν, integrated over a map, (A.7)

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

(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 (A.9)The abundance ratio can be calculated using the results given in Wilson & Rood (1994), where is the galactocentric distance in kpc:

(A.10)

Online material

All Tables

All Figures