FILAMENTARY ACCRETION FLOWS IN THE EMBEDDED SERPENS SOUTH PROTOCLUSTER

Understanding how clusters of stars form is a topic of considerable interest, since most stars appear to form within clustered environments (e.g., Lada & Lada 2003; Bressert et al. 2010). As noted in Myers (2011a), young stellar clusters are tightly packed with stars in their centers; the initial gas mass density inferred to create these stars is significantly higher than observed formation conditions. One simple solution to this problem is to feed in additional mass as the cluster is forming. Analysis of competitive accretion simulations indeed show that this can occur (Smith et al. 2009). Observationally, young embedded cluster-forming systems are often found to be associated with multiple filaments, in a hub-spoke like geometry (e.g., Myers 2009a). Recent Herschel observations have shown that filaments are ubiquitous in star-forming regions, and furthermore that prestellar cores and protostars tends to be associated only with those filaments that are gravitationally supercritical (André et al. 2010). Filaments therefore appear to be the prime candidate for providing a continuous source of material to be accreted onto forming clusters. Few observations exist, however, to test this idea, as most have focused on the column density structure of filaments and clusters (e.g., Arzoumanian et al. 2011).

Several decades ago, low resolution CO observations showed that some filaments have velocity gradients running along their length (on ∼10 pc scales, e.g., McCutcheon et al. 1982; Dobashi et al. 1992). Additionally, some filaments had hints of a velocity gradient across their short axis (e.g., Uchida et al. 1990; Dobashi et al. 1992). Higher resolution observations of dense gas species, however, have been lacking until recently. Hacar & Tafalla (2011) observed the low-mass star-forming region L1517 in Taurus, and showed it contains several coherent, roughly thermal filaments which are well-fit by an isothermal cylinder model. These filaments have small modulations in their (column) density and velocity profiles along the long axis, suggesting core accretion from the surrounding filament material. Hacar et al. (2013) conducted a similar study of the filaments in B213 in Taurus, finding evidence that the singular filament visible in column density maps is actually composed of many distinct velocity components of gas, each of which has a coherent and thermal velocity profile. In more massive filaments, spectral observations tend to indicate much more dynamic gas motions. Miettinen (2012) made pointed observations of C¹⁸O and C¹⁷O dense clumps embedded in the filamentary IRDC G304.74+01.32 and find strong signs of infall and strongly non-thermal line widths. If the non-thermal line widths can be interpreted as providing additional support, then the clumps and filament may be in virial equilibrium. In the DR21 filament, Schneider et al. (2010) find signatures of infalling motion over both the filament and clump scales, and kinematic signatures consistent with the filament forming from converging flows.

In order to study the question of filament accretion onto clusters, we present kinematic observations of the Serpens South embedded cluster and surrounding filaments. To the best of our knowledge, this is the first study to obtain sufficient complementary measurements of the gas velocity to provide reasonable certainty in attributing the motions to filamentary accretion.

The Serpens South cluster was first discovered in 2008 by Gutermuth et al. from Spitzer data. In this discovery paper, the cluster is clearly embedded in a dense filamentary structure which shows up in absorption at 8 μm. Gutermuth et al. (2008) identify approximately 90 young stellar objects (YSOs), and find that the ratio of class II to I sources in the region is low, around 0.7, which is much lower than the median value of 3.7 in nearby cluster-forming regions (Gutermuth et al. 2009), and suggests that the Serpens South embedded cluster is among the youngest examples known. As such, Serpens South may be a particularly good target to search for accretion flows in filaments and onto the central cluster, since there will have been less time for stellar feedback to disrupt or complicate the cluster-forming gas motions. Based on their Spitzer data, Gutermuth et al. (2008) estimate that stars are forming in the central cluster at a rate of 90 M_☉ Myr⁻¹, although more recent work by Maury et al. (2011) using a combination of Spitzer, Herschel, and MAMBO 1.2 mm dust continuum observations suggest that the rate may be several times smaller. Much of this discrepancy is due to the timescales assumed for protostellar evolution.

Additional recent observations of Serpens South include near-IR polarization measurements suggesting a relatively uniform magnetic field aligned perpendicular to the main filamentary structure (Sugitani et al. 2011), and a catalog of protostars associated with outflows suggesting that the energy injected in the central embedded cluster by the outflows may be enough to maintain the current level of non-thermal motion (Nakamura et al. 2011). X-ray observations from Chandra combined with the Spitzer data show that the YSOs have similar properties to those in Serpens Main and NGC 1333 (e.g., Winston et al. 2010), with analysis suggesting that the class III sources are likely a byproduct of early disk-stripping from inter-cluster dynamics, rather than indicating a population older than the class II sources (E. Winston et al., in preparation).

We obtained spectral line maps with Mopra⁵ across ∼100 arcmin² of the main filament in Serpens South in a variety of molecular tracers (discussed in the next section) in order to measure the motion of the dense gas in the filaments. We combine this data with measurements of the column density distribution (R. A. Gutermuth et al., in preparation) obtained with ASTE using the AzTEC camera (Wilson et al. 2008) to fully interpret the dynamical state of the filaments. Using these data, we find strong evidence that filamentary accretion flows are playing an important role in the ongoing evolution of the Serpens South cluster.

In Section 2, we discuss our Mopra observations, and the basic filamentary structure seen in Section 3. In Sections 4 and 5, we discuss analyze the spectral data and estimate the filamentary accretion rates. In Section 6, we perform some simple analysis of the column density structure of the filament, then discuss our main results in Section 7, and conclude in Section 8.

2.1. Data

Our observations were made with the ATNF Mopra 22 m telescope in Australia, in 2008 September, using the single pixel broadband spectrometer, MOPS.⁶ MOPS allows for multiple spectral windows to be observed simultaneously, increasing the mapping efficiency for multiple spectral lines; the spectral resolution in the zoom mode is ∼0.108 km s⁻¹ for the 90 GHz receiver. We used eight zoom windows to focus on different molecules of interest—the (1–0) transitions of N₂H⁺, HNC, HCN, H¹³CN, HCO⁺, H¹³CO⁺, plus SiO(2–1) and NH₂D(1,1). Most of these transitions trace dense gas; HCO⁺ is sensitive to outflows, while SiO traces primarily shocked emission (see, e.g., Sanhueza et al. 2012; Bergin & Tafalla 2007; Vasyunina et al. 2011 for a discussion on the physical conditions traced by these molecules).

