Multicomponent Kinematics in a Massive Filamentary Infrared Dark Cloud

Vlas Sokolov; Ke Wang; Jaime E. Pineda; Paola Caselli; Jonathan D. Henshaw; Ashley T. Barnes; Jonathan C. Tan; Francesco Fontani; Izaskun Jiménez-Serra

doi:10.3847/1538-4357/aafaff

1. Introduction

The formation of stars from the diffuse ISM is a multistage process of gas channeling and mass assembly. In the well-established low-mass star formation picture, individual stars are formed from dense cores, which become gravitationally unstable once enough mass has been accumulated from their surroundings. While this process is thought to be relatively well understood, the formation of high-mass stars and star clusters occurs more rapidly and is a more dynamic process, requiring efficient gas channeling onto the protostellar sources, and its initial conditions are still not very well understood (e.g., Tan et al. 2014).

Observations of the earliest stages of massive star and cluster formation are greatly hindered by the high extinction of, and large distances to, high-mass star-forming regions. Thus, resolving the physical scales of the massive starless cores typically requires interferometric observations in radio or submillimeter wavelengths. The initial conditions of this process, such as density and temperature profiles on the core scales (∼0.1 pc), turbulence levels, accretion from the surrounding envelopes, magnetic field strengths, or the angular momentum of the core, are inherited from the more extended scales of clouds or protostellar clumps (∼1 pc) and are often poorly constrained by observations. As these conditions directly regulate the subsequent fragmentation and infall of the star-forming material (Motte et al. 2018; Tan 2018), high angular resolution observations of the dense infrared dark gas with little evidence of disruptive feedback can provide improved constraints on the theoretical models and numerical simulations aiming to re-create the formation of high-mass stars and star clusters.

Large-scale infrared surveys of the Galactic plane (Perault et al. 1996; Egan et al. 1998) have revealed infrared dark clouds (IRDCs), which are the most promising candidates to host the earliest stages of massive star and cluster formation (Rathborne et al. 2006). Often appearing as dense extinction ridges, tracing higher-density spines of their parent giant molecular clouds (GMCs; Hernandez & Tan 2015; Schneider et al. 2015), IRDCs appear as cold (10–15 K), filamentary, high column density (10²³–10²⁵ cm⁻²) and high-mass (10^3–5 M_⊙; Kainulainen & Tan 2013) extinction features, typically hosting a number of high-mass cores (Rathborne et al. 2006). IRDCs vary greatly in their masses, densities, and star-forming activity (Kauffmann & Pillai 2010). Rathborne et al. (2006) have selected 38 darkest clouds representative of the earliest stages of massive star formation from the Simon et al. (2006) sample, which contains over 10,000 IRDCs, and have shown that the clouds are further fragmented into high-mass cores with a median mass of 120 M_⊙. Butler & Tan (2009) studied mid-infrared extinction of 10 IRDCs from the Rathborne et al. (2006) sample, selecting the clouds that are nearby, are massive, and show relatively simple diffuse emission. Butler & Tan (2012) identified a number of compact (R ∼ 0.1–0.2 pc) dense cores embedded in these 10 IRDCs and discussed their star-forming potential.

In this study, we investigate the physical properties of dense gas connecting the cloud to core scales (from the IRDC size of ∼6 pc down to 0.08 pc) for 1 of 10 clouds in Butler & Tan (2009, 2012), IRDC G035.39–00.33. It is seen as a highly filamentary feature at mid-infrared wavelengths (e.g., Butler & Tan 2009) and belongs to a larger GMC complex (Hernandez & Tan 2015). By combining available near- and mid-IR extinction measurements, Kainulainen & Tan (2013) have estimated the mass of G035.39 to be 16,700 M_⊙. The cloud temperature is relatively cold (∼15 K from dust emission measurements; Nguyen Luong et al. 2011), with high CO depletion and large deuteration fractions from single-dish observations indicating that the cloud is governed by cold chemistry with little evidence of stellar feedback (Hernandez et al. 2011; Jiménez-Serra et al. 2014; Barnes et al. 2016). Jiménez-Serra et al. (2010) have discovered extended SiO (2−1) emission throughout the northern part of the cloud, suggesting that its origin may arise from a slow shock from the formation of the IRDC (a possibility later explored by Henshaw et al. 2013; Liu et al. 2018) or from the stellar feedback driven by embedded YSOs (suggested to be only partially responsible for the SiO emission; Nguyen Luong et al. 2011). Despite its starless appearance on large, parsec scales, Nguyen Luong et al. (2011) find a population of protostellar cores within G035.39, nine of which are massive enough (>20 M_⊙) to form intermediate- to high-mass stars. On smaller scales, Henshaw et al. (2014) have studied the kinematics of the dense gas in the northern part of the IRDC, resolving narrow dense filaments in it with a velocity pattern suggesting a dynamical interaction between the filaments and the dense cores in the cloud. Previously, Sokolov et al. (2017) probed the temperature structure of the IRDC Green Bank Telescope (GBT) observations and found evidence for a collective gas heating from the embedded protostars in G035.39, and Sokolov et al. (2018) presented the Very Large Array (VLA) findings of narrow-line widths in the IRDC G035.39 and discussed their implications for the high-mass star-forming regime.

In this paper, we present a detailed analysis of the combined VLA and GBT ammonia observations of the whole G035.39 cloud. The paper is structured as follows. In Section 2 we summarize the data reduction, discuss our strategy for VLA and GBT data set combination, and present an overview of its outcome. In Section 3 we describe the separation of the data into individual components and present the spectral line fitting results for the largest among them. In Section 4 a quantitative comparison with similar studies of G035.39 is made, the implications of cloud kinematics for velocity flows around the protostellar and continuum sources are given, and the stability of the individual components is discussed. Finally, a summary of our main results is given in Section 5.

2. Data Reduction

The Karl G. Jansky VLA observations (project AW776; PI: Ke Wang) in compact D configuration were carried out on 2010 May 8, mapping the NH₃ (1, 1) and (2, 2) inversion lines in the K band across the whole G035.39 cloud in a five-point mosaic. The two ammonia lines were mapped in two consecutive sessions, to cover the full hyperfine line structure of both transitions with the first-generation VLA correlator at a channel separation of 15.625 kHz. The effective spectral resolution is 18.75 kHz, which is the FWHM of the sinc spectral response function of the correlator. The data were calibrated on the quasars J1851+005 (gain), J2253+1608 (bandpass), and 3C 48 (flux) within the CASA data reduction package. The data presented have been previously described in Sokolov et al. (2018).

