A Mass determination from C{}^{18}O emission

The standard model of low-mass star formation applied to massive stars: a multi-wavelength picture of AFGL 2591

Key Words.:
Radiative transfer – Techniques: interferometric – (Stars:) circumstellar matter – Stars: formation – Stars: massive – ISM: jets and outflows


Context:While it is currently unclear from a theoretical standpoint which forces and processes dominate the formation of high-mass stars, and hence determine the mode in which they form, much of the recent observational evidence suggests that massive stars are born in a similar manner to their low-mass counterparts.

Aims:This paper aims to investigate the hypothesis that the embedded luminous star AFGL 2591-VLA 3 (2.3 L at 3.33 kpc) is forming according to a scaled-up version of a low-mass star formation scenario.

Methods:We present multi-configuration Very Large Array 3.6 cm and 7 mm, as well as Combined Array for Research in Millimeter Astronomy CO and 3 mm continuum observations to investigate the morphology and kinematics of the ionized gas, dust, and molecular gas around AFGL 2591. We also compare our results to ancillary Gemini North near-IR images, and model the near-IR to sub-mm Spectral Energy distribution (SED) and Two Micron All Sky Survey (2MASS) image profiles of AFGL 2591 using a Monte-Carlo dust continuum radiative transfer code.

Results:The observed 3.6 cm images uncover for the first time that the central powering source AFGL 2591-VLA 3 has a compact core plus collimated jet morphology, extending 4000 AU eastward from the central source with an opening angle of at this radius. However, at 7 mm VLA 3 does not show a jet morphology, but instead compact ( AU) emission, some of which (0.57 mJy of 2.9 mJy) is estimated to be from dust emission. The spectral index of AFGL 2591-VLA 3 between 3.6 cm and 7 mm was found to be between 0.4 and 0.5, similar to that of an ionized wind. If the 3.6 cm emission is modelled as an ionized jet, the jet has almost enough momentum to drive the larger-scale flow. However, assuming a shock efficiency of 10%, the momentum rate of the jet is not sufficient to ionize itself via only shocks, and thus a significant portion of the emission is instead likely created in a photoionized wind. The CO emission uncovers dense entrained material in the outflow(s) from these young stars. The main features of the SED and 2MASS images of AFGL 2591-VLA 3 are also reproduced by our model dust geometry of a rotationally flattened envelope with and without a disk.