We made on-the-fly maps, scanning along R.A., and decl., with 11'' spacing between scan rows, and a scan speed of 371 s⁻¹. Square maps of 53 per side (with the inner 5' fully sampled) were made and later placed on a single grid using Gridzilla. An OFF position of (18^h16^m56^s, +03°02'19'') was observed roughly every three minutes, and a measurement of the system temperature taken every thirteen minutes. Every map required roughly an hour of integration time for each of the R.A. and decl. scans. We used Livedata⁷ to apply the bandpass calibration and remove a low order polynomial baseline. We then created uniformly gridded spectral cubes for each waveband using Gridzilla,⁷ using a final pixel size of 15''; Mopra's spatial resolution is ∼40'' at 3 mm (Ladd et al. 2005). We converted the spectra from $T_a^*$ to T_mb units using an antenna efficiency of 0.49 (Ladd et al. 2005). Due to the relatively large separation between the adopted OFF position and our mapping area, with this initial reduction, we found some spectra still showed significant ripples in their baseline. In these cases, CLASS⁸ was used to fit and remove a higher order polynomial baseline.

Figure 1 shows the region we mapped with Mopra, compared to a dust emission map recently obtained with AzTEC/ASTE by R. A. Gutermuth et al. (in preparation). In Figure 1, the grayscale image and contours show the dust emission, while the blue contour shows the area we mapped with Mopra. The yellow diamonds, orange open circles, and red filled circles show the YSOs identified by Gutermuth in a re-analysis of the data in Gutermuth et al. (2008) which now includes Spitzer 24 μm data, and adopts the classification scheme of Gutermuth et al. (2009). Clearly, the Mopra observations span the main filamentary structure in Serpens South. In total, we have more than 1500 spectra at each of the 8 wavebands, and cover an area of roughly 100 arcmin². The median value of the rms noise per channel in the spectra are ∼0.34 K km s⁻¹ for all eight molecular species. Table 1 summarizes the central frequency adopted for each line and the associated reference. For lines with hyperfine structure, the transition associated with the adopted frequency is also listed.

Zoom In Zoom Out Reset image size

Table 1. Basic Line Properties

	Frequency	HF Componenta	Ref.b	Int. Rangec	Integrated Noised
Molecule	(MHz)	(if applicable)		(km s−1)	(K km s−1)
N2H+(1–0)	93173.4669	F1, F = 2–1, 2–1	1	−10.5→7.5	4.38
H13CN(1–0)	86340.1840	F = 2–1	2	−8.5→6.5	3.87
HCN(1–0)	88631.8473	F = 2–1	3	−9→7	4.02
HNC(1–0)	90663.5560	F = 2–1	4	−2→2	2.04
H13CO+(1–0)	86754.2884	N/A	5e	−2→2	2.02
HCO+(1–0)	89188.5260	N/A	6	−5→5	3.23
SiO(2–1)	86846.9600	N/A	7	−1.5→1.5	1.76
NH2D(1,1)	85926.2630	N/A	8	−6.5→6.5	3.66

