HIERARCHICAL FRAGMENTATION AND JET-LIKE OUTFLOWS IN IRDC G28.34+0.06: A GROWING MASSIVE PROTOSTAR CLUSTER

There are different definitions of clump, core, and condensation in the literature when describing the spatial structure of dust continuum or molecular line emission. In the following discussion, we refer a clump as a structure with a size of ∼1 pc, a core as a structure with a size of ∼0.1 pc, and a condensation as a substructure of ∼0.01 pc within a core. A clump is capable of forming a cluster of stars, a core may form one or a small group of stars, and a condensation can typically form a single star or a multiple-star system. This nomenclature is consistent with the one adopted by Zhang et al. (2009).

Figure 1 presents the 0.88 mm continuum emission. The 069 × 064 resolution image resolves the clump into five groups of compact condensations, corresponding to the five dust cores discovered at 1.3 mm, namely SMA1, SMA2, SMA3, SMA4, and SMA5 (Zhang et al. 2009). Cores SMA1, SMA3, and SMA5 consist of one condensation. Cores SMA2 and SMA4 consist of three condensations, with one dominating over the other two in peak flux. In addition, a faint condensation is also revealed northeast of core SMA2. We tentatively assign this condensation to the SMA2 group. We name all the 10 condensations accordingly and label them in Figure 1. Refining the core positions using the positions of the dominant condensations, the five cores are spaced by a projected distance of 0.16 ± 0.02 pc (68 ± 10) and are well aligned at a position angle (P.A., east of north) of 48°. In cores SMA2 and SMA4, the three condensations are spaced by 0.03 ± 0.007 pc (13 ± 03).

Zoom In Zoom Out Reset image size

We estimate the mass of each condensation. Assuming optically thin dust emission, the dust mass can be estimated following

$$$ M_{\rm dust}=\frac{F_{\nu }d^2}{B_{\nu }(T_{\rm dust}){\kappa }_{\nu }}\,, $$$

where M_dust is the dust mass, F_ν is the continuum flux at frequency ν, d is the source distance, B_ν(T_dust) is the Planck function at dust temperature T_dust, and κ_ν = 10(ν/1.2 THz)^β cm² g⁻¹ is the dust opacity (Hildebrand 1983). We adopt the dust temperatures from Zhang et al. (2009), i.e., 16 K for SMA2 and 13 K for other cores. The dust opacity at millimeter/submillimeter wavelengths is uncertain. For dust in the interstellar medium, β = 2 (Draine & Lee 1984). Rathborne et al. (2010) derived β = 1.5 ± 0.3 for the entire P1 clump based on a global spectral energy distribution. We generated a dust opacity map by comparing the continuum image at 345 GHz and 230 GHz, and found β to vary from 1.3 to 2.5 across the map. Given the uncertainties, we adopt β = 1.5 for all the cores similar to Zhang et al. (2009) and discuss the effect of different β on the mass estimates.

Table 2 lists coordinates, size, peak flux, integrated flux, and estimated mass of all the 10 condensations. These parameters are determined by fitting a two-dimensional Gaussian function to the observed flux distribution. Fluxes are then corrected for primary beam attenuation. The masses of the condensations range from 1.4 to 10.6 M_☉ with an average of 4.8 M_☉. The mean size is 09 × 07, or 06 × 02 after deconvolution with the synthesized beam. The deconvolved size corresponds to about 3000 × 1000 AU at the source distance. The uncertainty in the mass estimation arises from several factors. First, the flux calibration contributes 15% in uncertainty. Second, the uncertainty of β contributes in a form of M_dust∝3.5^β. Third, the spatial filtering effect of radio interferometry preferentially picks up the compact structures and filters out the extended emission. Comparing images made from different array configurations, the EXT image recovered 20%–30% fluxes of that covered by COM and SUB images. This filtering effect explains why the masses reported in Table 2 are consistently lower than those in Zhang et al. (2009). The EXT image at 345 GHz represents the high-contrast structure in these cores—the small, compact condensations tracing the immediate surroundings of the protostars.