We have deconvolved the calibrated visibilities in CASA, using the tclean task, with the multiscale CLEAN algorithm (Cornwell 2008), following the Briggs weighting with the robust parameter set to 0.5. To achieve a similar synthesized beam for the two ammonia lines, we tapered the visibilities and applied a common restoring beam of 5 farcs 44. To recover the extended ammonia emission, we fill in the missing flux from the GBT data. A detailed description of the GBT data reduction can be found in Sokolov et al. (2017). Before merging the data, the GBT images were converted to spectral flux density units and convolved to the VLA channel separation of 0.2 from the original channel separation of 0.0386 . After regridding both GBT spectral cubes to match the VLA grid, we apply the VLA mosaic primary beam response to the GBT images. The VLA-only images for both NH₃ inversion lines were then combined with the processed GBT data using the feather task within the CASA package as an intermediate step, to obtain information about the spatial distribution of the extended ammonia emission. Finally, the VLA (1, 1) and (2, 2) images were deconvolved again with the tclean task, with the tclean mask determined from the aforementioned combined data. Constructing the clean mask from independently feathered images ensures that unbiased knowledge of the extended emission is incorporated into the tclean run. For specific details on the imaging strategy, we refer to Appendix A, where a full description of the imaging and data set combination is presented. The resulting NH₃ (1, 1) and (2, 2) spectral cubes, gridded into 1'' pixels, have a common restoring beam of 5 farcs 44. The typical rms value of the emission-free spectra in the resultant cubes is 14 mJy beam⁻¹ for NH₃ (1, 1) and 5 mJy beam⁻¹ for NH₃ (2, 2) inversion lines in a 0.2 channel (1.04 and 0.37 K, respectively).

2.1. Overview of the Data

Figure 1 presents the channel maps for NH₃ (1, 1) line, showing both the presence of multiple components throughout the IRDC and a line centroid change toward the redshifted regime as one follows the IRDC northward. The equivalent figure presenting the channel maps for the NH₃ (2, 2) line is shown in Figure 12 of Appendix B, showing similar morphology.

**Figure 1.** Channel maps of the ammonia (1, 1) line for the IRDC G035.39. Each panel shows the brightness temperature of the NH₃ (1, 1) spectral cube channel of 0.2 width, with text boxes at the bottom indicating the velocity centroid of the channel. The velocity range of 42–47 was chosen to capture the main hyperfine group of the inversion transition. The light-gray contours show the infrared extinction contours starting from A_V = 30 mag and progressing inward in steps of 30 mag (Kainulainen & Tan 2013).
Download figure:
Standard image High-resolution image

The parsec-scale ∼0.2 quasi-linear velocity gradient found in the GBT-only data (Sokolov et al. 2017) is present in our combined data: the IRDC kinematics reveal a gradual change into the blueshifted velocity regime as the filament is followed southward. Superimposed on this kinematics feature, additional velocity components are present. In the northern portion of the cloud, the western side of the IRDC splits off into a secondary, redshifted velocity component located around 47 , in agreement with the previous studies of the cloud (labeled F3 in Henshaw et al. 2013, 2014; Jiménez-Serra et al. 2014). A filamentary-like feature connects the northern region of the IRDC with the southern one, where it appears to split into two individual filaments that then join together at the southern part of G035.39. Additionally, the velocity feature between 42 and 43 , seen before as a clump-like structure in the ammonia maps of Sokolov et al. (2017), is now resolved into a distinct filament, orthogonally oriented to the rest of the IRDC.

2.2. Line Fitting

Molecular line emission tracing the dense gas within G035.39 is known to exhibit multiple line-of-sight velocity components (e.g., Henshaw et al. 2013, 2014; Jiménez-Serra et al. 2014). However, with no clear prior knowledge of the kinematical structure exhibited in the high-resolution ammonia emission, standard algorithms commonly used to fit the line profile are prone to converge away from the global minima or, alternatively, to diverge in the absence of a proper starting point of the algorithm. As our goal is to distinguish the physical properties of velocity-coherent structures, it is crucial to conduct a simultaneous fitting of all line-of-sight components, constraining their line parameters whenever a significant additional feature is present in the spectral profile. At the same time, overfitting has to be ruled out—that is, only statistically significant components have to be taken into consideration for the subsequent scientific analyses.

As the VLA spectral cubes consist of about 12,000 spectra with signal-to-noise ratio (S/N) higher than 3, the common practice of fitting multiple components by eye would require an excessive amount of human interaction, resulting in a low level of reproducibility, while methods for global nonlinear regression are prohibitively computationally expensive for such a high volume of data. Therefore, to fit multiple line-of-sight line components, we employ the same fitting strategy as in Sokolov et al. (2017), where an open-source Python package¹⁰ is used to assist the Levenberg–Marquardt fitter to find an optimal starting point for exploring the parameter space of the spectral model. The free parameters for the ammonia model (same formalism as in, e.g., Friesen et al. 2017), implemented within the pyspeckit package (Ginsburg & Mirocha 2011), are split into a multidimensional parameter grid, and a synthetic ammonia model is then computed for each. Subsequently, every spectrum in our VLA data is then cross-checked with every model computed to calculate the corresponding residual. The model parameters for a given spectrum that yield the best residual (i.e., the least sum of squared residuals) are taken to be the starting value for the Levenberg–Marquardt fitting algorithm.

Three independently derived sets of best-fitting parameters are obtained in such a fashion, for one-, two-, and three-component models.¹¹ The three fitting sets are then merged into one following a set of heuristical criteria that determine the number of line-of-sight components in a given pixel of the spectral cube. The heuristical procedure is described below:

1.
First, the three-component model is considered to be valid.
2.
If any of the three components are below S/N of 3—for both NH₃ (1, 1) and (2, 2) transitions—the model is marked as rejected. Similarly, should half of the sum of any neighboring components' FWHM be less than the separation between their velocity centroids (i.e., the FWHMs are not allowed to overlap), the model is rejected. This restriction is imposed to avoid line blending.
3.
A rejected model is replaced with a simpler one: the three-component model is replaced with the two-component one, and a rejected two-component model is replaced with a one-component model.
4.
The replaced models are required to meet the same heuristical criteria as in step 2.
5.
If none of the fitted models are valid, we consider the spectrum to be a nondetection.
6.
Steps 1–5 are repeated for all the spectra in the combined spectral cube, until a collection of ammonia fitting parameters is assembled in a position–position–velocity (PPV) space.

In the above procedure, the second step requires imposing constraints for both (1, 1) and (2, 2) lines. This is done to assure that the physical properties that depend on the population ratio of the two transitions (kinetic temperature and ammonia column density) are reliably constrained. Throughout this paper, this is the method we refer to as the strict censoring criterion for arriving at the PPV structure in the IRDC. However, we note that in principle, should only the kinematic information be needed from our ammonia data, either one of the inversion line detections would constrain it in our simultaneous line fitting. Therefore, we consider an additional method, which we call relaxed censoring, which naturally results in a PPV structure with a larger average number of line-of-sight components. Despite the differences in these two approaches, the average kinematic properties we will derive later will not be significantly changed if a relaxed censoring method is used (Appendix C provides quantitative comparison between the two approaches); therefore, we use the strict censoring approach because it allows us to discuss the temperature structure of the velocity components in addition to their kinematics.

3. Results

3.1. Velocity Components