Conclusions:The above results are consistent with a picture of massive star formation similar to that seen for low-mass protostars. However, within its envelope, AFGL 2591-VLA 3 contains at least four other young stars, constituting a small cluster. Therefore it appears that AFGL 2591-VLA 3 may be able to source its accreting material from a shared gas reservoir while still exhibiting the phenomena expected during the formation of low-mass stars.

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. hosokawa10), 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. yorke02) and ionization (e.g. keto02) may halt, decrease or alter accretion on to the star. However, in the earlier stages of protostellar evolution, hosokawa10 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 (mckee03), magnetic fields (e.g. hennebelle11; commercon11) or radiative heating and outflows (e.g. krumholz12), 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 (bonnell11; myers11), 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.5L at 3.33 kpc, lada84; henning90; rygl12). One of its most prominent features is a one sided conical reflection nebula observed in the near-IR (e.g minchin91; tamura92). Projected within the Cygnus-X star-forming complex, the distance to the source has recently been determined by trigonometric parallax measurements to be 3.330.11 kpc, more distant than previously assumed (between 1 and 2 kpc, e.g. poetzel92; hasegawa95; trinidad03; van-der-tak05). As will be described below, AFGL 2591 actually consists of several objects. However, as one source, AFGL 2591-VLA 3, dominates the 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 guertler91; van-der-tak99; mueller020 and de-wit09. Improving on these, preibisch03 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 trinidad03 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-tak99 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-wiel11 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 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, jimenez-serra12 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 campbell84; tofani95; trinidad03 and van-der-tak05 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 (trinidad03). A fourth radio continuum source has also been detected by campbell840 and tofani95 (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.

Wavelength Configuration Date of Program No. of Time Pointing centre Phase calibrator
observation antennas on-source (hr) R.A. Decl. flux density (Jy)
3.6 cm A 2007 Jul 26 AJ337 26 (22) 1.75 20 29 24.90 +40 11 21.00 (J2000) 1.49
3.6 cm B 2008 Jan 18 AJ337 26 (21) 1.96 20 29 24.90 +40 11 21.00 (J2000) 1.91
3.6 cm C 2008 Mar 9 AJ337 27 (23) 1.97 20 29 24.90 +40 11 21.00 (J2000) 2.13
3.6 cm D 2007 Apr 12 AJ332 26 (20) 0.37 20 29 24.90 +40 11 21.00 (J2000) 1.25
7 mm A 2002 Mar 25 AT273 27 (23) 1.97 20 27 35.95 +40 01 14.90 (B1950) 3.44
7 mm B 2008 Jan 18 AJ337 26 (24) 3.43 20 29 24.90 +40 11 21.00 (J2000) 3.15
7 mm C 2008 Mar 9 AJ337 27 (25) 4.18 20 29 24.90 +40 11 21.00 (J2000) 4.46
7 mm D 2007 Apr 28 AJ332 26 (17) 0.44 20 29 24.90 +40 11 21.00 (J2000) 1.64
Table 1: Summary of VLA observations

Knots of H and [SII] Herbig Haro objects have been detected toward AFGL 2591 (tamura92; poetzel92), suggesting the presence of shocked gas. These are coincident with an east-west bipolar outflow (e.g. lada84; hasegawa95), 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 minchin91 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 (yamashita87). At smaller scales, van-der-tak06 and wang12 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 jimenez-serra12, 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 (trinidad03; hutawarakorn05; sanna12). The Very Large Array (VLA)1 22 GHz water maser observations of trinidad03 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 sanna12 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), CO(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 CO, CO 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 Section 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-tak05. 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 for the 3.6 cm and 7 mm images respectively, and the synthesised beams were 0.43 by 0.40, P.A.=43 and 0.11 by 0.11, P.A.=43. These values, as well as the synthesised beam and map rms for each part of the CARMA observations described in Section 2.2, are summarised in Table 2.

2.2 Carma Co, CO and 3mm Continuum

CARMA observations at  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 CO(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.

The phase centre for AFGL 2591 was 202924.7 +401118.9 (J2000), and the total time on source was 2.5 hr. 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 (W. Kwon, private communication), 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 (D. Friedel, private communication), 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.

Array line/ Synthesised beam Map rms
cont. Size () P.A. () (mJy beam)
3.6 cm VLA-A to D cont. 0.430.40 43 0.030
7 mm VLA-A to D cont. 0.110.11 43 0.056
3 mm CARMA-D CO 4.43.7 96 100
3 mm CARMA-D CO 4.53.6 93 100
2.8 mm CARMA-D cont. 4.94.1 -4.8 2.2
2.7 mm CARMA-D cont. 4.33.5 95 2.1
Table 2: Summary of radio interferometric observations

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 (sloan03) 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 Section 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 Section 5.2.2, via a comparison of the SEDs of the resolved sources in the region.

Figure 1: The 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 Section 5.1. The errors shown are those reset to 10% if the error in the measured flux was 10%.
Wavelength Flux Density Aperture Description / Origin
(m) (Jy) radius
1.235 0.11 10% Irregular, 40 2MASS J band (1) 6
1.662 0.43 10% Irregular, 40 2MASS H band (1) 7
2.159 2.30 10% Irregular, 40 2MASS K band (1) 8
3.500 33.68 1.347 ISO apertures Sampled ISO-SWS spectrum (2)
3.822 68.98 2.759 ISO apertures Sampled ISO-SWS spectrum (2)
4.174 108.4 7.590 ISO apertures Sampled ISO-SWS spectrum (2)
4.558 127.2 8.906 ISO apertures Sampled ISO-SWS spectrum (2)
4.978 186.8 13.08 ISO apertures Sampled ISO-SWS spectrum (2)
5.436 233.8 16.37 ISO apertures Sampled ISO-SWS spectrum (2)
5.937 245.0 17.15 ISO apertures Sampled ISO-SWS spectrum (2)
6.483 273.7 19.16 ISO apertures Sampled ISO-SWS spectrum (2)
7.080 353.2 24.73 ISO apertures Sampled ISO-SWS spectrum (2)
7.732 406.2 28.43 ISO apertures Sampled ISO-SWS spectrum (2)
8.444 221.6 15.51 ISO apertures Sampled ISO-SWS spectrum (2)
9.221 78.69 5.508 ISO apertures Sampled ISO-SWS spectrum (2)
10.070 75.49 5.285 ISO apertures Sampled ISO-SWS spectrum (2)
10.997 154.3 10.80 ISO apertures Sampled ISO-SWS spectrum (2)
12.009 416.8 50.01 ISO apertures Sampled ISO-SWS spectrum (2)
13.115 583.4 70.01 ISO apertures Sampled ISO-SWS spectrum (2)
14.322 600.0 72.00 ISO apertures Sampled ISO-SWS spectrum (2)
15.641 556.8 66.81 ISO apertures Sampled ISO-SWS spectrum (2)
17.081 503.2 50.32 ISO apertures Sampled ISO-SWS spectrum (2)
18.653 542.3 54.23 ISO apertures Sampled ISO-SWS spectrum (2)
20.370 717.7 93.31 ISO apertures Sampled ISO-SWS spectrum (2)
22.245 886.1 115.2 ISO apertures Sampled ISO-SWS spectrum (2)
24.293 1118 145.3 ISO apertures Sampled ISO-SWS spectrum (2)
26.530 1391 180.9 ISO apertures Sampled ISO-SWS spectrum (2)
28.972 1728 293.8 ISO apertures Sampled ISO-SWS spectrum (2)
31.639 2118 465.9 ISO apertures Sampled ISO-SWS spectrum (2)
34.552 2413 530.9 ISO apertures Sampled ISO-SWS spectrum (2)
37.733 2810 618.3 ISO apertures Sampled ISO-SWS spectrum (2)
41.207 3384 744.5 ISO apertures Sampled ISO-SWS spectrum (2)
45.000 3879 853.5 ISO apertures Sampled ISO-SWS spectrum (2)
60 4830 290 6 IRAS 60m (3) 9
100 4390 290 6 IRAS 100m (3) 10
350 440 20% 60 CSO (4)
450 107 50% 60 JCMT (5) 11
850 15 20% 72 JCMT (5) 12

(given in parentheses in the column Description / Origin): (1) skrutskie06; (2) sloan03; (3) miville-deschenes05; (4) mueller02; (5) di-francesco08.

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 (skrutskie06) J, H, and K 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, H and J bands respectively. This Figure shows the conical near-IR reflection nebula (previously observed by e.g. kleinmann75; minchin91; tamura92), which is consistent with being a blue-shifted outflow lobe from the central source of AFGL 2591.

To obtain integrated J, H and K 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 202923.93 +401110.56 (J2000; easily visible as a point source in the Gemini North near-IR images presented in Section 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 bands respectively, with K 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 (P.A.=256, taken from sanna12, 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 Section 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.

Figure 2: a) Model three-colour J, H and K 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 K 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: K band, 130-300 MJy sr; green: H band, 95-150 MJy sr, blue: J band, 30-60 MJy sr.

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 release14. 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 (evans10), giving a positional accuracy of 0.6, and calibrated the flux level of the images in MJy sr 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 202924.86 +401119.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 Section 4. However, in Sections 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 robitaille11 and the genetic search algorithm used to model IRAS 20126+4104 in johnston11.

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 lucy99), 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. robitaille07; wood02a; whitney03). The axisymmetric three dimensional flared disk density is described between inner and outer radii and by


where =100 AU, is the cylindrical radius, is the height above the disk midplane and is set by the total disk mass , by integrating the disk density over , the cylindrical radius, and the azimuthal angle. The scaleheight increases with radius: where is the scaleheight at . We have assumed the parameter , derived for irradiated disks around low mass stars (dalessio98), and have taken , to provide a surface density of . This equation is essentially the same as equation (1) given in johnston11, however their second term is removed as this pertains to accreting -disks, whereas our disk is not accreting. Further, equation (1) is expressed in terms of instead of 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 (ulrich76; terebey84) with radius ,


where is the density scaling of the envelope, is the spherical radius, is the centrifugal radius, is the cosine of the polar angle (), and is the cosine of the polar angle of a streamline of infalling particles in the envelope as . The equation for the streamline is given by


which can be solved for . In our model, we assume the centrifugal radius is also the radius at which the disk forms .

We chose to describe the envelope density via a density scaling factor instead of the envelope accretion rate and stellar mass used in johnston11, 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 Section 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 , and a shape described by


where the shape of the cavity is determined by the parameter which we set to be , and is chosen so that the cavity half-opening angle at  AU is . The cavity was included by resetting the envelope density to inside the region defined by equation (4) only where the envelope density was initially larger.

The inner radius of the dust disk and envelope is expressed in terms of the dust destruction radius , which was found empirically by whitney04c to be


where the temperature at which dust sublimates is assumed to be K, and is the temperature of the star.

For the genetic algorithm, which is described in detail in Section 4.3 of johnston11, the size of the first generation was set to 1000, and the size of subsequent generations 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.

Figure 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.

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 vs. for the solar metallicity models of schaller92, 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 and a stellar surface gravity of were then used to select a model from the stellar atmosphere grid of castelli04. To sample the cavity half-opening angle , we first determined the relation between the half-opening angle and the inclination such that the equation of the cavity, equation (4), projected onto the plane of the sky produced the right opening angle observed in the Gemini near-IR images,


When sampling the half-opening angle, the inclination was first sampled and equation (6) was then used to determine an upper bound for , 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.

Parameter Description Value/Range Sampling
Stellar radius (R) 0.6 – 500 logarithmic
Stellar temperature (K) 5000 – 50,000 logarithmic
Envelope and disk inner radius (R) 1 – 100 logarithmic
Envelope outer radius (AU) logarithmic
Envelope density scaling factor (g cm) logarithmic
or Disk or centrifugal radius (AU) 10 – 10 logarithmic
Cavity half-opening angle at 10 AU (degrees) 0 – 90 linear
Inclination w. r. t. the line of sight (degrees) 0 – 90 linear
Cavity ambient density (g cm) logarithmic
Disk mass (M) 0.1 – 50 logarithmic
Disk scale-height at radius of 100 AU (AU) 0.2 – 20 logarithmic
Table 4: Assumed ranges for model parameters as input to the genetic search algorithm.

4.3 Model Fitting

As part of the genetic algorithm used here and described in johnston11, the models were fit to the data using the SED fitting tool of robitaille07, where , the external extinction, was a free parameter in the fit. The distance was set to be 3.33 kpc, and the visual extinction was allowed to vary between 0 and 40 magnitudes.

As described in robitaille07, 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, OSIA15) 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 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 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


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 Section 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, and , 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, will instead be referred to as the centrifugal radius of the envelope, , 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.