Our observations reveal fragmentation at different spatial scales in G28.34-P1. First, the 1 pc clump fragments into five cores (Zhang et al. 2009); second, two of the 0.1 pc cores fragment into even small condensations of 0.01 pc. The fragmentation at clump and core scales is consistent with a hierarchical fragmentation picture. The clump fragmentation is likely the result of initial physical conditions (density, temperature, turbulence, and magnetic fields), while the core fragmentation results from an increased density after the initial fragmentation. Based on the Very Large Array (VLA) NH₃ (1,1) and (2,2) spectra (Wang et al. 2008) and the IRAM 30m 1.3 mm continuum image (Rathborne et al. 2006), we measure an initial gas temperature of 13 K and density of 3 × 10⁵ cm⁻³, over a scale of 1 pc toward the P1 clump. These values yield a thermal Jeans mass of 0.8 M_☉. At the core scale, density increases to several 10 times the initial clump density. For instance, in SMA2 the core density measured from the 230 GHz image is 9.6 × 10⁶ cm⁻³, 32 times the initial clump density. This density, together with an increased gas temperature (16 K), leads to a thermal Jeans mass of 0.2 M_☉. These parameters are also listed in Table 3.

It is worth stressing that, in both clump and core fragmentation, the observed fragment (core/condensation) masses are significantly (∼10 times or more) larger than the expected thermal Jeans masses (Table 3), indicating that thermal pressure is not dominant in the fragmentation processes. Instead, turbulence may play an important role in supporting these large core/condensation masses, as discussed in Zhang et al. (2009) (see further discussion in Section 3.1.2). The 1σ rms sensitivity in Figure 1 is 0.8 mJy, corresponding to 0.2 M_☉, well below the thermal Jeans mass at the clump scale (0.8 M_☉). However, we did not detect any cores of 0.8 M_☉ in the filament. The non-detection of Jeans mass cores in the clump may suggest that low-mass cores are not significantly centrally peaked, thus they cannot be substantiated from the smooth emission in the clump.

The most intriguing feature in G28.34-P1 is the configuration of well aligned, regularly spaced cores. The alignment and regularity strongly suggest a mechanism that is responsible for shaping the molecular filament into what it is now. Similar features (but at different spatial scales) have been reported in a large number of nearby dark clouds (Schneider & Elmegreen 1979), in a few filamentary IRDCs (e.g., Jackson et al. 2010; Miettinen & Harju 2010), as well as in numerical simulations (Martel et al. 2006). It has been suggested that these fragments are most likely the results of gravitational collapse of a cylinder, as first proposed by Chandrasekhar & Fermi (1953) and followed up by Nagasawa (1987) (also see discussion in Jackson et al. 2010). The theory predicts that, under self-gravity, the gas in the cylinder will ultimately break up into pieces with a certain interval (λ_max) at which the instability grows the fastest and thus dominates the fragmentation process. Such "sausage" instability produces a chain of equally spaced fragments along the filament, with the spacing roughly the interval λ_max. The mass per unit length along the cylinder (or linear mass density) has a maximum value, $(M/l)_{\rm max} = 2v^2/G = 465\,(\frac{v}{\rm km\,s^{-1}})^2$ M_☉ pc⁻¹, where G is the gravitational constant. If the cylinder is supported by thermal pressure, v is the sound speed c_s; if, on the other hand, it is mainly supported by turbulent pressure, then v is replaced by the velocity dispersion σ.

We now compare the observational results with the theoretical predictions. In G28.34-P1, the FWHM line width is 1.7 km s⁻¹, measured from the NH₃ (1,1) spectrum (Wang et al. 2008; Zhang et al. 2009). The line width is equivalent to a velocity dispersion of  σ = 0.72 km s⁻¹ if the line profile is Gaussian. This velocity dispersion leads to a maximum mass density of 240 M_☉ pc⁻¹. The total core mass in P1 is 183 M_☉ over 0.8 pc along the filament (Zhang et al. 2009), implying a linear mass density of 230 M_☉ pc⁻¹. Taking into account the missing flux in the interferometer maps, the observations match the theoretical predictions well if G28.34-P1 is mainly supported by turbulence.