Notes. ^aIn the case of a species with hyperfine splitting, the transition used for centering. ^bReference for the adopted frequency: (1) Pagani et al. 2009; (2) Pearson et al. 1976, as recommended by JPL (http://spec.jpl.nasa.gov); (3) Ahrens et al. 2002, as recommended by CDMS (Müller et al. 2005); (4) Bechtel et al. 2006; (5) Schmid-Burgk et al. 2004, as recommended by CDMS; (6) Ulich & Haas 1976, as recommended by NIST (Lovas et al. 2009); (7) Manson et al. 1977, as recommended by CDMS; (8) Main frequency listed for no hyperfine splitting from Bester et al. (1983), as recommended by NIST (our sensitivity is too low to pick up the satellite components). ^cThe range in velocities, relative to the line centroid, used when calculating the integrated intensity. ^dThe median formal 1σ error on the integrated intensity in brightness temperature units, i.e., the typical spectral rms multiplied by the square root of the number of spectral channels integrated over. ^eNote that hyperfine splitting was detected by Schmid-Burgk et al. (2004), but with the components separated by only 0.133 km s⁻¹, this would not be visible in our Mopra data.

Download table as:  ASCIITypeset image

2.2. Line Fitting

For the emission lines without hyperfine structure (HCO⁺, H¹³CO⁺, and SiO),⁹ we fit each spectrum with a single Gaussian, using the least squares minimization code MPFIT (Markwardt 2009; see also http://purl.com/net/mpfit) to constrain the fit parameters to a reasonable range and calculate the errors associated with each fitting parameter. We similarly fit NH₂D with a single Gaussian, as the line was sufficiently weak that none of the hyperfine satellites were detected (in most cases, not even the central component was visible). For these four transitions, the parameters fit are the velocity centroid, velocity dispersion (in Gaussian sigma units), and line peak.

For the emission lines with hyperfine structure (N₂H⁺, HNC, HCN, and H¹³CN), we substitute the Gaussian model with a single uniform slab model with the appropriate hyperfine emission information (adapted from Rosolowsky et al. 2008), and again use MPFIT to find the best fitting line parameters—in this case, the velocity centroid, velocity dispersion (in Gaussian sigma units), optical depth, and excitation temperature.

For all of the results discussed in this paper, we only analyze spectra which are above a signal to noise level of 3.5—i.e., we require an integrated intensity to be 3.5 times that expected from the noise level (the spectral rms times the square root of the number of velocity channels), and also a minimum of two spectral channels showing emission about 3.5 times the spectral rms (using two channels instead of one decreases the likelihood of the emission being a noise spike). Calculation of the integrated intensity is discussed in Section 3.

Over nearly the entire area mapped, most of the emission lines detected show strong indications of self-absorption; only N₂H⁺, H¹³CO⁺, and, where detected, NH₂D appear single-peaked, while SiO has no unambiguous detections. Even the N₂H⁺ lines may be mildly self-absorbed in some locations, but our signal to noise and spectral resolution are insufficient to determine this definitively. Figure 2 shows the spectrum of all eight emission lines near the cluster center, at the location where the N₂H⁺ emission is strongest, (−18^h30^m01^s, −2°03'12''). In all cases, the spectra have been summed over the nearest 4 × 4 pixels (60'') to increase the signal to noise. The red vertical lines indicate the location of the lines based on the centroid velocity fit to the summed N₂H⁺ spectrum. Hyperfine components are indicated, along with their relative intensities expected for the optically thin case. As the figure clearly illustrates, many of the emission lines have high levels of self-absorption. We therefore focus primarily on the line fit parameters for the N₂H⁺ and H¹³CO⁺ emission. Note that since there are seven N₂H⁺ hyperfine components with different optical depths, the N₂H⁺ line widths and centroids fit tend to be satisfactory, even if a small amount of self-absorption is present.

Zoom In Zoom Out Reset image size

We first calculate the integrated intensity for the five molecular lines that do not exhibit self-absorption and hence are more likely to be optically thin to show the column density structure of the dense gas. The full MOPS bandwidths are very large (∼200 km s⁻¹); in order to minimize the contribution of noise to the calculation, we restrict our integration to velocity channels close to where emission is expected. The velocity range is sufficient to extend slightly beyond the maximum observed line extent; the range adopted for each emission line is given in Table 1. The velocity center for integration is, in order of preference, as the signal to noise level allows: the centroid velocity fit to the particular line, the centroid velocity fit to the N₂H⁺ line,¹⁰ or the approximate cloud-wide centroid velocity, 7.5 km s⁻¹. Figure 3 shows the integrated intensity maps for N₂H⁺, H¹³CO⁺, and NH₂D, i.e., the brightest species without self-absorption, along with SiO, which lacks any clear detections. The typical uncertainty is given in Table 1.

Zoom In Zoom Out Reset image size

Figure 3 shows that both the N₂H⁺ and H¹³CO⁺ integrated intensities have a similar structure to that seen in dust emission. The integrated intensity for NH₂D is much fainter, but appears to follow a similar trend. All three of these molecular species are expected to be good tracers of dense gas (see, e.g., the discussion in Sanhueza et al. 2012; Bergin & Tafalla 2007), which is confirmed by their good correspondence to the dust emission. H¹³CN (not shown in Figure 3) is also expected to be a good dense gas tracer (e.g., Hily-Blant et al. 2010); however, our data has too poor a signal to noise level to see this definitively. Summing our data over many pixels shows hints that the central filament is brighter than the surrounding gas, although even this is only seen at a marginal level. The H¹³CN spectrum shown in Figure 2 is one of the brightest spectra resulting from summing the data over 4 × 4 pixels squared.

SiO is not detected with certainty anywhere in the map. As a tracer of shocked gas (see, e.g., Caselli et al. 1997 and references therein), its emission structure could be expected to differ from the other four species. The general lack of emission may indicate a lack of strong (highly supersonic) shocks; however, with the noise levels in our observations, we cannot rule out weaker shocks.

3.1. Filament Geometry

In Figures 1 and 3, differing filament shapes or geometries are visible north and south of the main cluster. The southern filament is dominated by a single nearly straight filament extending to the southeast, with several fainter filaments (e.g., pointing due south) apparent primarily in the dust emission map. In the north, however, the filamentary structure is more complex. A careful visual examination of both the dust emission and molecular data to the northwest of the main cluster shows a seemingly complex structure suggestive of overlap along the projected line of sight. Intertwined filaments have recently been identified through a careful analysis of the quiescent B213 filament in Taurus (Hacar et al. 2013). Other observations of filaments, such as Orion (Johnstone & Bally 1999; Bally et al. 1987), show morphologies suggestive of similar phenomena in higher-mass systems. Unlike B213, the large line widths in the northern filament of Serpens South make the separation of velocity components difficult if not impossible, and for this paper, we restrict our analysis to the simpler southern filament.

We define the central or peak ridge of the southern filament by searching for the location of peak emission along every horizontal cut through the filament. For this detection, we use the dust continuum map; while the N₂H⁺ integrated intensity map shows a similar distribution, the dust map has a higher signal to noise level. This filament definition is shown in Figure 4. For comparison, we also show a similar outline for the more complex northern filaments. For our analysis of the southern filament, we consider the data within 1' of the peak ridge (i.e., three 15'' pixels on either side of the peak, or a diameter of 0.16 pc), and along a 0.33 pc extent. This area covers all of the strongest emission present in the southern filament and excludes any gas potentially associated with the much fainter filament to the west. Within this area, the southern filament has a mass of ∼20 M_☉, based on the dust emission map.¹¹ Maury et al. (2011) estimate a mass of 610 M_☉ over a 15' by 25' area (i.e., several times larger than the entire extent of our Mopra observations), based on 1.2 mm MAMBO dust continuum data. For both our estimate and that in Maury et al. (2011), a distance of 260 pc (e.g., Straižys et al. 2003), is assumed, based on extinction estimates in the Aquila rift.

Zoom In Zoom Out Reset image size

There is, however, some debate over the true distance to Serpens South. Recent very long baseline interferometry (VLBI) measurements of a YSO in the Serpens main core suggest that the distance to Serpens main is roughly 415 pc (Dzib et al. 2010). Complicating the picture is the fact that the gas velocities appear similar for Serpens main, Serpens South, and Aquila; see the discussion in Maury et al. (2011) for more details. Since most existing analyses of Serpens South assume the closer distance of 260 pc (e.g., Gutermuth et al. 2008; Maury et al. 2011; Arzoumanian et al. 2011), we also adopt that value here, but discuss in Section 7 how our results would differ if we instead assumed 415 pc.

Figure 4 shows the centroid velocity fit to our N₂H⁺ spectra (where the signal to noise ratio is 3.5 or more). A clear velocity gradient is visible along the southern filament (shown by the red dashed line), while the northern filaments have a structured, but more complex velocity field. If the southern filament is oriented with its southeast part further from us and the cluster center closer to us, then an accelerating accretion flow would produce a velocity gradient similar to the one observed. The infall line profiles discussed in the next section lend further weight to this choice of filament orientation.

We measure the velocity gradient for the southern filament following the procedure outlined in Goodman et al. (1993). We assume the N₂H⁺ centroid velocities follow a simple linear form:

$$\begin{equation} A \Delta \alpha + B \Delta \delta = C, \end{equation} \tag{ 1 }$$

where Δα and Δδ are the changes in right ascension and declination across the map. We found the best-fitting values to the constants A and B using MPFIT (Markwardt 2009), and then converted all of the angular velocity gradients to physical scale assuming a distance of 260 pc. This yields an overall velocity gradient of 1.4 ± 0.2 km s⁻¹ pc⁻¹, at an angle of ∼20° west of north. Note that while the H¹³CO⁺ centroid velocities (not shown) agree quite well with the N₂H⁺ centroid velocities, fewer of the spectra have a signal to noise level above 3.5, which leads to a less reliable measurement of the velocity gradient.

We can estimate the accretion rate implied by this velocity gradient by using a simple cylindrical model, shown in Figure 5. The cylinder has a mass M, length L, radius r, an inclination to the plane of the sky of angle α, and true motions of velocities V_∥ along the filament long axis, and V_⊥ infalling into the short axis (this final term will be of interest in the following section). The accretion rate, $\dot{M}_{\parallel }$, along the filament onto the central cluster is therefore the velocity along the filament times the density and the area perpendicular to the flow, or

$$\begin{equation} \dot{M}_{\parallel } = V_{\parallel } \times \left (\frac{M}{\pi r^2 L}\right) (\pi r^2), \end{equation} \tag{ 2 }$$

which simplifies to $\dot{M}_{\parallel } = V_{\parallel } ({M}/{L})$. Due to projection effects, we observe

$$\begin{equation} L_{{\rm obs}} = L \cos (\alpha) \quad {\rm and} \quad V_{\parallel ,{\rm obs}} = V_{\parallel }\ \sin (\alpha) \end{equation} \tag{ 3 }$$

with V_{∥, obs} = ∇V_{∥, obs}L_obs. Therefore, the accretion rate is given by

$$\begin{equation} \dot{M}_{\parallel } = \frac{\nabla V_{\parallel ,{\rm obs}} M}{\tan (\alpha)}. \end{equation} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

For the southern filament, we estimate a mass of 20.3 M_☉ along a length of 0.33 pc and width of 0.08 pc (see Section 3.1), which corresponds to an accretion rate of 28 M_☉ Myr⁻¹ for tan (α) = 1. Note that since we observe both a gradient along the filament and infall onto the filament (discussed in the next section), this already constrains α not to lie close to either of the possible extreme values (0° for lying parallel to the plane of the sky and 90° for lying directly along the line of sight). We will furthermore show that the accretion rates estimated from both the motion along and across the filament are relatively high, which suggests that α is likely close enough to 45°, so that the correction factor for each measure is only several times unity.

The velocity gradient used for the accretion rate estimate is also similar to values measured in other star-forming regions—for example, Bally et al. (1987) find a gradient of ∼0.7 km s⁻¹ pc⁻¹ along several pc of the Orion integral shaped filament; Schneider et al. (2010) find gradients of ∼0.8–2.3 km s⁻¹ pc⁻¹ along the DR21 filament, and Kirk et al. (2010) find a general relationship in non-filamentary molecular gas structures between size and velocity gradient suggesting values around 1 km s⁻¹ pc⁻¹ for size scales of a few tenths of a parsec. The similarity of the measured gradient further suggests that projection effects are not extreme for the filament, if Serpens South is not atypical.

If the gradient were instead interpreted as rotation, then the total rotational energy would be

$$\begin{equation} E_{{\rm rot}} = \frac{1}{2} \left(\frac{M L^2}{3}\right) (\nabla V)^2, \end{equation} \tag{ 5 }$$

where the second term is the moment of inertia for a cylinder rotating longways around its end, and the final term represents the angular velocity squared. For the southern filament, this yields E_rot ≃ 0.7 M_☉ km² s⁻². The gravitational binding energy is given by

$$\begin{equation} E_{{\rm grav}} = \frac{G M M_{{\rm clust}}}{L}, \end{equation} \tag{ 6 }$$

where M_clust is the mass of the central cluster, roughly 100 M_☉, calculated in the same manner as the filament mass (see Section 3.1). Thus, the gravitational binding energy is a factor of more than ten higher, at E_grav ≃ 26 M_☉ km² s⁻². With a free-fall time of just ∼0.3 Myr, several times less than the estimated age of the region, even if the filament was initially rotating, infall motions would rapidly dominate.

Finally, the velocity gradient could also be interpreted as outflow motion, if the orientation was the reverse of that assumed above. In this instance, stellar winds and outflows would be the most likely candidate source for driving the motions. Using a simple calculation to estimate the momentum and energy associated with such outflowing motions yields values of 9 M_☉ km s⁻¹ and 2 M_☉ km² s⁻² respectively. Nakamura et al. (2011) estimate that the total momentum and energy in outflows (measured using CO(3–2) data) in Serpens South is ∼8 M_☉ km s⁻¹ and ∼65 M_☉ km² s⁻¹ respectively, with a total mass of outflowing material of only 0.6 M_☉. While in principle, the YSO population could provide sufficient energy and momentum to drive the observed large-scale flow, the outflows observed by Nakamura et al. (2011) are multi-directional, not highly collimated, as would be required.

We next utilize information from an optically thick line showing self-absorption to measure gas accretion across (i.e., building up additional mass onto) the filament. In the classic case of two converging layers, the gas is assumed to have an increasing excitation temperature toward the center (where the density is also higher); self-absorption in optically thick lines occurs for the emission originating from the furthest pieces of gas at any given velocity. Therefore, emission from the (blueshifted) layer moving toward the observer is emitted primarily from the higher excitation, denser part of the slab, while the reverse is true for the (redshifted) layer moving away from the observer (e.g., Myers et al. 1996). Infall profiles therefore show brighter peaks blueward of the peak emission of an optically thin line, with a lower red shoulder, and sometimes a central dip. Using simplified assumptions about the gas properties, the emission profile of a line showing such a red-blue asymmetry can be used to estimate a characteristic infall velocity of the gas (e.g., Evans 1999; Myers et al. 2000).

In Serpens South, we observe self-absorption in several of our emission lines: HCO⁺, HCN, and HNC. Of these lines, the HNC is the least self-absorbed, and therefore is the brightest and has the highest signal to noise level. A high signal to noise level is required to obtain reasonable model fits to infall spectra—deVries & Myers (2005) find signal to noise levels >30 are required for good results. In order to further maximize our signal to noise level, we therefore spatially average spectra along the southern filament. As discussed in the previous section, there is a global velocity gradient running along the southern filament, which needs to be accounted for prior to averaging. We use the N₂H⁺ centroid velocities fit at each pixel to first shift all of the HNC spectra to a common velocity center prior to averaging, again considering spectra within 1' (or 3 pixels on either side) of the filament's peak ridge identified. Figure 6 shows the result of this summation for HNC, as well as for the two optically thin lines, N₂H⁺ and H¹³CO⁺. In the case of N₂H⁺, only the isolated component is shown for simplicity. The HNC spectrum shows a brighter blue and fainter red peak with strong central dip around the peak emission of N₂H⁺ and H¹³CO⁺, in agreement with expectations for infall. The HNC(1–0) transition does have hyperfine structure; however, this is on a much smaller scale than the infall profile seen, as shown in the right hand panel of Figure 6.

Zoom In Zoom Out Reset image size

We estimate the infall velocity using the Hill5 model (and line fitting code) described in deVries & Myers (2005). The Hill5 model considers gas with a linearly increasing excitation temperature toward the center, with each half of the gas moving toward the other at a fixed velocity. In their analysis, deVries & Myers (2005) found that the Hill5 model tended to provide the best fit to a range of both simulated and real spectra with infall motions; the inclusion of an increase in the excitation temperature allowed for much better match than the classic two-slab model of Myers et al. (1996). The variables in the Hill5 model are the infall velocity, V_in, the (Gaussian sigma) velocity dispersion, δV, total optical depth, τ, system velocity V_lsr, and peak excitation temperature, T_peak. Note that the system velocity is expected to be 0 in our case, since the individual HNC spectra were all shifted by the N₂H⁺ centroid velocity fit at the same pixel prior to the summation. Using the Hill5 model, we find the best-fitting parameters listed in Table 2 under "Fit-1."

In order to estimate the fitting errors, we ran the fitter on 1000 realizations of the spectrum with random noise added at the same level as observed.¹² This test yielded several unexpected results. In a small number of cases (26 of 1000), the "best" fit clearly got "stuck" in a local minimum of χ² space, with unphysical parameters (optical depths in excess of 20 and peak excitation temperatures often around several hundred). These models provided visually poor fits to the spectra, and so we discarded these cases. These fits likely could have been fixed by modifying the fitting tolerances in the Hill5 fitter, but this seemed unnecessary, given the small overall number. Of the remaining 974 fits to the data plus random noise, 686 were found to have values clustered very close to the original fit (these were used to calculate the formal errors in "Fit-1"), while the remaining 288 were clustered around a different set of fit parameters, listed under "Fit-2" in Table 2. Visually, this other set of fits also are a reasonable match to the data, and the distribution of χ² values are similar for both sets of fits. Therefore, these second fit parameters appear to be a viable second solution to the best-fitting Hill5 model, using only the information available from the HNC spectrum.

We use the additional information available from the optically thin spectra (H¹³CO⁺ and possibly N₂H⁺) to argue that Fit-1 best matches the full suite of data. We calculate the corresponding optically thin spectra for both Fit-1 and Fit-2, using the fitted parameters but a lower optical depth. For Fit-2, the optically thin spectrum predicted has two distinct velocity components separated by more than the line width. Effectively, Fit-2 represents a model of two independent layers of material moving toward each other. As shown in Figure 6, two distinct velocity components are not observed in either of the brightest optically thin tracers, i.e., the N₂H⁺ isolated hyperfine component or H¹³CO⁺. This implies that Fit-1 is the most consistent with the observations.

We ran several simple tests with the online non-LTE radiative transfer calculator, Radex¹³ (Van der Tak et al. 2007) to check that the fitted values are reasonable. We assume an HNC abundance of 6 × 10⁻¹⁰ typical of dense core gas (Hily-Blant et al. 2010), and adopt gas properties (density, column density, temperature, total velocity dispersion) based on the observations discussed earlier. The predicted optical depth ranges from roughly 6 to 8, for temperatures (and the associated derived densities and column densities) between 12 and 15 K. Given that a cold dense core HNC abundance is likely an overestimate for a filament, and the uncertainties in various parameters used to estimate the density and column density from the dust maps, the Hill5 model fit of an optical depth of 2.7 for "Fit-1" is in reasonable agreement; the optical depth of 1.6 from "Fit-2" is less consistent with the observational expectations.

Using the same toy model as shown in Figure 5, the accretion rate onto the filament can be estimated from the above infall velocity, i.e., the velocity onto the filament times the density and filament surface area of the short axes, or

$$\begin{equation} \dot{M}_\perp = V_{\perp } \times \left(\frac{M}{\pi r^2 L}\right) (2\pi r L), \end{equation} \tag{ 7 }$$

which simplifies to $\dot{M}_\perp = 2 V_{\perp } ({M}/{r})$. Due to projection effects, we observe

$$\begin{equation} V_{{\rm in}} = V_{\perp } \cos (\alpha), \end{equation} \tag{ 8 }$$

while the observed filament radius/short axis should be invariant under rotation. Therefore, the accretion rate onto the filament is given by

$$\begin{equation} \dot{M}_\perp = \frac{2 V_{{\rm in}} M}{r \cos (\alpha)}. \end{equation} \tag{ 9 }$$

This corresponds to an accretion rate of 130 M_☉ Myr⁻¹ for cos (α) = 1 and is therefore a lower limit to the true rate. Additionally, the infall velocity fit with the Hill5 model is expected to be a lower limit to the true characteristic infall speed: the infall was modeled as two infalling parallel slabs, while in a cylindrical geometry, some of the inward motions will appear smaller along the line of sight. Finally, note that using the alternative V_in of 0.54 km s⁻¹ from Fit-2 corresponds instead to a lower limit of ∼280 M_☉ Myr⁻¹.

Both the accretion rate onto the filament and from the filament onto the cluster, without accounting for projection effects, are within a factor of a few of the star formation rate in the cluster estimated by Gutermuth et al. (2009), and somewhat larger than the rate estimated by Maury et al. (2011). Any correction for projection effects made to the accretion onto the filament will serve to further increase that estimate, i.e., beyond 130 M_☉ Myr⁻¹, while the accretion rate from the filament onto the cluster could be larger or smaller than our estimated 28 M_☉ Myr⁻¹ after correcting for projection.

If the filament were projected at an angle of 30° to 60° relative to the plane of the sky, then the accretion onto the central cluster measured would be within a factor of 0.6 to 1.7 of its true value, while the accretion onto the filament measured would be a factor of 1.15 to 2 of its true value. Based on the combination of filament appearance, velocities measured, and accretion rates derived, we believe that 30° ⩽ α ⩽ 60° is a reasonable range of filament angles.

5.1. Other Lines

In addition to the self-absorbed HNC emission, both HCO⁺ and all three hyperfine components of HCN also show signs of self-absorption. We performed the same spectral averaging across the southern filament for these emission lines as described for HNC. We additionally shift each HCN (and associated H¹³CN) hyperfine component to the velocity it would have if it were the sole emission line. Figure 7 shows the resulting spectral line profiles. As can be easily seen from the figure, all four of the optically thick lines show such strong self-absorption that relatively little emission is seen. We attempted to model the emission lines with a Hill5 (or other) profile, but due to the low signal to noise, even with averaging over the filament, we were unsuccessful.

Zoom In Zoom Out Reset image size

The degree of self-absorption in HCO⁺ and HCN is surprisingly large—the HCO⁺ emission peaks at only half the value of the H¹³CO⁺ peak, whereas in an optically thin environment, the ratio would be somewhere around a factor of 40 (for non-PDR molecular clouds at a galactocentric separation, similar to the Sun, of 8 kpc; Savage et al. 2002). In H¹³CN, the line ratios for optically thin gas would be 1:5:3 (for the top right, bottom left, and bottom right panels in Figure 7), whereas HCN is clearly brightest in the (F = 0–1) or first/leftmost hyperfine component. This would be consistent with a picture wherein the highest levels of self-absorption occur in the optically thickest components, and the least in the optically thinnest.

Recent work suggests that the rate coefficients for HCN and HNC differ significantly (Dumouchel et al. 2010), implying sensitivity to different density ranges of material, and thus explaining why the two do not show similar levels of self-absorption. Alternatively, studies have shown that hyperfine anomalies in HCN(1–0) are relatively common in higher column density environments, and tend to lead to precisely the trend in peak brightnesses observed (Loughnane et al. 2012). The H¹³CN emission is extremely faint, but appears to be roughly consistent with optically thin emission with the predicted line strengths (with no significant emission detected in the first/leftmost/F = 0–1 component); this is also consistent with the hyperfine anomaly study of Loughnane et al. (2012), who find that H¹³CN does not have anomalous ratios in sources where HCN does.

We ran several additional checks to ensure that these severely self-absorbed spectra are real, and were not affected by any systematic effects. For the velocity centroids, all eight molecules were observed simultaneously with MOPS, suggesting that any error with either the measured line velocities, or the adopted central frequencies of each molecule, would have to be present in all of our observations. The Serpens South data presented here are part of a larger project which includes four other embedded cluster-forming regions. Examining our data in the other four regions shows very good agreement between the different molecular line centroid velocities where there is optically thin emission, suggesting that there are no systematic effects with our velocity scale. Similarly, the peak line emission in these other regions, and the relative values of peak intensity, agree well with expectations, suggesting no systematic effects there either. Finally, we note that similarly severely self-absorbed emission line profiles showing a complete absence of a red shoulder have also been observed in Serpens South with a different telescope and molecular line transition. CO(3–2) emission observed with the ASTE 10 m telescope show a similarly shaped profile in some areas of the cluster (see, e.g., the bottom panels in Figure 6 of Nakamura et al. 2011).

6.1. Mass per Unit Length

In this section, we examine some basic properties of the column density structure of the southern filament. The simplest property which can be measured for filaments is their mass per unit length. Assuming that a filament can be approximated as an infinitely extended, self-gravitating, isothermal cylinder with no magnetic support, then the maximum mass per unit length which can be in equilibrium is

$$\begin{equation} M_{{\rm line,crit}} = 2 c_s^2 \left /\right. G \end{equation} \tag{ 10 }$$

$$\begin{equation} M_{{\rm line,crit}} = 16.7 \left(\frac{T}{10 \,{\rm K}}\right)\, M_{\odot }\, {\rm pc}^{-1} \end{equation} \tag{ 11 }$$

(e.g., Ostriker 1964; Inutsuka & Miyama 1997). Observations of the ammonia (1, 1) and (2, 2) inversion transitions (R. Friesen et al., in preparation) indicate that the kinetic temperature across Serpens South is in the range of 10 to 15 K, implying M_{line, crit} is between 17 and 25 M_☉ pc⁻¹. If the observed non-thermal motions were also considered to contribute to supporting the filament (e.g., Fiege & Pudritz 2000a; Miettinen 2012), then with a typical velocity dispersion of ∼0.36 km s⁻¹ (Gaussian sigma units) for N₂H⁺ (mean or median value over the southern filament), then M_{line, crit} is roughly 80 M_☉ pc⁻¹.

In the southern filament, we find a mass of 20.3 M_☉ along 0.33 pc, implying M_line = 62 M_☉ pc⁻¹, i.e., several times the critical value, unless all of the non-thermal motion is contributing thermal-like support. The total amount of support available is therefore likely insufficient to keep the filament in equilibrium, implying that it should be radially contracting. Supercriticality is consistent with our velocity measurements showing that gas is globally infalling in the filament. Arzoumanian et al. (2011) similarly find that filaments with higher column densities, such as in Aquila, tend to have supercritical values of M_line. Filaments associated with massive star formation tend to be highly supercritical with M_line often ten times larger than that in nearby regions (Hennemann et al. 2012).

6.2. Radial Column Density Profile

In this section, we analyze the radial column density profile of the filament using the millimeter continuum emission. We use the peak ridge of the filament discussed in Section 3, and for each pixel in the original AzTEC/ASTE data, we measure the column density profile along lines of constant declination. This provides for more regular and better sampling than to measure points perpendicular to the local filament orientation; we later correct our radial separation measurements to account for an average filament angle of 45°.

All ground-based (sub)millimeter measurements are insensitive to the largest scale structures, due to the observing methods which must be employed. For the AzTEC/ASTE map we analyze, R. A. Gutermuth et al. (in preparation) have devised a novel technique which is able to recover flux out to about 5', a much larger scale than is typical. In summary, the technique iteratively subtracts a model of the astronomical signal from the bolometer timestream data and then recomputes the atmospheric filtering and regenerates the signal model map. The initial model is just the base pipeline-produced map masked to a signal to noise of 2, and convergence of the process is achieved when no new signal is found in an iteration down to a signal to noise of 3. Despite this marked improvement in map reconstruction, the very largest-scale structures with shallow or flat radial column density profiles will not be included in the map. Some amount of large-scale structure is expected to be present, as Serpens South is coincident with a roughly degree wide overdensity of dust extending to the well-known W40 active star-forming region. The 2MASS-based dust extinction map of Rowles & Froebrich (2009)¹⁴ show that visual extinctions in excess of 5 to 10 mag are common across this area. Some, but not all, of this extinction undoubtedly represents diffuse material unassociated with the individual star-forming regions like Serpens South. Any material that is associated with Serpens South, however, would effectively result in a missing linear offset to the derived column density profile; the Serpens South filament is sufficiently compact that the shape of the column density profile should be correct. The function we fit to the column density profile is dependent on absolute values, not just the shape, so we limit our analysis to column densities where the fractional errors due to any missed large-scale structure are reasonable. If a constant background of A_V = 1 mag is missing from the AzTEC/ASTE map, there would be a ⩽30% underestimate of the inferred column density at fluxes of 0.067 Jy bm⁻¹; we adopt this as our lower limit. An additional complication to the analysis arises from a second, fainter filament to the west of the one we analyze. In order to prevent this structure from confusing the radial column density profile, we additionally exclude pixels which lie westward of the minimum value between the two filaments.

After measuring the radial column density profile along the entire peak ridge, we calculate the mean and standard deviation at each radial separation, similar to Arzoumanian et al. (2011). The mean profile and standard deviation of the column density profile of the southern filament is shown in Figure 8. For comparison, the dotted red line shows the ASTE beam; clearly the filament is resolved.

Zoom In Zoom Out Reset image size

We compare the filament column density profile to the same suite of models as in Arzoumanian et al. (2011):

$$\begin{equation} \Sigma _p(r) = A_p \frac{\rho _c R_{{\rm flat}}}{[1 + (r/R_{{\rm flat}})^2]^\frac{p-1}{2}}, \end{equation} \tag{ 12 }$$

where Σ_p is the column density at radial separation r, ρ_c is the central (mass) density, R_flat characterizes the inner flat part of the profile, and

$$\begin{equation} A_p = \frac{1}{\cos \, i} \int _{-\infty }^{\infty } \frac{du}{(1+u^2)^{p/2}} \end{equation} \tag{ 13 }$$

is a constant, and the inclination angle along the line of sight, i, is assumed to be 0 for simplicity.

Assuming that the observed column density profile is additionally convolved with the ASTE beam size, we use the MPFIT minimization code (Markwardt 2009) to determine the best fitting parameters. Allowing all three free parameters to vary, we find ρ_c = 3.6 ± 2.2 × 10⁵ cm⁻³, R_flat = 0.010 ± 0.006 pc, and p = 1.9 ± 0.2. Our values for R_flat and p lie within the range found by Arzoumanian et al. (2011) for all filaments within the Aquila Rift (note their ρ_c values were not published), including a power law index which is shallower than for an isothermal cylinder. It is important to note, however, that unlike Arzoumanian et al. (2011), we are not able to make a strong a statement on the power law index of the filament. The restrictions discussed above on where we are able to measure the radial column density profile limit our analysis to relatively close to the central ridge. Arzoumanian et al.'s comparable radial column density profile extends a factor of nearly 30 larger in radius than ours. With the relatively small range in radial separations, our ability to constrain the power law index is weak. Figure 8 shows also the best fit obtained using a power law index of 4 (corresponding to an isothermal cylinder), demonstrating that an isothermal model also can provide a reasonable fit to our data. The temperature required for the isothermal fit is 43 K, which is similar to the temperature implied by the non-thermal velocity dispersion in N₂H⁺ (around 50 K for σ = 0.36 km s⁻¹). This is consistent with Section 6.1 where we showed that thermal pressure was insufficient to balance the gravity.

Recently, Juvela et al. (2012) created synthetic observations of filaments formed in turbulent simulations, to investigate the effect of instrumental noise and changing dust properties on the inferred filament column density profile and mass per unit length, especially as relates to Herschel data. They find that the mass per unit length is reasonably robust to these uncertainties, while the parameters fit to the filament column density profile can be more uncertain in cases of higher noise/poorer resolution or a dust opacity which varies with density.

The southern filament in Serpens South has similar column density properties to star-forming filaments recently studied with Herschel (Arzoumanian et al. 2011), with a mass per unit length several times larger than the critical value for thermal support, and a radial column density profile which is consistent with p ∼ 1.5–2.5. In fact, the Serpens South filament is included in the Arzoumanian et al. (2011) sample of "Aquila" filaments, although the individual filament locations/properties are not given in their paper to allow a more detailed comparison. The relatively large velocity dispersions seen in dense gas tracers in the southern filament could be interpreted as evidence for providing sufficient stability against gravitational collapse, however, our more detailed kinematic analysis suggests that the filament is in fact not in equilibrium. We find evidence for mass flows both along the filament onto the central cluster and onto the filament from the less dense surroundings, with accretion rates of order 28 and 130 M_☉ Myr⁻¹ respectively. These rates are uncertain by a factor of a few, for several reasons, including the angle of the filament with respect to the plane of the sky.

Although both accretion rates are uncertain to a factor of a few, it appears likely that the accretion rate onto the filament is larger than the accretion rate from the filament onto the central cluster. Formally, in order to have the latter rate be higher, the filament angle would need to satisfy

$$\begin{equation} \frac{28\, M_{\odot }\, {\rm Myr}^{-1}}{\tan (\alpha)} \ge \frac{130\, M_{\odot }\, {\rm Myr}^{-1}}{\cos (\alpha)} \end{equation} \tag{ 14 }$$

which is satisfied for α ⩽ 12°, i.e., a filament which is nearly parallel to the plane of the sky. Such a small inclination angle seems unlikely, but note that uncertainties in all observed quantities (mass, radius, velocity gradient, and infall velocity) could allow for a slightly larger range of angles. One must always use caution in interpreting velocity measures from different molecular tracers. The gas velocity could be some (non-constant) function of density, and using different molecular tracers for each accretion rate measure could in principle reflect that, rather than overall flow differences in the two directions. The only way to circumvent this potential issue entirely would be to use observations of a single self-absorbed molecule and its rarer isotopologue for both measures. We do not have such data available at high enough signal to noise; however, the consistency of the HNC self-absorption dip with the N₂H⁺ and H¹³CO⁺ peak emission suggests all are tracing similar volumes of gas, and that this should not be a major problem.

If the tentative difference in accretion rates is true, this could imply that in the future, the accretion rate onto the central cluster will also increase (given the larger amount of mass available in the filament), and thus lead to more vigorous star formation in the cluster center. It is also possible that the increase in mass in the filament could promote local instabilities in the filament, leading to core formation, and eventually to a series of star formation sites along the filament. Dense cores are, in fact, already observed in the AzTEC/ASTE data along the filament, as well as in 8 μm absorption. A small number of class 0 and I YSOs are also seen, which may be an indication of past instabilities in the filament. The small cluster of class II YSOs just past the southern end of the southern filament suggests a previous episode of star formation, perhaps also associated with the filament.

At the current estimated accretion rate of 28 M_☉ Myr⁻¹, the cluster's membership would roughly double in 1 Myr with the addition of about 60 YSOs, assuming a mean mass of 0.5 M_☉. An increasing accretion rate would lead to a larger number of YSOs. From these estimates, we can see that our estimated accretion rate from the filament onto the cluster is at a comparable level to previous observational estimates of the rate of mass supply required to form the current generation of YSOs in the cluster center. We can therefore conclude that filamentary accretion is an important, if not the dominant, factor in governing the formation and evolution of stars in the central cluster.

Throughout this analysis, we have adopted a distance to Serpens South of 260 pc (e.g., Straižys et al. 2003), for consistency with previous work, although recent VLBI observations suggest that a distance of ∼415 pc (Dzib et al. 2010) applies to at least the main Serpens cluster. If this distance is also appropriate for Serpens South, then the quantities we estimate in our analysis here become a factor of several larger. The mass in the southern filament would increase by a factor of 2.5, the velocity gradient would decrease by a factor of 1.6, and all other quantities of interest (lengths, accretion rates, mass per unit length) would increase by a factor of 1.6; the accretion rates become 45 M_☉ Myr⁻¹ parallel to the filament (onto the cluster) and 208 M_☉ Myr⁻¹ perpendicular to the filament. None of these changes effect our conclusions qualitatively, and in fact, they tend to strengthen our finding that filamentary accretion is important in Serpens South.

Filament-cluster systems such as Serpens South present a unique opportunity to gain a deeper understanding of cluster formation. With more sensitive observations, it should be possible to measure the spatial variation in the infall velocity across the filament, which could be used to model the full three-dimensional velocity profile of material in the filament. This in turn can provide a straightforward test of the relative importance of magnetic fields, gravity, and turbulence in the ongoing cluster accretion, as discussed in Balsara et al. (2001) and Heitsch (2013).

Using a combination of dust column density information derived from AzTEC/ASTE continuum maps (R. A. Gutermuth et al., in preparation) and dense gas kinematics derived from Mopra spectral maps, we study the flow of gas along the southern filament of the Serpens South embedded cluster. We find that the filament appears typical in terms of column density properties (e.g., a supercritical mass per unit length). The multiple molecular emission lines observed in our Mopra survey allow us to infer the possible presence of an accretion flow from the filament onto the central cluster (if the filament is in front of the cluster), and additionally we find that material continues to be accreted onto the filament. We estimate accretion rates of 28 M_☉ Myr⁻¹ for accretion onto the central cluster, and a lower limit of 130 M_☉ Myr⁻¹ for infall onto the filament (assuming a distance to Serpens South of 260 pc). Both of these are similar in magnitude to the estimated current rate of star formation in the central cluster.

Filamentary accretion flows therefore appear to be an important mechanism for supplying the material necessary to form the YSOs in this cluster. As more embedded clusters and their associated filaments are identified in current surveys, it will be interesting to measure the associated gas kinematics to determine whether other embedded clusters behave similarly to Serpens South.

We thank the referee for suggestions which helped to clarify aspects of this paper. We thank Steve Longmore for assistance with the Mopra observing setup and data reduction. We thank Phil Edwards at ATNF for allocating additional time to complete our Mopra observations, and all the ATNF Narrabri staff for their support. We thank Rachel Friesen for sending her ammonia temperature results prior to publication, and for several interesting discussions on Serpens South. We thank Elaine Winston and Scott Wolk for sharing the results of their Chandra survey of Serpens South prior to publication. H.K. thanks Paola Caselli for a valuable discussion on astrochemistry, Chris de Vries for answering queries on his infall fitter on numerous occasions, and Fabian Heitsch for a critical read-through of an earlier draft of this paper, as well as sending a copy of his filamentary accretion paper pre-publication. H.K. acknowledges funding from a Smithsonian Scholarly Studies Program grant. This research was supported in part by the National Science Foundation under grant number 0708158 (T.L.B.).