Parameter Description Value for envelope Value for envelope Value for envelope with
with disk model without disk model disk control model
Stellar radius (R) 90 200 52
Stellar temperature (K) 13,000 8900 18,000
Envelope and disk inner radius (R) 2.9 3.4 4.3
Envelope outer radius (AU) 180,000 420,000 790,000
Envelope density scaling factor (g cm) 2.4 1.1 5.6
or Disk or centrifugal radius (AU) 35,000 9500 2700
Cavity half-opening angle at 10 AU (degrees) 16 20 49
Inclination w. r. t. the line of sight (degrees) 59 46 54
Cavity ambient density (g cm) 2.1 2.0 1.0
Disk mass (M) 14 16
Disk scale-height at radius of 100 AU (AU) 11 11
Stellar luminosity (L) 2.3 2.3 2.4
Envelope density parameter (g cm AU) 1.6 1.0 8.0
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 value after convergence), are also given in Table 5.

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 (preibisch03).

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 and or are degenerate, as they both determine the envelope density, which can be seen from equation (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 tending to smaller values, tending to intermediate or large values, and the inclination preferring intermediate values. The surface for the cavity half-opening angle 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 Section 4.2. The shape of the surface for the cavity density is variable at higher values, but always has its minimum at low values, between and g cm. For the two runs which include a disk, the surfaces of the remaining parameters and are almost flat across the sampled ranges, although they slightly prefer values close to the upper limit. However, for the rises steeply toward the upper limit of the parameter range after its minimum at 11 AU.

Figure 4: Plots of the 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 for the envelope and disk inner radius . Using equation (5), we find = 34-39 AU for the three best-fitting models, and hence we expect values of  AU. This value is in agreement with the radius of the inner cavity adopted by jimenez-serra12, who in turn took this value from the disk of emission observed by preibisch03, likely to be tracing the inner rim of the envelope and disk. At 1 kpc, preibisch03 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-wiel11, 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 to be poorly constrained by our modelling, although we determined an upper limit for the cavity half-opening angle of in Section 4.2 (defined at 10 AU). wang12 adopted an opening angle of 60, and hence a half-opening angle of 30, in their modelling. The models used in jimenez-serra12 and van-der-wiel11 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 wang12, who in turn used the value determined by van-der-tak06. However, jimenez-serra12 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 Section 6.1, as a compromise between these two values.

We found low cavity densities were preferred for all runs, with best-fitting values between 10 and  g cm. In wang12 and van-der-tak99, 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 wang12 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 wang12 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 wang12 is more robust.

It was not possible to compare our best-fitting values for the disk scale-height to the model of wang12, as they instead used a 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 L for both the envelope with and without disk models respectively, which is in agreement with previous results (2.1-2.5L at 3.33 kpc, lada84; henning90; rygl12).

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 equation (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 minh08 is 210 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 minh08, 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 jimenez-serra12. Within this radius, which is half the determined size of their envelope 5400 AU, the total mass of the model of jimenez-serra12 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: 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 = 30Jy beam. Greyscale: -0.03 to 3.77 mJy beam (1.2 peak value). The synthesised beam is shown in the bottom left corner: 0.43 0.40, P.A. = 43. b) Close-up of the 3.6 cm continuum emission towards VLA 3. Contours, greyscale and beam as in panel a).
Figure 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. Greyscale: -0.06 to 2.14 mJy beam (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, P.A. = 43.
Source Peak position Peak flux Integrated flux Deconvolved P.A.
name (cm) (J2000) density (mJy beam) density (mJy) source size () (degrees)
VLA 1 3.6 20 29 24.59 +40 11 15.0 3.14 0.07 80 1.6 4.9 4.0 120
0.7 20 29 24.580 +40 11 14.50 0.46 0.06 64 3.2 1.8 1.6 90
VLA 2 3.6 20 29 24.52 +40 11 20.1 0.45 0.03 5.2 0.11 2.5 2.3 110
VLA 3 3.6 20 29 24.88 +40 11 19.5 0.72 0.03 1.52 0.03 1.9 0.84 90
0.7 20 29 24.882 +40 11 19.45 1.8 0.11 2.9 0.14 0.35 0.26 50
VLA 4 3.6 20 29 24.32 +40 11 28.0 0.39 0.03 0.99 0.02 1.1 0.78 90
VLA 5 3.6 20 29 24.00 +40 11 10.3 0.15 0.03 10.9 0.22 4.4 2.9 0
Table 6: Measured properties of the observed 3.6 cm (8.4 GHz) and 7 mm (43 GHz) continuum sources

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 campbell840 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 tofani95 and trinidad03. 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 202924.98 +401119.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. hasegawa95), and therefore these observations confirm the proposal of trinidad03 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  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-tak05, 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, sanna12 found that this emission most likely arises from the outflow cavity walls, which the masers they observed were also found to trace.

Table LABEL:vlafluxtable 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 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: 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; green, H band: 150 - 500 MJy sr; blue, J band: 40 - 80 MJy sr.

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 preibisch03. VLA 2 also appears to be coincident with the close binary discovered by preibisch03, and VLA 4 with a source in the Gemini North image at 202924.3 +401128 (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 202923.96 +401109.25 (J2000; which is the same source as 2MASS 20292393+4011105, previously mentioned in Section 3.2). The sources VLA 1 and VLA 3 are discussed in detail below.

The Central Source of AFGL 2591, VLA 3

Source Wavelength Integrated flux Ref.
name density (mJy)
VLA 1 6 cm 79 2.0 (1)
3.6 cm 80 1.6 (2)
7 mm 64 3.2 (2)
3.4 mm 87 1.4 (3)
2.83 mm 71 1.2 (3)
2.81 mm 76 15 (20%) (2)
2.7 mm 72 14 (20%) (2)
2.6 mm 93 2.5 (3)
18.0 m 7.96 0.8 10 (4)
12.5 m 1.59 0.2 10 (4)
11.7 m 1.84 0.2 10 (4)
VLA 2 6 cm 3.61 0.72 (1)
3.6 cm 5.2 0.11 (2)
VLA 3 6 cm 0.4 0.1 (1)
3.6 cm 1.52 0.03 (2)
1.3 cm 1.57 0.4 (5)
7 mm 2.9 0.14 (2)
3.4 mm 29.5 0.8 (3)
2.83 mm 38.7 0.7 (3)
2.81 mm 71 14 (20%) (2)
2.7 mm 79 16 (20%) (2)
2.6 mm 52.9 1.5 (3)
1.3 mm 151 4.5 (3)
18.0 m 7.53 0.75 10 (4)
12.5 m 7.45 0.75 10 (4)
11.7 m 4.40 0.44 10 (4)
VLA 4 6 cm 0.65 0.13 (1)
3.6 cm 0.99 0.02 (2)
VLA 5 3.6 cm 10.9 0.22 (2)
Table 7: Overview of observed fluxes for VLA 1 to VLA 5. References are: (1) campbell840, (2) this work (3) van-der-tak99, (4) marengo00, (5) trinidad03
Figure 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 marengo00 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 Sections 5.2.1 and 5.2.2). The green dashed line shows the fit to the long-wavelength fluxes for VLA 1.

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 LABEL:othersfluxtable (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.30.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-tak05 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 Section 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 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 sanna12 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 ) 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.430.03 mJy for 3.6 cm and 3.30.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.50.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.

Figure 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 respectively. The greyscales extend from to 65.5 and 54.0 mJy beam. Contours are at . The crosses show the positions of VLA 1 to 3. The synthesised beams are 4.94.1 P.A. and 4.3 3.5 P.A. respectively.

Vla 1

Figure 8 shows the SED of VLA 1 as green squares. This source was detected in the mid-IR by marengo00, as well as at several radio frequencies, for which the fluxes are listed in Table LABEL:othersfluxtable. The spectrum of VLA 1 appears to be flat at radio wavelengths, with a fitted spectral index of . 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, tofani95 and trinidad03 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 ().

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 was fit to the fluxes measured by marengo00. The best-fitting greybody has a temperature of 11114 K and peaks at 263.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 (%) at other wavelengths, and therefore should not significantly affect the SED and image modelling results given in Section 5.1. The mid-IR fluxes for only the resolved central source of AFGL 2591 from marengo00 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, CO and 3 mm Continuum

This section details the properties of the dense molecular gas surrounding AFGL 2591 inferred from the CO 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-tak99, van-der-tak06, and resolved further by wang12. 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 R.A. and 0.18 in declination to agree with the 7 mm position of MWC 349 A reported in rodriguez07. 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 LABEL:othersfluxtable.

Figure 10 presents CO 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 the Online Material as Figs. 11 and 12.

Figure 13 presents channel maps of the CO emission towards AFGL 2591 between -9.0 and -2.7 km s (where -5.7 km s is the rest velocity of the cloud), and Fig. 14 shows the CO 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. This can be seen in the CO 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 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.

In Fig. 15, integrated CO 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 preibisch03). 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 P.A. 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 hasegawa95 using single dish CO observations with 14.3 resolution, which extends in velocity from -45 to 35 km s. Therefore, the observations presented here trace the higher density but comparatively lower-velocity gas, within a few km s of the line centre. To minimise confusion due to overlap in velocity of various components, the contours in Fig. 15 only show integrated CO emission from -4.3 to -3.3 km s (red) and -8.0 to -7.0 km s (blue).

Figure 10: CO (black line) and CO (grey line) spectral profiles measured at the position of AFGL 2591-VLA 3: 202924.88 +401119.5 (J2000).


Figure 11: CO channel map at 0.3 km s resolution between -12.0 and -5.7 km s. The rms in the map is =0.1 Jy beam. The map peak flux is 1.3 Jy beam. Contours are at -3, 3, 5, 10, 15 . The synthesised beam is shown in the bottom left-hand corner (4.4 3.7, P.A. 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.


Figure 12: CO channel map at 0.3 km s resolution between -5.7 and 0.6 km s. The rms in the map is =0.1 Jy beam. The map peak flux is 1.3 Jy beam. Contours are at -3, 3, 5, 10, 15 . The synthesised beam is shown in the bottom left-hand corner (4.4 3.7, P.A. 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.

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 CO 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 hasegawa95. 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 202926.7 +401123 (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 jimenez-serra12 and the blue shifted outflow cone observed by sanna12, 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 (trinidad03), or emission from a spherical gravitationally confined HII region (van-der-tak05). 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 reynolds86, which describes the emission from a partially optically thick ionized jet (his equation 19):


where is the velocity of the jet, is the ionization fraction; is the mean particle mass per hydrogen atom of the ionized material, given by ; is the observed flux at the frequency ; is the spectral index; is the distance to the source; 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; is the temperature of the ionized gas; is the inclination measured from the line-of-sight and is a function of the spectral index and the dependance of optical depth on radius (see reynolds86, equation 17).

Assuming an isothermal, uniformly ionized jet with a density gradient of ;  km s for the velocity of the jet, found for the blue-shifted HH objects towards AFGL 2591 (poetzel92); (commonly found for low-mass sources, e.g. bacciotti99);  GHz; 0.21.2, the ratio of the maximum deconvolved width to the maximum radius;  K ; a spectral index between 0.4 and 0.5; and , taken from the results of the SED and near-IR image profile modelling presented in Section 5.1; 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.

Figure 13: CO channel map at 0.3 km s resolution between -9.0 and -2.7 km s. The rms in the map is =0.1 Jy beam. The map peak flux is 1.1 Jy beam. 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, P.A. 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.

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 (e.g. podio06). 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 hasegawa95 in CO(J=3-2), which corresponds well to the position and direction of the outflow lobe suggested by the ionized jet, red-shifted CO emission and bow shocks discussed in Section 5.3. For their optically thick case, hasegawa95 find a mass loss rate of  Myr at 1 kpc corresponding to  Myr 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 Myrkm s, which is very similar to the momentum rates determined by hasegawa95 for both the large and small scale red-shifted CO outflows (7.7 and 7.4 Myrkm s at 3.33 kpc with , scaled from 7.5 and 7.2 Myrkm s at 1 kpc with ). 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. curiel89; anglada96). The momentum rate of the jet can be expressed as:


where is the shock efficiency or ionization fraction, found to be 0.1 for low-mass sources (e.g. bacciotti99), is the initial velocity of the stellar wind or jet, taken to be 500 km s, and is the optical depth of the emitting gas. Here, the flux is measured at 3.6 cm. The optical depth can be determined using equation 6 of anglada98:


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 Myrkm s. Therefore the required momentum rate in the jet is larger than the value of 3.9 - 5.2 Myrkm s 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.

Figure 14: Shown in colour-scale: intensity-weighted first moment map of the CO 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 CO integrated intensity map from -2.3 to -8.7 km s. The peak flux of the integrated intensity map is 2.5 Jy beam km s. Contours are at -3, 3, 4, 5, 6, 7 =0.33 Jy beam km s. The synthesised beam is shown in the lower right-hand corner (4.5 3.6, P.A. 93). The rest velocity of the cloud is -5.7 km s.

Turning to the CO 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. lada96) 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 CO 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) instead of VLA 2, where it is approximately 6.3 km s. Thirdly, the blue-shifted emission seen in CO agrees very well with the high-velocity large-scale emission seen in Fig. 3 of hasegawa95. Therefore, if it were the case that the blue-shifted CO emission was created by VLA 2, this source would also be responsible for the large-scale blue-shifted outflow seen by hasegawa95. 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 CO 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),

Figure 15: Gemini North three-colour JHK image of AFGL 2591 overlaid with contours of CO emission integrated over -4.3 to -3.3 km s (red) and -8.0 to -7.0 km s (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 (P.A. 100), and yellow diamonds show roughly the peak positions of the H knots observed in tamura92. Red-shifted contours are -3,3,4,5,6 0.1 Jy beam km s and blue-shifted contours are -3,3,4,5,6,7 0.1 Jy beam km s. The synthesised beam is shown in the lower left-hand corner (4.5 3.6, P.A. 93). Near-IR image stretch as in Fig. 7.

where is the abundance ratio of to , is the temperature of the gas assuming local thermodynamic equilibrium, 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 Section 5.2.1, and the abundance ratio was calculated using the results given in wilson94, where is the galactocentric distance in kpc,


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

It should be noted that the shape of the 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 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 and kinetic energy in the flow, the final term of equation (11) was replaced by or respectively, where is the velocity of each channel, -5.7 km s is the rest velocity of the cloud, and is the inclination of the flow relative to the line of sight, taken to be from the results of the SED modelling presented in Section 5.1.

Parameter Value
Velocity range (km s) -8 to -5.3
Deprojected outflow length, (AU) 61,000
Outflow mass, (M) 133
Outflow momentum, (M km s) 164
Outflow energy, (M km s) 178
Outflow energy, (erg) 3.54 10
Average velocity relative to , (km s) -1.2
Dynamical timescale, (yr) 2.4 10
Table 8: Derived properties of the inner part of the blue-shifted outflow traced by the CO emission

The mass weighted average velocity of the flow was calculated as the momentum over the mass in the flow: , 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 47300 AU at 3.33 kpc. Therefore the actual length of the flow was calculated as . The dynamical timescale of the flow was determined using the equation , 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, beuther02), the timescale for the inner part of the flow seen in CO is at the upper end of these values. Scaling the blue-shifted dynamical timescale of  yr for the CO outflow observed by hasegawa95 to the new distance of 3.33 kpc, we obtain a dynamical timescale of  yr, similar to the value we find for the smaller-scale CO 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 Section 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-wiel11, as well as in the 80 or 270,000 AU diameter disk-like structure seen in the CS(1-0) maps of yamashita87. 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 johnston09, 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 Section 6.1: 0.77 - 1.0, the assumed wind velocity: 500 km s, and Table 4 of rodriguez83. For the above spectral type estimates, we note that the fluxes in Table LABEL:vlafluxtable 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 sanna12, 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 panagia73 and the theoretical HR diagram of maeder89. Thus if this cluster has a fully populated IMF, extrapolating a Krupa IMF (kroupa01) 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 panagia73 of , and a distance of 3.33 kpc, the spectral type of the calculated exciting star of VLA 1 corresponds to a 2MASS K band magnitude of 11.1. The detection limit of the 2MASS point source catalog in K band is 14.3 mag (skrutskie06), which indicates that there must be more than  mag of extinction at K band for this source to be obscured, corresponding to  mag of extinction in the V band (cardelli89). 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 202924.6 +401116.4 (J2000), suggested to be an edge-on disk by preibisch03, 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 relative to the ambient cloud. Further, it is interesting to note that the lack of CO 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 CO 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 ( 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 Section 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 Section 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 Section 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 at this radius. This jet is also roughly parallel to the outflow seen at larger scales by hasegawa95, 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 ( 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 preibisch03, 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 at wavelengths between 8 and 100m, 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 Section 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 Myr. 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 Section 5.3, the observed CO 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 hasegawa95. Thus although we cannot rule out that the that this emission is an outflow launched by VLA 2, there is substantial evidence that the CO 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 CO 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. arce07). 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.

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.

Appendix A Mass determination from CO emission

Starting from equation A1 from the appendix in scoville86:


where is the Boltzmann constant, is the rotational constant for the molecule, is the dipole moment of the molecule in e.s.u. (1 Debye= e.s.u.), is the planck constant, is the rotational quantum number for the lower energy level, 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


Multiplying equation 13 by the ratio of the LHS to RHS of equation (14) and simplifying the terms containing gives


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


where is the H to CO abundance ratio, is the CO to CO abundance ratio, is the mean atomic weight of the gas, is the molecular mass of , is the angular diameter of a uniform disk source, and is the distance to the source.