Figure 2 shows the PPV space with the best-fit line-of-sight components, which demonstrates the complex kinematics of the G035.39 cloud. As different parts of the cloud may represent unique, noninteracting physical structures, it is necessary to identify the components that are linked in both position and velocity. The inability to separate such unique structures would hinder quantitative analyses on the gas motions within the IRDC.

Since nonhierarchical algorithms have problems with extended emission (e.g., Pineda et al. 2009), we use a hierarchical clustering of discrete PPV data from the G035.39 data set. For this task, we make use of the Agglomerative Clustering for ORganising Nested Structures (ACORNS), a clustering algorithm to be fully described in J.D. Henshaw et al. (2018, in preparation; see Barnes et al. 2018). In simple terms, ACORNS links PPV points together into clusters based on a hierarchical agglomerative clustering. The algorithm starts with the most significant data point, consequently agglomerating surrounding points in the order of descending significance¹² until a hierarchy is established. After the first cluster has been assigned in such a fashion, the most significant point in the unassigned data set is taken, and the linking is repeated. The procedure runs until no data points are left unassigned, and a dendrogram-like cluster hierarchy, assigning clusters to branches and trees (in a manner similar to Rosolowsky et al. 2008b), is established.

The parameters that defined the linking criteria of the ACORNS algorithm were the following. The clusters were restricted to be at least 9 pixels in size to avoid spurious components, and the link between neighboring data points was set to be smaller than or equal to the velocity resolution of the data (0.2 ). Additionally, no links with line widths of twice the velocity resolution (0.4 ) were allowed to be established to prevent spatially close yet distinct components from being merged. In our ACORNS run, the PPV points below S/N significance of 3 were set to be excluded.

We obtain two sets of results from the ACORNS clustering for the two masking strategies introduced in Section 2.2. The strict censoring results in 13 hierarchical clusters, 7 of which cover an area in excess of five synthesized VLA beams (i.e., have physical properties constrained in more than 116 pixels). The overall view of the structures identified is shown in Figure 3, while the separated maps of the seven largest components will be presented in the following sections. The relaxed masking results in 24 components, of which 10 are spanning a significantly large area (a comparison with the adopted method is shown in Appendix C).

**Figure 3.** PPV diagram of the fitted velocity components within IRDC G035.39 along the R.A. projection. The coordinates are given in arcsecond offset relative to the α(J2000) = 18^h57^m08^s, δ(J2000) = +2°10'30'' coordinate. All the data are plotted in black, similarly to Figure 2, and individual velocity components are marked in different colors. The data not found to be associated with any clusters are plotted in gray. The figure shows the clustering obtained with the strict masking criteria (introduced in Section 2.2); the equivalent figure for relaxed censoring can be found in Appendix C. In addition to the R.A. projection, a projection along decl. is shown in the inset axis for the F1 filament.
Download figure:
Standard image High-resolution image

3.1.1. Continuum Sources within the Velocity Components

Throughout this study, we will relate the continuum sources discovered toward the IRDC, in particular the mass surface density cores of Butler & Tan (2012) and the 70 μm sources of Nguyen Luong et al. (2011), to the continuous ammonia-emitting regions identified as velocity-coherent structures in PPV space. In order to establish whether a given continuum source is related to the ACORNS-identified components, we use the following heuristics. First, we require at least 24 pixels with significant ammonia spectra detections (corresponding to an effective area larger than that of one VLA beam) to be enclosed in a three-beam diameter around the coordinate of interest. Second, to assert the direct proximity of the continuum source to the velocity component under consideration, we require at least one significant ammonia spectra detection within one beam radius around the position of interest. These heuristics are computed on every continuum source of interest and for every velocity component detected. Sources that match the criteria above will be referred to as sources related to the ACORNS component throughout this study.

Table 1 gives a summary of the velocity components we identify in this section and lists their corresponding related sources, along with the average properties of the components. The component names follow their respective colors in Figure 3, and only those spanning in excess of five beam sizes are listed. In case a source appears to be related to more than one component, it is listed several times in the table.

Table 1. The Velocity Components Identified with ACORNS, Listed by Decreasing Spatial Extent

Component ID^a	Size^b	^b	σ^b	T_kin^b	Related Sources
	(pixels)	()	()	(K)	(Butler & Tan 2012)	(Nguyen Luong et al. 2011)
F2 (Blue)	2981	45.62 ± 0.41	0.28 ± 0.14	12.29 ± 2.18	H5, H6	#12, #28, #6, #20, #5
F4 (Red)	1377	45.46 ± 0.30	0.31 ± 0.12	11.69 ± 1.70	H2, H4	#22, #10, #16
F5 (Green)	927	44.96 ± 0.32	0.32 ± 0.14	12.40 ± 2.03	H1	#7
F1 (Orange)	313	42.52 ± 0.17	0.17 ± 0.06	11.22 ± 1.29	H3	#11
F6 (Brown)	198	45.86 ± 0.10	0.33 ± 0.25	11.20 ± 1.76	⋯	⋯
F3 (Purple)	270	46.79 ± 0.18	0.24 ± 0.06	12.12 ± 1.90	⋯	#18
F7 (Olive)	124	45.14 ± 0.09	0.15 ± 0.05	13.47 ± 1.94	H6	#6

Notes.

^aThe color labels directly correspond to the component colors in Figure 3. ^bThe mean followed by the standard deviation of the values belonging to the given component.

Download table as: ASCII Typeset image

3.2. Velocity Structure across the IRDC

Various processes in star-forming clouds, such as turbulence, gravitational acceleration of infalling material, feedback from the embedded sources, or rotation, can create line-of-sight velocity gradients. In order to perform a qualitative analysis of the velocity gradients, we have to interpolate the slope of discrete measurements on the plane of sky. To achieve this, we use a variation of the moving least-squares method (Lancaster & Salkauskas 1981) for constructing continuous surface from a sample of discrete data. The method, often used in computational geometry and computer vision (e.g., Levin 1998; Schaefer et al. 2006), was independently developed in astrophysics by Goodman et al. (1993) to find velocity gradients from line-of-sight velocity fields. The automated local velocity gradient fit, which applies the Goodman et al. method across the cloud to study the internal velocity field, was developed by Caselli et al. (2002a) and has seen considerable usage since (e.g., Barranco & Goodman 1998; Caselli et al. 2002b; Tafalla et al. 2004; Kirk et al. 2013; Henshaw et al. 2014, 2016).

The local line-of-sight velocity gradient is fit to a 2D scalar field given by the nearby PPV points, approximated by a first-degree bivariate polynomial f(α, δ) = v + aΔα + bΔδ, where the R.A. and decl. offsets Δα and Δδ are confined to a circle of radius R around a given set of sky coordinates and (v, a, b) are the free parameters of the least-squares fit (v is the centroid velocity, and is the velocity gradient magnitude). For each velocity component identified and for every position we construct a set of valid pixels as follows. First, only the pixels with a well-defined value (i.e., the relevant spectral component having a peak S/N > 3) are included in the velocity gradient analysis, and the rest are masked. Second, the moving least-squares method operates only on the values within a three-beam diameter (i.e., R = 1 farcs 5 × 5 farcs 44) and only if more than 46 centroid velocity measurements (i.e., filling an effective area larger than two VLA beams) are defined within it. Finally, only spatially continuous regions are considered, and regions within radius R that are not connected with position with valid pixels are masked as well. For every valid point with a fitted value, an optimal set of parameters is obtained by a weighted least-squares fit to a collection of nearby pixels:

where is a vector on the plane of sky, and the weights of the fit at a given position, were set to be the fitting uncertainties on the .

Following Goodman et al. (1993), we derive the velocity gradient magnitude and orientation θ for all the velocity components. Additionally, we propagate each pair of best-fit (a, b) values into uncertainties on and θ. The dominant majority (96%) of the values derived are statistically significant (S/N > 3), with a mean value of uncertainties on of 0.06 . Figure 4 shows the for the four largest velocity components (F2, F4, F5, and F1), overlaid with the velocity gradients derived for them.

**Figure 4.** Map of derived V_lsr values for the biggest velocity components (left panel: F2; right panels, top to bottom: F1, F5, F4), overlaid with the velocity gradient arrows. The directions of the arrows point to the steepest velocity field change in the redshifted direction, while the arrow lengths represent the relative vector magnitudes of the gradient. The open and filled red stars denote the positions of 70 μm *Herschel* sources from Nguyen Luong et al. (2011) below and above 20 M_⊙, respectively, and the red diamonds indicate the position of cores from Butler & Tan (2009, 2012). The overlaid contours indicate the highest extinction contours from Kainulainen & Tan (2013), ranging from A_V = 40 to 120 mag, progressing in steps of 20 mag. The relative position on the plane of sky of the components shown is displayed in Appendix E. The data used to create this figure are available.
Download figure:
Standard image High-resolution image

The overall velocity gradient magnitude distribution peaks at 0.35 (Figure 5; values with S/N < 3 are not shown) and shows an asymmetrical tail toward larger values. To study the potential effect of continuum cores dynamically impacting the magnitude of large-scale dense gas motions in the IRDC, we have selected a subsample of the velocity gradient magnitudes that we consider to be spatially relevant to the Butler & Tan (2009, 2012) cores. While the distribution of the core-adjacent values appears to be narrower than the overall distribution, possibly due to smaller dynamic ranges of density probed, the core-adjacent values are representative for the whole cloud with no statistically significant enhancements.

**Figure 5.** Kernel density estimate (KDE) of the velocity gradient magnitudes for all the velocity components derived. In addition to the total distribution of shown in blue, values within three-beam-diameter separation from Butler & Tan (2009, 2012) cores (following the selection criteria outlined in Section 3.1.1) are shown in orange. The KDE bandwidth selection was performed following Scott's rule (Scott 1992). A histogram of statistically significant values (S/N > 3) of velocity gradient magnitudes in low-mass dense cores is shown for comparison (adapted from Goodman et al. 1993).
Download figure:
Standard image High-resolution image

When compared to studies of nearby low-mass dense cores probing similar spatial scales, the gradient magnitude values we find are consistently smaller than those of Caselli et al. (2002a), who have performed a similar velocity gradient analysis on a large sample of low-mass dense cores with a dense gas tracer N₂H⁺, finding typical velocity gradients of 0.8–2.3 . The difference can be attributed to the higher spatial resolution used in the study (corresponding to about 0.04 pc, calculated by taking a median distance to the core sample). A closer comparison can be made with the Goodman et al. (1993) results, which are obtained both with the same molecular transition, NH₃ (1, 1), and with similar physical scales probed (0.07 pc from the median distance to the Goodman et al. core sample), yielding the velocity gradient magnitudes of 0.3–2.5 , which are broadly consistent with our findings in G035.39 (Figure 5).

3.3. Temperature Maps

Following the line-fitting and component separation procedures outlined in Sections 2.2 and 3.1, respectively, we arrive at the highest angular resolution temperature maps for G035.39. Furthermore, the sensitivity of our combined observations allows us to resolve the temperatures for the line-of-sight components in the IRDC. Figure 6 shows an overall distribution of the kinetic temperature values, while the spatial distribution of the ammonia-derived gas kinetic temperature values is presented in Figure 7 for the four largest IRDC components (F2, F4, F5, and F1). The mean of the overall temperatures for all the spectra analyzed in G035.39 is 11.8 K, with a standard deviation of 2.9 K, and the corresponding statistics for individual components are listed in Table 1.

