A MASSIVE PROTOSTAR FORMING BY ORDERED COLLAPSE OF A DENSE, MASSIVE CORE

Massive stars impact many areas of astrophysics. In most galactic environments they dominate the radiative, mechanical, and chemical feedback on the interstellar medium, thus regulating the evolution of galaxies. Many low-mass stars form in clusters along with massive stars, and their protoplanetary disks can be affected by this feedback also. There is some evidence that our own solar system was influenced in this way (Adams 2010). Despite this importance, there is no consensus on the basic formation mechanism of massive stars. Theories range from Core Accretion, i.e., a scaled-up version of low-mass star formation (e.g., the Turbulent Core Model of McKee & Tan 2003, hereafter MT03), to Competitive Accretion at the centers of forming star clusters (Bonnell et al. 2001; Wang et al. 2010), to Stellar Collisions (Bonnell et al. 1998).

The Core Accretion theory predicts the existence of an envelope-fed, central accretion disk and relatively ordered and collimated bipolar outflows powered by accretion around a massive protostar. The other formation mechanisms predict a much less ordered geometry of gas and dust surrounding the protostar.

Collimated bipolar outflows have been observed from massive protostars (e.g., Beuther et al. 2002). Outflows may limit the formation efficiency (mass ratio of final star to initial gravitationally bounded protostellar core) to ∼0.5 (Matzner & McKee 1999), since they expel material along polar directions. The resulting low-density cavities have been proposed to significantly affect the appearance of the protostar in the mid-IR (MIR) (De Buizer 2006; De Buizer et al. 2012), and this is seen in the radiative transfer (RT) calculations of Zhang & Tan (2011, hereafter ZT11) and Zhang et al. (2013, hereafter ZTM13) using the RT code developed by Whitney et al. (2003).

The Turbulent Core Model also relates core properties to the surrounding self-gravitating star-cluster-forming clump. For a marginally unstable pre-stellar core of a given mass, its size, density, and subsequent accretion rate all depend on the pressure set by the mean mass surface density of the clump, Σ_cl. Observed values of Σ_cl are ∼0.1–1 g cm⁻² (MT03; Butler & Tan 2012), leading to high predicted accretion rates, ∼10⁻⁴–10⁻³ M_☉ yr⁻¹, potentially important for overcoming the high radiation pressure from the massive protostar.

G35.20−0.74 (hereafter G35.2) is a massive protostar associated with a well-defined outflow. At a distance of 2.2 kpc (Zhang et al. 2009), radio continuum emission shows that it contains an ultracompact H ii region at the center in the form of a collimated jet along the north–south (N–S) direction (Heaton & Little 1988; Gibb et al. 2003), which may be an example of an outflow-confined H ii region (Tan & McKee 2003). The outflow is also observed in the near-IR (NIR) (Fuller et al. 2001) with the northern cavity brighter than the southern one. Ground-based 10 and 18 μm observations only reveal the northern cavity (De Buizer 2006), suggesting it is inclined toward us. Outflow is also seen in CO (Gibb et al. 2003; Birks et al. 2006), including a wider angle flow in the NE–SW direction. It is unclear if this results from a separate driving source or is related to the N–S jet, but misaligned by interaction with the surrounding core/clump. In previous studies, an early B star was thought to be present at the center (Dent et al. 1985), surrounded by an envelope with mass ∼500 M_☉ (Little et al. 1998) to ∼3800 M_☉ (Paron & Weidmann 2010).

Here we present the first SOFIA-FORCAST images of G35.2. Combined with other multiwavelength data, we construct the spectral energy distribution (SED). We present a simple exploration of parameter space of RT models of massive protostars and use the observed SED and intensity profiles along the outflow axis to constrain the models.

The SOFIA-FORCAST (Herter et al. 2012) observations of G35.2 were performed on 2011 May 24 and 26 (UT) in the 11.3, 19.7, 31.5, and 37.1 μm filters at an altitude of 43,000 feet. The chopping secondary of SOFIA was driven at 4 Hz, with a matched chop and nod throw of 60'', and with a nod performed every 30 s. The final effective on-source exposure times were 909 s at 11.3 μm, 959 s at 19.7 μm, 4068 s at 31.5 μm, and 4801 s at 37.1 μm. The fluxes were calibrated by the SOFIA data reduction pipeline, which during the Basic Science period has an estimated 3σ ⩽ 20% calibration error in all filters (Herter et al. 2012).

The 11 and 31 μm images were observed simultaneously with the dichroic mode of FORCAST; and similarly for the 19 and 37 μm images. The relative offsets between the short- and long-wavelength arrays are known to be better than 1σ = 015. The best registration of the 19 μm image with respect to the 18.3 μm image from De Buizer (2006) was obtained using a χ² minimization algorithm. The latter was convolved with a Gaussian profile to reach the same resolution of the FORCAST image. Several features in the two images match up well after the registration increasing confidence in the technique. We adopt the pixel size (0768) as our conservative 1σ uncertainty in the relative offset between the 18.3 and 19.7 μm images. Considering the absolute astrometric uncertainties of ∼02 from the 18.3 μm image, we get an absolute astrometric uncertainty of 079 from the 19.7 μm image. The absolute astrometric uncertainty of the 37.1 μm image is then estimated to be better than 08. However, the 11.3 μm FORCAST image suffered from such low signal-to-noise ratio that it could not be used to perform similar registration for the 31.4 μm image. We had to use the χ² minimization algorithm to find the best registration of the 31.4 μm image and the astrometrically calibrated 37.1 μm image, noticing that they look remarkably similar morphologically. The absolute astrometric uncertainty in the 31.4 μm image is estimated to be 089.

Using standard stars observed throughout several Basic Science flights, an average FWHM for each wavelength was determined: 39 at 19.7 μm, 41 at 31.5 μm, and 45 at 37.1 μm. The images were then deconvolved using the maximum likelihood method (Richardson 1972; Lucy 1974), which yielded an approximately factor of two better resolution than the pre-deconvolved images. The deconvolved images compare favorably to simple unsharp masking of the original images, and in the case of the 19.7 μm image, it compares well to the Gemini 18.3 μm image, and thus the substructures revealed in the deconvolved images are believed to be real.

We also retrieved calibrated level 2 data (AOR 1342241157, HOBYS GT program) from Herschel-SPIRE and a single long integration SCUBA 850/450 μm jiggle-map from the James Clerk Maxwell Telescope (JCMT) archive at CADC.⁹ The SCUBA data were reduced and calibrated in a standard way using SURF (Sandell et al. 2001).

Figure 1 shows the FORCAST images at 31.5 and 37.1 μm along with the 11 and 18 μm images by Gemini-T-ReCS (De Buizer 2006) and 3.6 and 8 μm images by Spitzer-IRAC (GLIMPSE Survey; Churchwell et al. 2009). The airborne FORCAST images have a dynamic range of ∼4 orders of magnitude, similar to the space-based IRAC observations and much higher than the ground-based T-ReCS observations. The origin of the images is chosen to be the radio source G35.2N (R.A.= 18^h58^m13033, decl. = +1^d40^m3614). While the 11 and 18 μm images only show the northern side of the outflow, the FORCAST 31 and 37 μm images also reveal the fainter southern outflow cavity. This cavity is also seen in the IRAC images, again being fainter than the northern side.

Zoom In Zoom Out Reset image size

This trend agrees with the results of the RT simulations by ZT11 and ZTM13, which showed that at ∼40 μm, the near-facing outflow cavity is the dominant feature but, because of lower extinction, the structure along the outflow axis is more symmetric than at shorter wavelengths. In these models, the far-facing outflow appears in the NIR mainly due to the scattering by dust in and around the outflow cavity.

With their high dynamic range, the FORCAST images also show extended emission on larger scales, which is missing in the ground-based MIR images. The morphology of this extended emission has a resemblance to that seen with IRAC. This extended emission seems to be related to the cavity wall and/or a wider-angle part of the outflow. It is also possible that some, probably small, fraction of this emission is produced by other luminosity sources.

3.1. SED and Protostellar Model Global Parameters

We integrated over a 40'' × 50'' region to find total 31 and 37 μm fluxes of 527 and 922 Jy, respectively. No significant emission was seen outside of this area in either band. We constructed the SED using data from Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), GLIMPSE, Gemini-T-ReCS, IRAS-IGA (Cao et al. 1997), Herschel-SPIRE and JCMT-SCUBA (this paper), and IRAM 30 m Telescope (Mooney et al. 1995). Total 2MASS and IRAC fluxes were integrated from background-subtracted images in the same region as the FORCAST bands. The 10 and 18 μm fluxes (3.06 and 46.87 Jy; De Buizer 2006) were integrated with a 30'' aperture (again, there are no significant sources in the immediate vicinity). For IRAS bands, the beams are so large that the 40'' × 50'' region does not contain all the flux from this source, so we use larger apertures (65 × 3') that encircle all the emission from the source while avoiding fluxes from significant nearby sources. We obtained fluxes at 60 and 100 μm of 2.37 × 10³ and 4.35 × 10³ Jy. They may be considered as upper limits. For SPIRE and SCUBA images, we have carefully subtracted the backgrounds with two or more component Gaussian fitting. The fluxes are 546, 420, 123, and 16 Jy at 250, 350, 450, and 850 μm, respectively, with source sizes of 27'' × 18'', 25'' × 38 at 250 and 350 μm, and 22'' × 12'' in SCUBA bands. The SED is shown in Figure 2.

Zoom In Zoom Out Reset image size

We use the SED to constrain the properties of the protostar, its core, and the surrounding clump by fitting with the RT models of ZT11 and ZTM13. The models are developed self-consistently for a core embedded in a high-pressure clump forming a massive star via core accretion (MT03), including a treatment of rotating infall, an active accretion disk, and a disk wind filling the outflow cavity. We have further extended the ZTM13 models to also include emission from the surrounding clump (see Figure 3). Since we have not yet constructed a full model grid, we limit ourselves to achieving only a moderately good fit by exploring the parameters manually and with a quite coarse sampling of parameter space. We show four models that have relatively good matches to the SED in Figure 2 with their parameters listed in Table 1.

Zoom In Zoom Out Reset image size

Table 1. Parameters of the Four Fitted Models

Model	Core Mass	Mean Surface Density of	Protostellar	Outflow Opening	Inclination	Foreground	Luminosity
	Mc (M☉)	the Clump Σcl (g cm−2)	Mass, m* (M☉)	Angle θw, esc	θviewa	Extinction AV (mag)	Lbol (L☉)b
Model 1	240	1	34	51°	58°	8	2.2 × 105
Model 2	240	0.7	26	45°	51°	8	1.2 × 105
Model 3	240	0.7	22	40°	44°	15	9.0 × 104
Model 4	240	0.4	20	35°	43°	0	6.6 × 104

Notes. ^aInclination between the line of sight and the outflow axis. ^bTrue bolometric luminosity after allowing for foreground extinction and flashlight effect.

Download table as:  ASCIITypeset image

Model 1 (Figure 3) contains a core of initial mass M_c = 240 M_☉, and a protostar of mass m_* = 34 M_☉ forming at its center. The core is embedded in a clump with Σ_cl = 1 g cm⁻², thus setting its radius R_c = 0.11 pc (MT03). We distribute the clump material between 1 and 10 R_c with a power-law density profile (with index of −1.75; e.g., Mueller et al. 2002) that joins the core boundary smoothly and reaches Σ_cl ≃ 1 g cm⁻² by 10 R_c (counting two sides along a central sight line). Then the original clump mass (assuming spherical symmetry, i.e., before some material is swept up by the outflow) is 4.8 × 10³ M_☉. The size of the disk is set by assuming a ratio of rotational to gravitational energy of 2%, yielding an outer radius of 940 AU. Its mass is assumed to be one-third of the stellar mass. The accretion rate is 5.1 × 10⁻⁴ M_☉ yr⁻¹. The total luminosity is 2.2 × 10⁵ L_☉ including protostar, accretion boundary layer, and disk. An outflow cavity, filled with a disk wind, including dust along streamlines launched from the cooler disk outside of 7 AU, is included. Following ZT11 and ZTM13, the cavity opening angle θ_{w, esc} = 51°, to yield a final star formation efficiency from the core of 50%. The viewing angle between the sight line and the outflow axis is θ_view = 58°. At this viewing angle, if the received flux was assumed to be isotropic then a bolometric luminosity of 3.6 × 10⁴ L_☉ would be implied. After considering the foreground extinction as a free parameter to better fit the NIR data (A_V = 8 mag for Model 1) and applying an appropriate aperture size (∼50''), the inferred isotropic bolometric luminosity based on the received flux becomes 3.3 × 10⁴ L_☉. Three other models with smaller values of Σ_cl, m_*, and θ_view are also shown in Figure 2 (see also Table 1).

Overall, good agreement is achieved between these model SEDs and the observations. Small discrepancies in the NIR and IRAC bands may be due to patchy extinction and/or poorly modeled polycyclic aromatic hydrocarbon (PAH) features and/or thermal emission from transiently heated small dust grains. At 60 μm the observed IRAS flux is lower than the model (after convolving with the relatively broad filter response of IRAS 60 μm band) by ∼25%. These results show the SED fitting has a certain degeneracy in model parameter space, especially when the opening angle is not well constrained (ZTM13 found the outflow, especially the bright part, at 10 to 20 μm appears narrower than the actual cavity, so the MIR morphology may not be a reliable tracer of the total outflow opening angle). Models with narrower outflow cavities reprocess larger fractions of their luminosities and thus have smaller bolometric flashlight effect corrections for the considered viewing angles. Thus lower luminosity protostars are able to reproduce the observed FIR peak of the SED. Then, lower values of Σ_cl and θ_view are needed to match the MIR fluxes. Observations at 10 μm with higher spectral resolution will provide strong diagnostics for testing these models. However, even with a large uncertainty in the outflow opening angle, the protostellar mass is still estimated to be ∼20–34 M_☉ and Σ_cl ∼ 0.4–1 g cm⁻², i.e., G35.2 is a massive protostar forming from a massive core embedded in a high surface density clump.

3.2. Resolved Intensity Profiles along the Outflow Axis

The images of Model 1 are shown in Figure 3 for reference. The other three models selected from SED fitting show slightly narrower outflow cavities. We compare the intensity distribution along the outflow axis predicted by these models with observations in Figure 4. The axis direction of G35.2 is chosen via the radio continuum and MIR morphology (dashed line in Figure 1 with P.A. of 6°). The model profiles are all convolved with the corresponding instrument beams. At each offset from the center, we average over a perpendicular width of 2'' to estimate profiles for both model and data. All model and observed profiles are normalized to the average values over the near-facing regions that have significant emission. We also add constant background ambient intensities, which may be either due to instrumental noise (i.e., in the T-ReCS data) or from additional ambient interstellar material. Note that we have not attempted to use the lateral width of the cavities to constrain the models.

Zoom In Zoom Out Reset image size

We emphasize that a detailed model search to fit these profiles separately from the SED has not been performed, but still Model 1 agrees very well with the observations at 37, 31, 18, and 10 μm, producing the right peak positions and asymmetries of the two sides of the outflow. At 10 μm, Model 1 predicts higher extinction toward the center. A possible explanation could be the strength (depth and width) of the silicate absorption feature adopted in the dust model used in the RT calculations (ZTM13). In the IRAC bands, while a reasonable fit is achieved for the near-facing outflow cavity, the intensities of the far-facing side predicted by the model are too low compared to those observed. Possible reasons for this include (1) non-uniform, possibly patchy, extinction from the clump or the adopted uniform "foreground" extinction that may actually be relatively local to the source (e.g., a foreground extinction of A_V = 15 mag lowers 3.6 μm flux by a factor of two), (2) enhanced emission from PAHs and transiently heated small dust grains in and around the outflow cavity. Models 2, 3, and 4 also have good matches with the data in the near-facing and far-facing wings at 37, 31, 18, and 10 μm, but their profiles become too peaked toward the center due to the lower extinction of their lower Σ cores and clumps.

These results, especially the intensity profiles for Model 1 from ∼10–40 μm (for which the modeling has potentially fewer problems associated with small grain and PAH emission), again support the paradigm that G35.2 is a massive protostar forming from a high Σ core and clump. The match of the SED and intensity profiles with observations becomes less good when Σ_cl becomes lower. We estimate that Σ_cl should be no lower than 0.4 g cm⁻². The profiles also suggest that the accretion flow and outflows are relatively well ordered and symmetric.

SOFIA-FORCAST provides imaging with high dynamic range similar to space-based instruments at a unique wavelength region ∼30–40 μm, where lower extinction allows us to search for the predicted (ZT11; ZTM13) far-facing outflow cavity from a massive protostar forming from a high surface density core. At longer wavelengths the emission is predicted to become even more symmetric, being dominated by cooler dust in the core/clump. Our SOFIA-FORCAST observations at 31 and 37 μm did reveal emission from the far-facing outflow cavity of G35.2, which was too faint to detect by ground-based T-ReCs 11 and 18 μm observations.

We compiled the NIR to mm SED of G35.2. RT modeling of a massive protostar forming from a massive core bounded by a high Σ clump gave good agreement with this SED for four models. Depending on the outflow cavity opening angle (35°–50°), we found m_* ∼ 22–34 M_☉, L_bol ∼ (0.7–2.2) × 10⁵ L_☉, and Σ_cl ∼ 0.4–1 g cm⁻². Model 1 also produced intensity profiles along the outflow axis that fit the observations well at 10, 18, 31, and 37 μm, without need for extensive fine tuning. These results indicate G35.2 is a massive protostar, forming from high surface density core and clump, via relatively ordered, symmetric collapse and accretion. Powerful bipolar outflows are being launched and have cleared wide-angle cavities, which are also relatively symmetric.

A protostar with the luminosity estimated here (∼1 × 10⁵ L_☉) is expected to drive a CO outflow with momentum flux of ∼0.1 M_☉ yr⁻¹ km s⁻¹ (Richer et al. 2000), which is much larger than the observed value: even assuming that the larger NE–SW CO outflow is also driven by this source, the total momentum flux is still only ∼0.003 M_☉ yr⁻¹ km s⁻¹ (Gibb et al. 2003; Birks et al. 2006). This may indicate that the outflow from G35.2N, being partly ionized as indicated by the observed radio continuum jet, may be relatively deficient in CO emission. Note that the modest misalignment of the larger CO outflow with the radio and MIR jet may result from the interaction of the wider angle part of the outflow with a core/clump that does have some moderate lateral asymmetries.

We conclude that, at least in the case of G35.2, massive star formation appears to be well described by predictions of the Turbulent Core Model (MT03), even with a highly idealized realization of the model utilizing perfect axisymmetry and smooth density structures (see Cunningham et al. 2011 for a numerical simulation of this model including outflows). Such ordered accretion and outflow is not expected in models involving massive star formation via competitive accretion or stellar collisions, nor indeed in cores that are highly turbulent. This may indicate that magnetic fields are dynamically important in setting the structure of the core.

This work is based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NAS2-97001, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. Y.Z. acknowledges support from a Graduate School Fellowship from the University of Florida. J.C.T. acknowledges support from NSF CAREER grant AST-0645412 and NASA/USRA grant in support of SOFIA Basic Science observations. C.F.M. acknowledges support from NSF grants AST-0908553 and AST-1211729 and NASA grant NNX09AK31G.