The separation between two adjacent fragments, λ_max, depends on the nature of the cylinder. For an incompressible fluid cylinder in the absence of magnetic field, λ_max = 11R, where R is the unperturbed cylindrical radius. And the instability manifests at about two free-fall timescales (see Section 3.2.2). In G28.34-P1, the average separation between cores is λ_max = 0.16 pc, implying an initial radius of R = 0.014 pc, several times smaller than the core diameter of 0.1 pc (Zhang et al. 2009). If magnetic field is involved, the strength is of the order of $B = 4\pi \rho R\sqrt{G}$ (see Equation (93) in Chandrasekhar & Fermi 1953). With n = 3 × 10⁵ cm⁻³ and R = 0.014 pc, we obtain a magnetic field strength of 0.16 mG. This is relatively small compared to magnetic field strengths observed in massive star formation regions (e.g., Girart et al. 2009).

For an isothermal gas cylinder, λ_max = 22H, where H = v(4πGρ_c)^−1/2 is the scale height, whereas ρ_c is the gas density at the center of the cylinder. The equation can be simplified as

$$\begin{equation*} \lambda _{\rm max} = \left\lbrace \begin{array}{@{}l} \dsty 0.15\, {\rm pc}\, \left(\frac{c_s}{\rm 0.21\, km\,s^{-1}}\right) \left(\frac{n}{3\times 10^5 \,{\rm cm^{-3}}}\right) \\ \qquad\dsty \mbox{\rm for thermal support,} \\\ms \dsty 0.13\, {\rm pc}\, \left(\frac{\sigma }{\rm 0.72\, km\,s^{-1}}\right) \left(\frac{n}{5\times 10^6 \,{\rm cm^{-3}}}\right) \\ \qquad \dsty \mbox{ for turbulent support.} \end{array} \right. \end{equation*}$$

In the above equations, the sound speed c_s of 0.21 km s⁻¹ is calculated using a gas temperature of 13 K; the volume density 3 × 10⁵ cm⁻³ is the averaged clump density, while 5 × 10⁶ cm⁻³ is the average of clump density (3 × 10⁵ cm⁻³) and core density (9.6 × 10⁶ cm⁻³) given in Table 3. We see that both scenarios are in agreement with the observed spacing of λ_max = 0.16 ± 0.02 pc if the assumptions are reasonable. However, since the density should be the central density of the cylinder, we suspect that the turbulent support is more likely dominant in G28.34-P1.

Comparison between observational results and theoretical predictions suggests that the fragmentation in the G28.34-P1 filament is well represented by a cylindrical collapse. The filament is likely supported mainly by turbulence rather than thermal pressure.

We measure the core density profiles in unresolved cores SMA1, SMA3, and SMA5. The measurements are performed in the visibility domain to avoid defects in image deconvolution. We combine the SUB, COM, and EXT data for analysis for this purpose. The combined visibility data sample the G28.34-P1 region over baseline lengths ranging from about 5 to 200 kλ, offering a fairly well sampled (u, v) measurement.

For a spherical core, the observed flux distribution over radial distance F(r) represents the radial density profile modified with temperature gradient. The dust temperature scales as T_dust∝r^−a (Scoville & Kwan 1976) with a = 0.33, if the core is internally heated. As long as Rayleigh–Jeans approximation holds, this temperature gradient contributes to the total flux in the form of ∝r^−a. Assuming a power-law density profile, ρ∝r^−b, and optically thin dust emission, the observed dust continuum flux integrated along the line of sight, F∝∫ρT_dustds, or

$$$ F \propto r^{-(a + b - 1)}\,, $$$

when (a + b) > 1. In the visibility domain, this is Fourier transformed into a form of

$$$ A_{uv} \propto S_{uv}^{(a + b - 3)}\,, $$$

where A_uv is the visibility amplitude and $S_{uv} = \sqrt{u^2 + v^2}$ is the (u, v) distance (Looney et al. 2000; Zhang et al. 2009).

Figure 2 plots the amplitude versus (u, v) distance for SMA1, SMA3, and SMA5. The amplitude is a vector average of the visibility data in a concentric annulus defined by a bin of 8 kλ. Least-squares fitting yields b = 2.09 ± 0.09 for SMA1, 1.97 ± 0.08 for SMA3, and 1.84 ± 0.10 for SMA5. The results are consistent with Zhang et al. (2009) measured from the 230 GHz continuum data, where they derived b = 2.1 ± 0.2 for SMA4. The density profiles are similar to that of an isothermal sphere in hydrostatic equilibrium, where the radial density scales as ∝r⁻². We note that the large error bars and discrepancy at long (u, v) distances reflect sensitivity limits and may be also partly due to unresolved weak sources in the vicinity.

Zoom In Zoom Out Reset image size

Among the entire 8 GHz SMA band, only CO (3–2) is detected above 3σ (1 σ ≈ 0.1 K at 1 km s⁻¹ resolution). Nevertheless, CO (3–2) reveals high-velocity outflowing gas in P1. Figure 3 presents the CO (3–2) channel maps, and Figure 4 plots the integrated blueshifted emission in blue contours and redshifted emission in red contours, superposed on the continuum emission. One can find in Figure 3 that the extended CO emission close to the systemic velocity (78.4 km s⁻¹) is filtered out. This effect actually helped us identify outflows at high velocities. It is intriguing that all the five dust cores are associated with CO outflows, the majority of which are bipolar, jet-like outflows.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The outflows are centered on the dust cores and are generally oriented across the major axis of the filament. The most prominent outflow is centered on core SMA2, with a blue lobe emanating toward southeast (SE; see channels of 51–71 km s⁻¹ in Figure 3) and a red lobe toward northwest (NW; 86–116 km s⁻¹). In addition, a minor blue lobe is also seen at NW (66–76 km s⁻¹). The outflow centered on SMA3 is in an east–west orientation, with its blue lobe emanating westbound (51–71 km s⁻¹) and its red lobe eastbound (86–96 km s⁻¹). SMA4 lies in the center of a NW–SE outflow, with a SE blue lobe (66–76 km s⁻¹) and a NW red lobe (86–91 km s⁻¹). SMA1 is associated with a quadrupolar outflow. Its two blue lobes are seen at 66–71 km s⁻¹, emanating toward SE and north, respectively. Its two red lobes are seen from 86 through 91 km s⁻¹ toward NW and south, respectively. The outflow centered on SMA5 has an orientation almost parallel to the SMA4 outflow. This outflow is weak compared to others, and its SE blue lobe and NW red lobe can be barely seen at 66–71 km s⁻¹ and 86–91 km s⁻¹, respectively. Orientations of all the outflows are sketched as arrows in Figure 4.

The driving sources of these outflows can be identified geometrically. Cores SMA1 to SMA5 are roughly located at the geometric centers of the outflows, respectively. While SMA1, SMA3, and SMA5 are unresolved, SMA2 and SMA4 are resolved into three condensations. In these two cases, we assign the strongest and/or central condensation as the driving source of the relevant outflow. Hence, the outflows are likely driven by condensations SMA1, SMA2a, SMA3, SMA4a, and SMA5.

The CO outflows shown in Figures 3 and 4 consist of pairs of knots which are geometrically symmetric with respect to the central protostar. As illustrated in Figure 4, outflows SMA2a and SMA3 have four pairs of knots, and outflow SMA4a has two pairs. The knots are equally spaced by about 0.16 pc. These knots in outflows may represent outbursts due to disk variability arising from episodic, unsteady accretion (e.g., HH 211, Lee et al. 2010; HH 80-81, Qiu & Zhang 2009). With the spacing and maximum outflow velocity v_max, we estimate the period between two outbursts to be about 4000 yr. This value lies in the range of roughly 10³–10⁴ yr used in episodic accretion models for low-mass stars (Baraffe & Chabrier 2010; Zhu et al. 2010).

The outflows show slight bending (Figure 4). In the outflow SMA3, the red lobe bends 5° toward the south, and the blue lobe bends 6° toward the south, making the outflow C-shaped. In outflows SMA2a and SMA4a, and probably SMA1, the knots seem to trace the propagation of a periodically wiggled, S-shaped molecular jet, most evident at the closest knotty pairs near the protostar. Part of this effect may be contaminated by multiple outflows (e.g., in outflow SMA2a). Nevertheless, the mirror symmetric (C-shaped) wiggles and point symmetric (S-shaped) wiggles may indicate different bending mechanisms like jet precession, orbital motion of the jet source, and Lorentz forces (Fendt & Zinnecker 1998; Masciadri & Raga 2002; Raga et al. 2009). Similar bending effects have been observed in a number of Herbig–Haro jets (e.g., HH 211, Lee et al. 2010) and some of the most collimated high-mass molecular jets (e.g., Su et al. 2007). Because the wiggles are not prominent in CO (3–2), we do not discuss this effect further in this paper.

We calculate  the physical parameters of each outflow, including mass, momentum, energy, as well as dynamical age and outflow rate. Assuming local thermodynamic equilibrium and optically thin CO emission in the line wings, we first derive the CO column density following Garden et al. (1991):

$$$ N_{\rm CO} ({\rm cm^{-2}}) = 4.81 \times 10^{12}\, (T_{\rm ex}+0.92)\, {\mathrm{exp}}\left(\frac{33.12}{T_{\rm ex}}\right)\int T_{\rm B}\,dv\,\,, $$$

where dv is the velocity interval in km s⁻¹, and T_ex and T_B are excitation temperature and brightness temperature in K, respectively. We take the gas temperatures to be the excitation temperatures, i.e., 16 K for the SMA2a outflow and 13 K for other outflows. The outflow mass, momentum, energy, dynamical age, and outflow rate are then given by

$$\begin{eqnarray*} M &=& d^2\left[\frac{\rm H_2}{\rm CO}\right]\,\overline{m}_{\rm H_2} \int _{\Omega }N_{\rm CO}(\Omega ^\prime)d\Omega ^\prime,\\ P &=& Mv,\\ E &=& \frac{1}{2}Mv^2,\\ t_{\rm dyn}&=&\frac{L_{\rm flow}}{v_{\rm max}},\\ \dot{M}_{\rm out} &=& \frac{M}{t_{\rm dyn}}, \end{eqnarray*}$$

where d is the source distance, $\overline{m}_{\rm H_2}$ is the mean mass per hydrogen molecule assumed to be 2.33 atomic units, Ω is the total solid angle that the flow subtends, v is the outflow velocity with respect to the systemic velocity of 78.4 km s⁻¹, v_max is the maximum outflow velocity, and L_flow is the flow length. We adopt an empirical $[\frac{\rm H_2}{\rm CO}]$ abundance ratio of 10⁴ (Blake et al. 1987).

Table 4 lists derived physical parameters for each outflow. These outflows show symmetry in blueshifted and redshifted lobes. The energetics are small compared to outflows observed in high-mass protostellar objects (Beuther et al. 2002c; Zhang et al. 2005), but are similar to outflows emanating from nearby (≲ 1 kpc) intermediate-mass hot cores (IMHCs: IRAS 22198+6336, Sánchez-Monge et al. 2010; NGC 7129-FIRS 2, Fuente et al. 2005; IC 1396 N, Neri et al. 2007). Particularly, the outflow rate $\dot{M}_{\rm out}$, which amounts to (4–47) × 10⁻⁶ M_☉ yr⁻¹, and the flow momentum P, which ranges from 1 to 9 M_☉ km s⁻¹, are of the same order as all the three IMHC outflows known to date. The dynamic timescales are (1.5–3.4) × 10⁴ yr for all the outflows. These values are comparable to massive outflows but are one order of magnitude higher than the IMHC outflows. Because observations are limited by noise, and gas at higher velocities is possible, t_dyn is a lower limit of the outflow age and the accretion history of the dominant protostar. Theoretical models suggest that the massive star formation process takes a few times the free-fall timescale (t_ff$=\sqrt{\vphantom{A^A}\smash{\hbox{${\frac{3\pi }{32G\rho }}$}}}=\frac{3.66\times 10^7 \rm {yr}}{\sqrt{n(\rm {cm}^{-3})}}$) of the core. Beuther et al. (2002c) observed a rough equality between t_dyn and t_ff. For the P1 cores, we obtain t_ff of (1.2–6.7) × 10⁴ yr based on the densities given in Table 3. We see that t_ff approximates t_dyn well, and this supports the idea that flow ages are good estimates of protostar lifetimes (Beuther et al. 2002c).

Table 4. Derived Outflow Parameters

Parametera	SMA1		SMA2a		SMA3		SMA4a		SMA5
	Blue	Red	Blue	Red	Blue	Red	Blue	Red	Blue	Red
v (km s−1)b	[67,76]	[85,93]	[41,76]	[84,116]	[50,71]	[85,99]	[64,72]	[85,107]	[66,73]	[86,93]
M (M☉)	0.20	0.07	0.38	0.23	0.22	0.18	0.06	0.12	0.05	0.06
P (M☉ km s−1)	1.12	0.71	5.09	4.04	3.17	2.05	0.52	1.50	0.35	0.63
E (M☉ km2 s−2)	3.85	3.62	53.46	46.18	26.96	13.06	2.57	11.15	1.40	3.40
Lflow (pc)c	0.40	0.27	0.51	0.48	0.56	0.51	0.44	0.43	0.34	0.39
tdyn (104 yr)	3.43	1.81	1.33	1.25	1.93	2.42	2.99	1.47	2.68	2.61
$\dot{M}_{\rm out}$ (10−5 M☉ yr−1)	0.58	0.41	2.86	1.84	1.14	0.74	0.20	0.82	0.19	0.23

Notes. ^aParameters are not corrected for inclination angle. ^bRange of velocities used to derive outflow parameters. ^cProjected length of the longer lobe in the case of two lobes.

Download table as:  ASCIITypeset image

The SMA2a, SMA3, and SMA4a outflows are so well collimated that they resemble "molecular jets" commonly found in low-mass star-forming regions (e.g., L1157, Zhang et al. 2000; OMC-1S, Zapata et al. 2005; HH 211, Lee et al. 2007). Recently, interferometric observations have revealed collimated outflows in high-mass star-forming regions (e.g., AFGL 5142, Hunter et al. 1999; Zhang et al. 2007; IRAS 20126+4104, Cesaroni et al. 1999; Shepherd et al. 2000; IRAS 05358+3543, 19410+2336, 19217+1651, Beuther et al. 2002a, 2003, 2004; HH80-81, Qiu & Zhang 2009). To measure the collimation, we take the ratio of outflow total length (blue + red lobes) over width measured in Figure 4 as the collimation factor. We obtain 14, 10, and 7 for SMA2a, SMA3, and SMA4a outflows, over scales of 0.5–1 pc. These values are comparable to those of the most collimated low-mass outflows and are higher than those of the high-mass outflows discovered by single dishes and previous interferometers (Wu et al. 2004; Beuther et al. 2002c; Qiu & Zhang 2009). Beuther et al. (2007) reported one outflow in IRAS 05358+3543 with a collimation factor of ∼10 and claimed that it was the first massive outflow observed with such a high degree of collimation at a scale of 1 pc. In G28.34-P1, the SMA2a and SMA3 outflows are among the best collimated molecular outflows yet discovered in high-mass star-forming regions over parsec scales. We note that these collimation factors should be considered as lower limits because the outflow lengths and widths are not corrected for inclination angle and are not deconvolved with synthesized beam.