In the GBT-only analysis, Sokolov et al. (2017) identified localized heating at the location of massive protostellar sources. This effect, however, was only seen as a collective contribution due to limited angular resolution of their temperature map. Here the increased angular resolution and reliable detection of both ammonia inversion lines enable a discussion on the individual contribution of the 70 μm sources to the temperature enhancements seen in the IRDC. Figure 6 shows a comparison between the T_kin distribution and the temperature values within one VLA beam around individual 70 μm sources. Compared to the overall distribution of the kinetic temperature, which peaks at 11.7 K and has typical values in the 9–15 K range, the T_kin values in the vicinity of the 70 μm sources show slight temperature enhancements. Temperatures consistently lower than the mean (by about 2 K) are found toward the neck-like structure connecting the north and south parts of the IRDC (left subfigure on Figure 7), which is largely devoid of protostellar sources, with the exception of one 70 μm source on the western edge of the filament, previously classified as a low-mass dense core (source #12 in Nguyen Luong et al. 2011). The source shows moderate temperature enhancements around it, peaking at about 15 K relative to the overall mean of the filament of 10.5 K. The excellent correspondence of the position of the 70 μm source and the localized temperature enhancement in an otherwise cold filament shows that the protostellar source is indeed embedded in the filament, ruling out the possibility of the source being a line-of-sight projection on it.

**Figure 6.** The box plots show the systematic temperature enhancements seen within one VLA beam around each of the 70 μm sources (following the selection criteria outlined in Section 3.1.1). The source number follows Nguyen Luong et al. (2011), with the low-mass dense cores labeled alongside the massive (M > 20 M_⊙) dense cores in bold. The overall distribution (KDE; solid line) of the kinetic temperature values for all the velocity components derived and its 25th, 50th, and 75th percentiles (dotted lines) are shown on the right-hand side of the figure. The KDE bandwidth selection was performed following Scott's rule (Scott 1992).
Download figure:
Standard image High-resolution image

The protostellar source in the northeastern part of the IRDC (F2 component; left panel of Figure 7), classified as a low-mass dense core (9 ± 1 M_⊙; Nguyen Luong et al. 2011), is one of the warmest regions traced with our temperature maps. It is one of the most active among the protostellar sources in the IRDC, manifesting 8, 24, and 70 μm point-source emission. The gas kinetic temperature increases by 8 K in its vicinity, with the peak of 20 K coinciding with the location of the 70 μm point-source emission (Figure 7). With the velocity centroid of the ammonia line of 45.1 continuously morphing with the kinematics of the rest of the cloud (Figure 4), the corresponding temperature enhancement is an argument in favor of the source being embedded in the dense filament traced with our ammonia observations. The puzzling nature of the source is emphasized by the ammonia line widths found toward it. As shown in Sokolov et al. (2018), it coincides spatially with an extended region of subsonic turbulence and features some of the narrowest nonthermal velocity dispersions in the whole data set.

The velocity feature at 42.5 , labeled as the F1 filament in our study (expanded into a subfigure in Figure 3), forms a discontinuity in PPV space with respect to the rest of the IRDC kinematics, despite appearing to trace the southernmost end of Filament 1 of Jiménez-Serra et al. (2014). This blueshifted component shows a velocity swing of about 0.5 toward its center, where an extinction core (Butler & Tan 2009, 2012) is located, and a 70 μm source is found next to it (Nguyen Luong et al. 2011). The angular separation between the 70 μm source and the Butler & Tan (2009, 2012) extinction core is 4 farcs 3, which is smaller than both the beam of our ammonia observations and the 70 μm band of the SPIRE instrument. The two sources are thus likely to be related. The 70 μm source appears as a temperature enhancement in the temperature map, which suggests that star formation is ongoing in the F1 component.

We find the H6 core in the north of the cloud to be associated with two ACORNS components, F2 and F7, due to the VLA observations resolving two continuous overlapping line-of-sight regions. We report a 7 K increase in gas kinetic temperature in the F7 component, coinciding with the extinction peak of Butler & Tan (2012) and of similar size to the H6 core (Figure 8). The excellent correspondence between the temperature enhancement in the F7 components and the location of the core allows us to pinpoint the H6 location in the PPV space. The exact physical nature of the temperature enhancement is less certain, and we suggest that heating from embedded protostars, either as direct radiation or in the form of unresolved outflows, might stand behind the heating. Despite the temperature enhancement seen in the F7 component, we do not observe significant broadening of the line at the H6 position, and the line width, in fact, points to subsonic motions as will be seen in the next section.

**Figure 8.** Zoom-in of the H6 core location for two overlapping line-of-sight velocity components, F2 (left panel) and F7 (right panel), showing derived kinetic temperatures. The black circle marks the fitted size of the H6 core in Butler & Tan (2012).
Download figure:
Standard image High-resolution image

In G035.39, the point-like 70 μm emission, signifying protostellar activity, is found in all but one among the Butler & Tan (2012) cores. Located at the peak of the extinction map in the south of the F5 filament, the H1 core appears as a cold feature in the VLA temperature map (Figure 7), with a mean kinetic temperature of 10.7 K within the VLA beam, and has a complex kinematic morphology (Figure 4). Recent observations by Liu et al. (2018) find evidence for a pinched magnetic field morphology toward the core, indicating that the line-of-sight velocity gradients on Figure 4 could be an argument for accretion flows along the filament. The mean kinetic temperature of the H1 core agrees well with that of the low-mass dense cores (e.g., 11 K for a sample of 193 cores in Rosolowsky et al. 2008a).

4. Discussion

4.1. Subsonic Line Widths in G035.39

Turbulence, dominating molecular clouds on large scales (>1 pc), is found to dissipate toward smaller scales (Larson 1981), allowing for bound prestellar core formation and, ultimately, star formation. Turbulence dissipation has been detected via high-J rotational transitions of CO in quiescent clouds and IRDCs (Pon et al. 2014, 2015, 2016). Subsonic turbulence is observed throughout low-mass star-forming regions and is a natural consequence of turbulent energy dissipating toward smaller scales (Kolmogorov 1941). Observationally, it is seen both as a transition to coherence on small scales (within 0.1 pc) as coherent cores (e.g., Pineda et al. 2010, 2011) and as large, filament-scale subsonic regions in Taurus, Musca, and Serpens South (Tafalla & Hacar 2015; Friesen et al. 2016; Hacar et al. 2016).

Massive star formation, on the other hand, is thought to be governed by nonthermal motions. However, in IRDC G035.39 Sokolov et al. (2018) find that substantial regions of the cloud are predominantly subsonic (i.e., with sonic Mach number ). While the implications for the turbulence dissipation in other IRDCs are discussed in Sokolov et al. (2018), we cite the fraction of subsonic spectra identified in the ACORNS velocity components below. Table 2 lists the average Mach numbers, the spread of their distribution (one standard deviation), and the overall fraction of subsonic components relative to the total in each of the velocity components spanning in excess of five VLA beams. In addition, details of the Mach number derivation, as well as the propagation of uncertainty on matching that of Sokolov et al. (2018), are included as an appendix item in Appendix D.

Table 2. Subsonic Motions within the G035.39 Velocity Components

Component ID^a	^b	Fraction (%)
F2 (Blue)	1.24 ± 0.70	44
F4 (Red)	1.41 ± 0.62	32
F5 (Green)	1.42 ± 0.67	30
F1 (Orange)	0.74 ± 0.30	81
F6 (Brown)	1.54 ± 1.41	55
F3 (Purple)	1.00 ± 0.32	48
F7 (Olive)	0.54 ± 0.25	98

Notes.

^aThe color labels directly correspond to the component colors in Figure 3. ^bThe mean followed by the standard deviation of the values belonging to the given component.

Download table as: ASCII Typeset image

We report no significant trends in the Mach numbers against temperature or density (either column density or the extinction-based mass surface density). However, we note that in two velocity components that are associated with star formation based on both the presence of 70 μm emission and kinetic temperature increase (F1 and F7; Figures 7, 8) we observe a larger fraction of subsonic motions than in any other component (Table 2). We suggest that future VLA observations are needed to fully resolve the boundary over which the transition to coherence in the IRDC occurs, its relation to star formation, and the degree of similarity to the transition to coherence in low-mass star-forming clouds (e.g., Goodman et al. 1998; Pineda et al. 2010), which could further probe smaller-scale fragmentation (Pineda et al. 2015).

4.2. Comparison with N₂H⁺

The region of G035.39 mapped by Henshaw et al. (2014) with their Plateau de Bure Interferometer (PdBI) observations of N₂H⁺ (1−0) line transition is fully covered in our VLA mosaic. As the spectral setups of both data sets were consistent (final channel separation of 0.2 and 0.14 for VLA and PdBI data, respectively, with no extra smoothing applied), we are well positioned to perform a quantitative comparison of the two tracers of dense gas kinematics. Assuming that the NH₃ and N₂H⁺ molecules are tracing similar gas, the turbulent (or nonthermal) components of their line widths should be equal. Figure 9 shows a direct comparison between the two nonthermal line widths. As the N₂H⁺ molecule is heavier than NH₃, we have subtracted the thermal component from both transitions in making the figure for a more direct comparison, taking T_kin from the mean ammonia kinetic temperature for both molecules to assure the uniform comparison between the two data sets. In other aspects, the derivation of the nonthermal N₂H⁺ line width follows our calculations for the ammonia line widths (Appendix D), with substitutions made for molecular mass where appropriate. Only the data from the commonly mapped region are shown in Figure 9, regardless of the number of line-of-sight components found in either data set.

**Figure 9.** Comparison between the FWHM line widths derived in Henshaw et al. (2014) and those derived in our work. As the VLA observation mosaic covers a larger area, only the values overlapping with Henshaw et al. (2014) coverage are included.
Download figure:
Standard image High-resolution image

Surprisingly, we find the nonthermal line widths derived from our VLA observations to be about 20% narrower than those of the PdBI N₂H⁺ (1−0) from Henshaw et al. (2014). This is unexpected, as Henshaw et al. (2014) data resolution is comparable to our common VLA restoring beam size of 5 farcs 44. We also note that both data sets have also had the missing flux filled in from complementary single-dish observations. In addition, the critical densities of the NH₃ (1, 1) and (2, 2) transitions are an order of magnitude lower than that of N₂H⁺ (1−0) (Shirley 2015), so if the higher-density gas is less turbulent, a broader line width should be seen in the ammonia spectra, contrary to our results. However, ammonia can selectively trace high-density regions, as its numerous hyperfine components allow us to trace gas also at densities significantly larger than the critical density and its abundances increase toward the core centers (Tafalla et al. 2002; Crapsi et al. 2007; Caselli et al. 2017); thus, NH₃ is expected to follow material traced by N₂H⁺ (1−0) as found in low-mass star-forming regions (Benson et al. 1998). Finally, the presence of magnetic fields is expected to affect the line widths of the ion species, trapping the ion molecules and narrowing their line profiles (Houde et al. 2000a), a trend opposite to what we find here. However, the phenomenon requires large turbulent flows to occur in directions misaligned with respect to the magnetic field lines, and it was not found to come into play in dark clouds primarily supported by thermal motions (Houde et al. 2000b).

Two possible explanations of the line width discrepancy above are differences in hyperfine structure modeling and varying chemical differentiation. In a detailed study of two starless cores, Tafalla et al. (2004) found a similar discrepancy in the line widths of NH₃ and N₂H⁺ and attributed the differences in nonthermal line widths to non-LTE effects of the N₂H⁺ lines (see Caselli et al. 1995). Because the PdBI data from Henshaw et al. (2014) were fit with an isolated hyperfine component only, modeling it as a simple Gaussian component with τ ≪ 1, the relative broadening of the N₂H⁺ line with respect to ammonia could be an optical depth effect. However, Henshaw et al. (2014) discuss this scenario and make several arguments in favor of the isolated hyperfine transition being optically thin. Another scenario is a potential chemical differentiation between N₂H⁺ and NH₃ across the IRDC, preferentially tracing regions of lower turbulence in ammonia emission. Large differences in abundance profiles of different molecules are commonly found toward low-mass starless cores (e.g., Tafalla et al. 2002; Spezzano et al. 2017) and IRDCs (Feng et al. 2016). For the particular case of NH₃ and N₂H⁺, several observational works have proposed that their abundance ratio could be an evolutionary indicator in dense cold cores (Palau et al. 2007; Fontani et al. 2012a, 2012b). Chemical differentiation could then play a dominant role in setting the line width differences, outweighing the other arguments listed above that would favor a broader ammonia line width. This, however, would require a significant difference in the abundance profiles for the two molecules, a difference that can be quantified by comparing the spatial distribution in the integrated intensities of NH₃ and N₂H⁺. However, both tracers appear to follow the dense gas closely, showing strong linear correlations with mass surface density from Kainulainen & Tan (2013): for NH₃ (1, 1) (Pearson's r = 0.73), NH₃ (2, 2) (r = 0.71), and N₂H⁺ (1−0) (r ∼ 0.7; Henshaw et al. 2014). Therefore, the origin for the discrepancy remains unclear, and it is likely that higher angular resolution observations are needed to understand it.

4.3. Complex Gas Motions in G035.39

What physics processes drive the observed velocity gradients in G035.39? The multicomponent velocity field we derive in this study has a complex structure and is likely to be driven by a combination of turbulence, cloud rotation, global gas motions, and/or gravitational influence of material condensations. To test whether the gravitational influence is the driving factor for the velocity gradients, we perform a simple analysis, taking the mass surface density of Kainulainen & Tan (2013) as a proxy for gravitational potential. In a simple picture of gas infalling into a gravitational well, one would expect an acceleration of the infalling gas toward the mass condensations forming the well, and therefore the observed line-of-sight velocity gradient magnitudes would be enhanced toward such clumps. However, we find no correlation between the local velocity gradient magnitude and the mass surface density from Kainulainen & Tan (2013); furthermore, we find no significant enhancement in the velocity gradient magnitudes toward the locations of mass surface density peaks (Figure 5), indicating that the velocity field we observe has a more complex origin. In light of this, the discussion of the nature of the velocity gradients is limited to regions of IRDC with continuum sources where we see ordered velocity gradients.

In the vicinity of two massive dense cores (MDC 7 and MDC 22) from Nguyen Luong et al. (2011), a smooth velocity gradient pattern is seen across the cores. Observationally, a pattern like this can be modeled as a solid-body rotation (Goodman et al. 1993). Assuming that the linear velocity gradient represents such a rotation, we can take its magnitude as the angular velocity of the rotating core. Both sources show similar velocity field patterns, with typical velocity gradient magnitudes of ∼1 within three VLA beam diameters: the mean value of is 1.05 for MDC 22 (1.01 for MDC 7), while the 1σ spread in the values is 0.15 (0.14 for MDC 7). The uncertainties on the velocity gradient magnitudes are typical of the rest of the IRDC, with mean errors of 0.03 and 0.05 within three VLA beam diameters of MDC 22 and MDC 7, respectively. Taking the deconvolved FWHM sizes of the 70 μm sources from Nguyen Luong et al. (2011) for both sources (0.12 and 0.13 pc for MDC 7 and MDC 22, respectively) as values for R, we arrive at the specific angular momenta of J/M ∼ 2 × 10²¹ cm² s⁻¹, where we consider the individual differences between the two sources to be lost in the uncertainties of the calculation. The latter consideration stems from a large number of assumptions taken prior to calculation of the specific angular momentum, such as implicitly assuming a spherical radius for the cores. Furthermore, we have assumed the dimensionless parameter p = 0.4 for uniform density and solid-body rotation (Goodman et al. 1993) in our J/M calculation, although significant (factor of 2–3) deviations are expected for turbulent, centrally concentrated cores (Burkert & Bodenheimer 2000). Finally, our implicit interpretation of the physical nature of the velocity gradients can be biased as well. The observed ammonia kinematics trace the entire IRDC, not just its cores, and the velocity gradient vector field we recover with the moving least-squares method is sensitive not only to ∼0.1 pc core motions but also to gas motions on larger scales. The seemingly ordered velocity field around MDC 7 and MDC 22 does not necessarily represent rotation of the cores. It could be a manifestation of the larger, filamentary rotation, an effect of ordered gas flows on larger scales, or a turbulent component at a larger scale. As low-mass dense cores have been found to have gradient directions that differ from the gradients at larger, cloud scales (Barranco & Goodman 1998; Punanova et al. 2018), the motions we see could be dominated by these large scales. Whether this is indeed the case would have to be verified with higher angular resolution observations of higher gas density tracers.

Despite the uncertainties discussed above, it is interesting to see how the specific angular momentum compares to the low-mass star formation. In low-mass star-forming cores, rotational motions have been shown to be insufficiently strong to provide significant stabilizing support against gravity (Caselli et al. 2002a). Compared to their lower-mass counterparts, cores in massive IRDCs are denser and are born in a relatively more pressurized and turbulent environment. If the origin of angular momentum of the star-forming material is related to the degree of nonthermal motions, the specific angular momentum inherited by the IRDC cores could be boosted (see Tatematsu et al. 2016, who find systematically larger J/M values in Orion A cores). Despite this line of reasoning, in IRDC G035.39 we arrive at values of J/M that are fairly typical of the low-mass dense cores (e.g., Crapsi et al. 2007) and align well with the specific angular momentum–size relation in Goodman et al. (1993). The agreement with the low-mass core J/M is consistent with our previous results, where we have shown that G035.39 has a smaller degree of nonthermal motions (Sokolov et al. 2018) than other IRDCs, and implies that, given the larger mass of the IRDC cores, the rotational support is even less dynamically important than for the low-mass dense cores.

5. Conclusions

Our main results can be summarized as follows:

1.
In the first interferometric mapping of the entire IRDC G035.39 cloud, we have mapped the NH₃ (1, 1) and (2, 2) inversion transitions, combining the VLA and GBT observations.
2.
Our data reveal the convoluted kinematics in the cloud, bridging the cloud and core scales. We identify seven principal structures in PPV space and find velocity gradients within them that are comparable to the dynamics of the low-mass star-forming cores imaged with similar resolution.
3.
Contrary to expectations from theoretical observational perspectives, we find a systematic line width difference between NH₃ and N₂H⁺, with ammonia line width being about 30% smaller, a discrepancy also found toward low-mass dense cores.
4.
We construct a gas temperature map for the IRDC, with the highest angular resolution available to date for this source and, wherever deemed significant, for separate line-of-sight velocity components. We report significant heating from the embedded protostellar sources that has not been seen individually from previous studies of gas and dust temperature in G035.39. The temperature of the seemingly starless core in the IRDC is typical of low-mass dense cores.
5.
We have identified a number of sources, of both protostellar and starless nature, that exhibit complex gas motions around them. A few of these sources show an ordered velocity field in their vicinity and appear to influence the dynamics of the gas around them. For the selected sources, we calculate the corresponding specific angular momentum, assuming solid-body rotation of a uniform-density core, to be about 2 × 10²¹ cm² s⁻¹, similar to that of the low-mass dense cores.

V.S., J.E.P., and P.C. acknowledge the support from the European Research Council (ERC; project PALs 320620). K.W. acknowledges support by the National Key Research and Development Program of China (2017YFA0402702), the National Science Foundation of China (11721303), and the German Research Foundation (WA3628-/1). A.T.B. would like to acknowledge the funding from the European Union's Horizon 2020 research and innovation program (grant agreement no. 726384). J.C.T. acknowledges NASA grant 14-ADAP14-0135. This publication made use of ACORNS, a Python package used to link PPV data.

Facilities: VLA - Very Large Array, GBT. -

Software: scipy (Jones et al. 2001), astropy (Astropy Collaboration et al. 2013, 2018), pyspeckit (Ginsburg & Mirocha 2011), APLpy (Robitaille & Bressert 2012), CASA (McMullin et al. 2007), ACORNS (J.D. Henshaw et al. 2018, in preparation).

Appendix A: VLA and GBT Data Combination Strategies

This section expands on the subtleties of the deconvolution process we use to arrive at the final combined VLA+GBT spectral cubes. The CLEAN algorithm (Högbom 1974) deconvolves images by representing them as a number of point sources, iteratively subtracting fractions of a dirty beam from the peak intensity regions, until the stopping threshold has been reached. While the method is best suited for collections of bright point sources, imaging fainter and more extended emission is not a trivial task. Using CLEAN to recover extended emission with surface brightness comparable to that of noise artifacts would often lead only to an amplification of the noise level in the resultant image; hence, a careful approach to the problem is needed in such cases. A common strategy to recover dim extended emission involves restricting the CLEAN algorithm by a spatial mask where a significant level of emission is known to exist, thus preventing iterations of CLEAN to run on noise features. The overall logic behind the procedure we use to image ammonia data in G035.39 is similar but occurs in a few steps—we rely on the knowledge of the small-scale emission first to set up the initial CLEAN mask, subsequently gradually extending emission-based CLEAN masks into the deconvolution process as it progresses. This approach assures the certainty of the extended emission cleaning and will be shown to be more reliable than a simpler version of it. A summarized description of this procedure was presented in Sokolov et al. (2018), as well as in Section 2.

As a first step in the imaging process, we deconvolve the VLA-only visibilities on their own, reliably targeting the spatial distribution compact emission. To do this, we use the multiscale generalization of the CLEAN algorithm (Cornwell 2008) implemented within CASA (McMullin et al. 2007) as the TCLEAN task, using Briggs weighting and the robust parameter of 0.5. In three consequent runs of TCLEAN, we progressively clean the emission down to multiples of the noise level (5σ, 3σ, and, finally, down to 1σ) for both NH₃ (1, 1) and (2, 2) spectral cubes. For the first run of the algorithm, the CLEAN mask is calculated¹³ for each spectral channel based on a dirty image (with a 3σ noise level threshold). For the subsequent TCLEAN runs, we expand the clean mask by using the output image from the previous run at the 2σ noise level threshold.

Second, we feather the cleaned VLA images with the GBT data, using the FEATHER task within the CASA package. We found the default settings for the dish diameter (∼100 m) and the default low-resolution scale factor of unity to produce reasonable results. Prior to feathering, all the GBT spectra are smoothed with a Gaussian kernel to 0.2 channel separation of the VLA data and interpolated onto the spatial map used in the VLA imaging (1'' pixels in a 240 × 384 grid). Cubic interpolation was used, as the nearest-neighbor interpolation was found to produce artifacts resulting from a large relative size of the GBT pixels (∼10''). While these data are already sufficient for conducting scientific analyses, we found that they can be significantly improved through a procedure outlined below. In order to demonstrate the improvement, the feathered data will be compared with the final imaging results in the text concluding this section.

Finally, armed with the reliable information about the extended emission spatial distribution in the feathered spectral cubes, we perform a final deconvolution of the VLA visibilities. For this final imaging step, we use the rescaled GBT images as a starting model. The multiscale CLEAN proceeds in the same fashion as in the previous runs, with the exception of the initial clean mask, which is now being calculated on the feathered cube from step 2.

Given the inherently uncertain nature of the VLA imaging process, it is anything but trivial to judge the degree of improvement in final image quality. While overall the results are similar to a simpler feathering approach, there are important differences that make us prefer the new method. A comparison of the integrated intensities between the new and old imaging results, overlaid with significantly detected (S/N = 3) contours, is shown in Figure 10. While the overall morphology of the two approaches is similar, the method we use results in fewer spurious noise-driven features. The latter claim is justified by a lower number of isolated S/N > 3 closed contours off the main filament of G035.39, shown for comparison in Figure 10.

**Figure 10.** Integrated intensities of the NH₃ (1, 1) and (2, 2) lines for the two imaging setups we consider: the chosen method is shown in the left panel, while a simpler feathered CLEAN run is shown in the right panel. Solid black contours mark S/N = 3 detection in the integrated intensity. Overlaid in white is the mid-IR extinction contour at A_V = 25 mag, arbitrarily chosen to represent the cloud border (Kainulainen & Tan 2013).
Download figure:
Standard image High-resolution image

To quantify the degree of improvement, we use the fact the (1, 1) and (2, 2) lines of NH₃ were imaged independently in separate sessions. Therefore, if we assume that the integrated emission of the two lines originates from the same physical region, tracing the dense cold gas making up the IRDC, we can use the correlation between the two as a proxy for testing how well the deconvolution has reconstructed the ammonia emission. Note that while the temperature differences across G035.39 will reduce the overall correlation, such a bias is expected to come into play for both the imaging setups. Figure 11 shows a comparison of how the correlation in integrated intensities varies for the different S/N cuts (i.e., on the y-axis we check for correlation of all pixels with S/N > x). While significant (r > 0.5) linear correlation is present between (1, 1) and (2, 2) integrated intensities for all the S/N ranges, the correlation is consistently higher for the imaging setup we use in this work.

**Figure 11.** Pearson's r correlation coefficient between the integrated intensities of the (1, 1) and (2, 2) lines. The deconvolution setup we use in this work (in orange) results in a higher correlation between independent data cubes than a simpler approach (in blue), indicating that our approach is more reliable.
Download figure:
Standard image High-resolution image

Appendix B: Channel Maps of the NH₃ (2, 2) Line

Channel maps for the NH₃ (2, 2) line, presented in the same manner as the NH₃ (1, 1) channel maps in the main text (Figure 1), are presented in Figure 12. Overall, both the integrated intensity structure and the S/N levels agree with those of the NH₃ (1, 1) channel maps.

Appendix C: Relaxed Censoring PPV Structures

The strict censoring approach (Section 2.2) we take when deriving the PPV structure of the cloud could selectively filter out weaker secondary line-of-sight components. If a relaxed censoring strategy is used instead, the PPV structure would naturally result in more line-of-sight components. Because we require both NH₃ lines to be significantly detected to derive the kinetic temperature and constrain the column density, the relaxed censoring is not considered as a reliable approach throughout this paper. Nonetheless, if only the kinematic structure of the IRDC is to be discussed, the relaxed censoring could better capture the velocity patterns in its filaments. In this context, if we consider either (1, 1) or (2, 2) line detection to constrain the gas kinematics, how would this change the basic kinematic properties we derive throughout this work?

For the same ACORNS clustering criteria used, the hierarchical procedure outlined in Section 3.1 results in more distinct velocity structures: 24 components, of which 10 exceed the collective area of five synthesized VLA beams. Figure 13 shows the ACORNS components identified in a manner similar to Figure 3. The velocity gradient magnitude distribution, if derived in the same manner as in Section 3.2, has a mean value of 0.72 , close to 0.75 obtained in Section 3.2. The spread in the distributions of velocity gradient magnitudes is similar as well: with standard deviations of 0.42 and 0.45 for relaxed and strict censoring strategies, respectively.

**Figure 13.** ]PPV diagram of the fitted velocity components within IRDC G035.39 along the R.A. projection. The coordinates are given in arcsecond offset relative to the α(J2000) = 18^h57^m08^s, δ(J2000) = +2°10'30'' coordinate. All the data are plotted in black, similarly to Figure 2, and individual velocity components are marked in different colors. The data not found to be associated with any clusters are plotted in gray. The figure shows the clustering obtained with the relaxed masking criteria (introduced in Section 2.2). In addition to the R.A. projection, a projection along decl. is shown in the inset axis for the F1 filament.
Download figure:
Standard image High-resolution image

Given the similarity of the derived kinematic properties, we argue that the subsequent discussion on IRDC kinematics in Section 4.3, namely, on the specific angular momenta, is unaffected by the method chosen to arrive at the full PPV structure of G035.39.

Appendix D: Nonthermal Line Widths and Propagation of Uncertainties in Mach Numbers

To derive the degree of nonthermal gas motions and gas turbulence in G035.39 from the ammonia data, we use the following relations. First, accounting for the channel width broadening effects, we subtract the channel width from the observed FWHM maximum of the line:

where is the fitted line width from the observed line profile and is taken to be the channel separation of our VLA setup at 0.2 . To find the degree of nonthermal gas motions, we then remove the component caused by the thermal broadening of the line, obtaining the nonthermal velocity dispersion from the line width (Myers et al. 1991):

where is the molecular weight of the ammonia molecule.

A useful way to quantify the degree of the nonthermal motions in the gas is to compare those against the local sonic gas speed , where , the average mass of the particle in the gas medium, is taken to be 2.33 amu. Following this, a sonic Mach number is computed as the ratio of the two quantities: .

In order to enable a statistically sound discussion on the Mach numbers, one needs to first estimate the uncertainty on their values. In our analyses, we assume the uncertainty to come from two factors: the temperature of the gas, used to calculate the thermal contribution to line width and the sonic sound speed in the medium, and the error on the line width itself. Both errors are estimated within the pyspeckit package, which provides 1σ uncertainties on T_kin and σ. Assuming these uncertainties to be uncorrelated, we propagate the errors as follows:

where the notations and are the uncertainties on the gas kinetic temperatures and observed velocity dispersion, respectively.

Substituting appropriate terms into Equation (1), we arrive at the following expression:

where the terms and are introduced for clarity.

The expression above is used to quantify the uncertainties on the values throughout this work (where applicable), as well as in Sokolov et al. (2018).

Appendix E: Relative Orientation of the Velocity Components

The relative arrangement of the identified velocity components in PPV space can be misleading when projected along one of the axes. As throughout this work the velocity components we identify are normally visualized in the decl.–velocity plane, it might not be completely clear as to how exactly the components, plotted individually, appear relatively to each other on the plane of sky. To address this potential issue, Figure 14 shows the relative sky position of the velocity components plotted as subpanels in Figures 4 and 7.

**Figure 14.** Integrated intensity of the NH₃ (1, 1) line, outlining the overall structure of G035.39, overplotted with rectangles that each refer to the actual on-sky position and extent of every velocity component displayed in Figures 4 and 7. The different colors are identical to the ones used in Figure 3.
Download figure:
Standard image High-resolution image

Multicomponent Kinematics in a Massive Filamentary Infrared Dark Cloud

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction