Multicomponent kinematics in a massive filamentary IRDC

TL;DR

This study investigates the complex velocity structures and motions within dense filaments of a massive infrared dark cloud, revealing similarities with low-mass star formation processes and providing detailed kinematic measurements.

Contribution

It introduces a detailed analysis of multicomponent kinematics in a high-mass IRDC using combined VLA and GBT data, highlighting the presence of coherent structures and comparing their motions to low-mass cores.

Findings

Identified seven extended velocity-coherent components in the IRDC.

Found velocity gradients consistent with core rotation and similar to low-mass cores.

Discovered ammonia line widths are narrower than N2H+ lines, contrary to expectations.

Abstract

To probe the initial conditions for high-mass star and cluster formation, we investigate the properties of dense filaments within the infrared dark cloud G035.39-00.33 (IRDC G035.39) in a combined Very Large Array (VLA) and the Green Bank Telescope (GBT) mosaic tracing the NH3 (1,1) and (2,2) emission down to 0.08 pc scales. Using agglomerative hierarchical clustering on multiple line-of-sight velocity component fitting results, we identify seven extended velocity-coherent components in our data, likely representing spatially coherent physical structures, some exhibiting complex gas motions. The velocity gradient magnitude distribution peaks at its mode of 0.35 km/s/pc and has a long tail extending into higher values of 1.5 - 2 km/s/pc, and is generally consistent with those found toward the same cloud in other molecular tracers and with the values found towards nearby low-mass dense cloud cores at the same scales. Contrary to observational and theoretical expectations, we find the non-thermal ammonia line widths to be systematically narrower (by about 20%) than those of N2H+ (1-0) line transition observed with similar resolution. If the observed ordered velocity gradients represent the core envelope solid-body rotation, we estimate the specific angular momentum to be about 2 x 10^21 cm^2/s, similar to the low-mass star-forming cores. Together with the previous finding of subsonic motions in G035.39, our results demonstrate high levels of similarity between kinematics of a high-mass star-forming IRDC and the low-mass star formation regime.

Tables

Table 1. Table 1: The velocity components identified with ACORNS, listed by decreasing spatial extent.

Component IDaaThe color labels directly correspond to the component colors on Fig. 3.	SizebbThe mean followed by the standard deviation of the values belonging to the given component.	$V_{\mathrm{lsr}}$ bbThe mean followed by the standard deviation of the values belonging to the given component.	$\sigma$ bbThe mean followed by the standard deviation of the values belonging to the given component.	$T_{\mathrm{kin}}$ bbThe mean followed by the standard deviation of the values belonging to the given component.	Related sources
	(pixels)	($\mathrm{km}{s}^{-1}$)	($\mathrm{km}{s}^{-1}$)	(K)	(Butler & Tan, 2012)	(Nguyen Luong et al., 2011)
F2 (Blue)	2981	$45.62\pm 0.41$	$0.28\pm 0.14$	$12.29\pm 2.18$	H5, H6	#12, #28, #6, #20, #5
F4 (Red)	1377	$45.46\pm 0.30$	$0.31\pm 0.12$	$11.69\pm 1.70$	H2, H4	#22, #10, #16
F5 (Green)	927	$44.96\pm 0.32$	$0.32\pm 0.14$	$12.40\pm 2.03$	H1	#7
F1 (Orange)	313	$42.52\pm 0.17$	$0.17\pm 0.06$	$11.22\pm 1.29$	H3	#11
F6 (Brown)	198	$45.86\pm 0.10$	$0.33\pm 0.25$	$11.20\pm 1.76$	-	-
F3 (Purple)	270	$46.79\pm 0.18$	$0.24\pm 0.06$	$12.12\pm 1.90$	-	#18
F7 (Olive)	124	$45.14\pm 0.09$	$0.15\pm 0.05$	$13.47\pm 1.94$	H6	#6

Table 2. Table 2: Subsonic motions within the G035.39 velocity components.

Component IDaaThe color labels directly correspond to the component colors on Fig. 3.	$ℳ$bbThe mean followed by the standard deviation of the values belonging to the given component.	$ℳ\le 1$ fraction (%)
F2 (Blue)	$1.24\pm 0.70$	44
F4 (Red)	$1.41\pm 0.62$	32
F5 (Green)	$1.42\pm 0.67$	30
F1 (Orange)	$0.74\pm 0.30$	81
F6 (Brown)	$1.54\pm 1.41$	55
F3 (Purple)	$1.00\pm 0.32$	48
F7 (Olive)	$0.54\pm 0.25$	98

Equations

$$\mathbf{p}(\mathbf{r})=\underset{a,b}{\mathrm{argmin}}\left\{\sum _{{\mathbf{r}}^{\prime }:\Vert \mathbf{r}-{\mathbf{r}}^{\prime }\Vert \le R}\frac{(v_{\mathrm{lsr}}({\mathbf{r}}^{\prime })-f({\mathbf{r}}^{\prime },a,b){)}^2}{w({\mathbf{r}}^{\prime }{)}^2}\right\},$$

$$\Delta V_{FWHM}=\sqrt{\Delta\mbox{$V_{\mathrm{obs}}$}^{2}-\Delta\mbox{$V_{\mathrm{chan}}$}^{2}},$$

$$\mbox{$\sigma_{\rm nt}$}=\sqrt{\frac{\Delta\mbox{$V_{\mathrm{obs}}$}^{2}}{8\ln{2}}-\frac{\mbox{$k_{\rm B}$}\mbox{$T_{\mathrm{kin}}$}}{\mbox{$m_{\rm NH_{3}}$}}},$$

$$\mbox{$\sigma_{\mathcal{M}}$}=\sqrt{\mbox{$\sigma_{\sigma_{\rm obs}}$}^{2}\big{[}\frac{\partial\mathcal{M}}{\partial{\mbox{$\sigma_{\rm obs}$}}}\big{]}^{2}+\mbox{$\sigma_{T_{\rm kin}}$}^{2}\big{[}\frac{\partial\mathcal{M}}{\partial{\mbox{$T_{\mathrm{kin}}$}}}\big{]}^{2}},$$

$$\mbox{$\sigma_{\mathcal{M}}$}=\frac{\sqrt{\mbox{$m_{\rm gas}$}\mbox{$m_{\rm NH_{3}}$}}}{2}\sqrt{\frac{\mbox{$\sigma_{T_{\rm kin}}$}^{2}(\mbox{$\sigma_{\rm obs}$}^{2}-\mbox{$\sigma_{\rm chan}$}^{2})^{2}+4\mbox{$\sigma_{\sigma_{\rm obs}}$}^{2}\mbox{$\sigma_{\rm obs}$}^{2}\mbox{$T_{\mathrm{kin}}$}^{2}}{\mbox{$k_{\rm B}$}\mbox{$T_{\mathrm{kin}}$}^{3}(\mbox{$m_{\rm NH_{3}}$}(\mbox{$\sigma_{\rm obs}$}^{2}-\mbox{$\sigma_{\rm chan}$}^{2})-\mbox{$k_{\rm B}$}\mbox{$T_{\mathrm{kin}}$})}},$$

Full text

Using agglomerative hierarchical clustering , we identify seven extended velocity-coherent components in our data, likely representing spatially coherent physical structures The velocity gradient magnitude distribution peaks at its mode of 0.35 $\mathrm{km~{}s^{-1}~{}pc^{-1}}$ and has a long tail extending into higher values of $1.5-2$ $\mathrm{km~{}s^{-1}~{}pc^{-1}}$, and is generally consistent with those found toward the same cloud in other molecular tracers and with the values found towards nearby low-mass dense cloud cores at the same scales. Contrary to observational and theoretical expectations, we find the line widths to be systematically narrower (by about 20%) than those of N2H+ ($1-0$) line transition observed with similar resolution. If the observed ordered velocity gradients represent the core envelope solid-body rotation, we estimate the specific angular momentum to be about $2\times 10^{21}~{}\mathrm{cm^{2}~{}s^{-1}}$, similar to the low-mass star-forming cores. Together with the previous finding of subsonic motions in G035.39, our results demonstrate high levels of similarity between kinematics of a high-mass star-forming IRDC and the low-mass star formation regime.

ISM: kinematics and dynamics — ISM: clouds — stars: formation — ISM: individual objects: G035.39-00.33

††journal: ApJ††facilities: VLA, 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 (Henshaw et al. in prep).

\turnoffedit

1 Introduction

The formation of stars from the diffuse ISM is a multi-stage 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 (${\sim}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 (${\sim}1$ pc), and are often poorly constrained by observations. As these conditions directly regulate the subsequent fragmentation and infall of the star-forming material (Tan, 2018; Motte et al., 2018), high-angular resolution observations of the dense IR-dark gas with little evidence of disruptive feedback can provide improved constraints on the theoretical models and numerical simulations aiming to recreate 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, or 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 GMCs (Hernandez & Tan, 2015; Schneider et al., 2015), IRDCs appear as cold ($10-15$ K), filamentary, high-column-density ($10^{23}-10^{25}$ cm*-2*) and high-mass ($10^{3-5}$ $M_{\sun}$, 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 shown that the clouds are further fragmented into high-mass cores with a median mass of 120 $M_{\sun}$. Butler & Tan (2009) studied mid-infrared extinction of 10 IRDCs from the Rathborne et al. (2006) sample, selecting the clouds that are nearby, massive, and show relatively simple diffuse emission. Butler & Tan (2012) identified a number of compact ($R\sim 0.1-0.2$ pc) dense cores embedded in these ten 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 ${\sim}6$ pc down to pc) for one of ten 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 NIR and MIR extinction measurements, Kainulainen & Tan (2013) have estimated the mass of G035.39 to be 16700 $M_{\sun}$. The cloud temperature is relatively cold (${\sim}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 appearence on large, pc-scales, Nguyen Luong et al. (2011) find a population of protostellar cores within G035.39, nine of which are massive enough ($>20$ $M_{\sun}$) 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) have 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 its implications for the high-mass star-forming regime.

In this paper, we present a detailed analysis of the combined Very Large Array (VLA) and Green Bank Telescope (GBT) ammonia observations of the whole G035.39 cloud. The paper is structured as follows. In §2 we summarize the data reduction, our strategy for VLA and GBT dataset combination, and present an overview of its outcome. In §3 we describe the separation of the data into individual components, and present the spectral line fitting results for the largest among them. In §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 made, and the stability of the individual components is being discussed. Finally, a summary of our main results comprise section §5.

2 Data Reduction

The Karl G. Jansky Very Large Array (VLA) observations (project AW776; PI: Ke Wang) in compact D-configuration were carried out on May 8, 2010, mapping the NH3 (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 of 15.625 kHz. The data were calibrated on the quasars J1851+005 (gain), J2253+1608 (bandpass), and 3C48 (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 multi-scale CLEAN algorithm (Cornwell, 2008), following the Briggs weighting with 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.44″. To recover the extended ammonia emission, we fill in the missing flux from the Green Bank Telescope (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 of 0.2 $\mathrm{km~{}s^{-1}}$ . 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 (1,1) and (2,2) images were deconvolved again with tclean task, with the tclean mask . 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 dataset combination is presented. The resulting NH3 (1,1) and (2,2) spectral cubes, gridded into 1″ pixels, have a common restoring beam of 5.44″. The typical rms value of the emission free spectra in the resultant cubes is 14 $\mathrm{mJy~{}beam^{-1}}$ for NH3 (1,1) and 5 $\mathrm{mJy~{}beam^{-1}}$ for NH3 (2,2) inversion lines in a 0.2 $\mathrm{km~{}s^{-1}}$ channel (1.04 K and 0.37 K, respectively).

2.1 Overview of the data

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

The parsec-scale ${\sim}0.2$ $\mathrm{km~{}s^{-1}~{}pc^{-1}}$ quasi-linear velocity gradient, found in the GBT-only data (Sokolov et al., 2017) is present in our combined data 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, red-shifted, velocity component located around 47 $\mathrm{km~{}s^{-1}}$, in agreement with the previous studies of the cloud (labelled 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\mbox{--}43$ $\mathrm{km~{}s^{-1}}$, 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 a dozen thousand spectra with S/N ratio higher than three, the common practice of fitting multiple components by eye would require an excessive amount of human interaction resulting in low level of reproducibility, while methods for global non-linear 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 package111https://github.com/vlas-sokolov/multicube 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 pyspeckit package (Ginsburg & Mirocha, 2011), are split into a multi-dimensional 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 models222Upon close inspection, no spectra containing more than three velocity components were identified.. The three fitting sets are then merged into one following a set of heuristical criteria that determines the number of line-of-sight components in a given pixel of the spectral cube. The heuristical procedure is described below:

First, the three component model is considered to be valid. 2. 2.

If any of the three components are below signal-to-noise of three — for both NH3 (1,1) and (2,2) transitions — the model is marked as rejected. Similarly, should the be less than the separation between their velocity centroids (i.e., the FWHM are not allowed to overlap), the model is rejected. This restriction is imposed to avoid line blending. 3. 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. 4.

The replaced models are required to meet the same heuristical criteria as in step 2. 5. 5.

If none of the fitted models are valid, we consider the spectrum to be a non-detection. 6. 6.

The 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 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 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 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 position-position-velocity (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, non-interacting 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 non-hierarchical 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 dataset. For this task, we make use of the Agglomerative Clustering for ORganising Nested Structures (ACORNS), a clustering algorithm to be fully described in Henshaw et al. (in prep.) (cf. 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 significance333For our combined ammonia data, we set the “significance” of the datapoints in ACORNS to be the peak signal-to-noise ratio of the combined NH3 (1,1) and (2,2) spectra. until a hierarchy is established. After the first cluster has been assigned in such fashion, the most significant point among unassigned dataset 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 datapoints was set to be smaller or equal to the velocity resolution of the data (0.2 $\mathrm{km~{}s^{-1}}$). Additionally, no links with line widths of twice the velocity resolution (0.4 $\mathrm{km~{}s^{-1}}$) were allowed to be established to prevent spatially close yet distinct components from being merged. In our ACORNS run, the PPV points below signal-to-noise significance of three were set to be excluded.

We obtain two sets of results from the ACORNS clustering for the two masking strategies introduced in §2.2. The strict censoring results in thirteen hierarchical clusters, seven 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 on Fig. 3, while the separated maps of the seven largest components will be presented in the following sections. The relaxed masking results in twenty-four components, of which ten are spanning a significantly large area (a comparison with the adopted method is shown in Appendix C).

3.1.1 Continuum sources within the velocity components

Throughout this study, we will relate the continuum sources discovered towards the IRDC, in particular the mass surface density cores of Butler & Tan (2012) and the 70 $\mathrm{\mu 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. Firstly, 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. Secondly, 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 on Fig. 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.

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 $V_{\mathrm{lsr}}$ 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(\alpha,\delta)=v+a\Delta\alpha+b\Delta\delta$, where the R.A. and Dec. offsets $\Delta\alpha$ and $\Delta\delta$ are confined to a circle of radius $R$ around a given set of sky coordinates $\mathbf{r}=(\alpha,\delta)$, and $(v,a,b)$ are the free parameters of the least squares fit ($v$ is the centroid velocity and $\sqrt{a^{2}+b^{2}}$ is the velocity gradient magnitude). For each velocity component identified and for every position $\mathbf{r}$, we construct a set of valid pixels as follows. First, only the pixels with a well defined $V_{\mathrm{lsr}}$ value (i.e. the relevant spectral component having a peak $\mathrm{S/N}>3$) are included in the velocity gradient analysis, and the rest are masked. Secondly, the moving least squares operates only on the values within three-beam diameter (i.e., $R=1.5\times 5.44\arcsec$) 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 $V_{\mathrm{lsr}}$ regions within radius $R$ that are not connected with position $\mathbf{r}$ with valid $V_{\mathrm{lsr}}$ pixels are masked as well. For every valid point with fitted $V_{\mathrm{lsr}}$ value, an optimal set of parameters $\mathbf{p}(\mathbf{r})=a_{\mathbf{r}},b_{\mathbf{r}}$ is obtained by a weighted least-squares fit to a collection of nearby pixels:

[TABLE]

where $\mathbf{r^{\prime}}=(\alpha+\Delta\alpha,\delta+\Delta\delta$) is a vector on the plane of sky, and $w(\mathbf{r^{\prime}})$, the weights of the fit at a given position, were set to be the fitting uncertainties on the $V_{\mathrm{lsr}}$.

Following Goodman et al. (1993), we derive the velocity gradient magnitude \|\nabla\mbox{V_{\mathrm{lsr}}}\| and orientation $\theta$ for all the velocity components. Additionally, we propagate each pair of best-fit $(a,b)$ values into uncertainties on \|\nabla\mbox{V_{\mathrm{lsr}}}\| and $\theta$. The dominant majority (96%) of the values derived are statistically significant ($\mathrm{S/N}>3$), with a mean value of uncertainties on \|\nabla\mbox{V_{\mathrm{lsr}}}\| of 0.06 $\mathrm{km~{}s^{-1}~{}pc^{-1}}$. Fig. 4 shows the $V_{\mathrm{lsr}}$ for the four largest velocity components (F2, F4, F5, and F1), overlaid with the velocity gradients derived for them.

The overall velocity gradient magnitude distribution peaks at 0.35 $\mathrm{km~{}s^{-1}~{}pc^{-1}}$ (Fig. 5, values with ${\rm S}/{\rm N}<3$ are not shown) and shows an asymmetrical tail towards 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 sub-sample 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 \|\nabla\mbox{V_{\mathrm{lsr}}}\| 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.

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 N2H+, finding typical velocity gradients of $0.8-2.3$ $\mathrm{km~{}s^{-1}~{}pc^{-1}}$. 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, NH3 (1,1) as well as with similar physical scales probed (0.07 pc from the median distance to Goodman et al. core sample), yielding the velocity gradient magnitudes of $0.3-2.5$ $\mathrm{km~{}s^{-1}~{}pc^{-1}}$, which is broadly consistent with our findings in G035.39 (Fig. 5).

3.3 Temperature maps

Following the line fitting and component separation procedures outlined in §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. 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 the 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 enables a discussion on the individual contribution of the 70 $\mathrm{\mu m}$ sources to the temperature enhancements seen in the IRDC. Fig. 6 shows a comparison between the $T_{\mathrm{kin}}$ distribution and the temperature values within one VLA beam around individual 70 $\mathrm{\mu 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_{\mathrm{kin}}$ values in the vicinity of the 70 $\mathrm{\mu m}$ sources show slight temperature enhancements. Temperatures consistently lower than the mean (by about 2 K) are found towards the neck-like structure connecting the north and the south parts of the IRDC (left subfigure on Fig. 7), which is largely devoid of protostellar sources, with the exception of one 70 $\mathrm{\mu m}$ source on the western edge of the filament, previously classified as 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 $\mathrm{\mu 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.

The protostellar source in the north-eastern part of the IRDC (F2 component, left panel on Fig. 7), classified as a low-mass dense core ($9\pm 1$ $M_{\sun}$, 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-$\mathrm{\mu 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 $\mathrm{\mu m}$ point source emission (Fig. 7). With the velocity centroid of the ammonia line of 45.1 $\mathrm{km~{}s^{-1}}$ continuously morphing with the kinematics of the rest of the cloud (Fig. 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 towards it. As shown in Sokolov et al. (2018), it coincides spatially with an extended region of subsonic turbulence, and features some of the narrowest non-thermal velocity dispersions in the whole dataset.

The velocity feature at 42.5 $\mathrm{km~{}s^{-1}}$, labelled as F1 filament in our study (expanded into a subfigure in Fig. 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 blue-shifted component shows a velocity swing of about 0.5 $\mathrm{km~{}s^{-1}}$ towards its center, where an extinction core (Butler & Tan, 2009, 2012) is located, and a 70 $\mathrm{\mu m}$ source is found next to it (Nguyen Luong et al., 2011). The angular separation between the 70 $\mathrm{\mu m}$ source and the Butler & Tan (2009, 2012) extinction core is 4.3″, which is smaller than both the beam of our ammonia observations as well as the 70 $\mathrm{\mu m}$ band of the SPIRE instrument. The two sources are thus likely to be related. The 70 $\mathrm{\mu 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 as the H6 core (Fig. 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 a 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.

In G035.39, the point-like 70 $\mathrm{\mu 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 (Fig. 7), with a mean kinetic temperature of 10.7 K within the VLA beam, and has a complex kinematic morphology (Fig. 4). Recent observations by Liu et al. (2018) find evidence for a pinched magnetic field morphology towards the core, indicating that the line of sight velocity gradients on Fig. 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 towards smaller scales (Larson, 1981), allowing for bound prestellar core formation, and, ultimately, for 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 towards 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), as well as large, filament-scale subsonic regions in Taurus, Musca, Serpens South (Tafalla & Hacar, 2015; Hacar et al., 2016; Friesen et al., 2016).

Massive star formation on the other hand, is thought to be governed by non-thermal 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 $\mathcal{M}\leq 1$). 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 numbers derivation, as well as the propagation of uncertainty on $\mathcal{M}$ matching that of Sokolov et al. (2018) are included as an appendix item in Appendix D.

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 $\mathrm{\mu 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 JVLA observations are needed to fully resolve the boundary over which 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 N2H+

mapped by  with their Plateau de Bure Interferometer (PdBI) observations of N2H+ (1-0) line transition is fully covered in our VLA mosaic. , we are well-positioned to perform a quantitative comparison of the two tracers of dense gas kinematics. Assuming that the NH3 and N2H+ molecules are tracing similar gas, the turbulent (or non-thermal) components of their line widths should be equal. Figure 9 shows a direct comparison between the two non-thermal line widths. As the N2H+ molecule is heavier than NH3, we have subtracted the thermal component from both transitions in making the figure for a more direct comparison, taking $T_{\mathrm{kin}}$ from the mean ammonia kinetic temperature for both molecules to assure the uniform comparison between the two datasets. Only the data from the commonly mapped region are shown on Fig. 9

Surprisingly, we find the non-thermal line widths derived from our VLA observations to be about $20\%$ narrower than those of the PdBI N2H+ (1-0) from Henshaw et al. (2014). This is unexpected, as Henshaw et al. (2014) data resolution are comparable to our common VLA restoring beam size of 5.44″. We also note that both datasets have also had the missing flux filled in from complementary single-dish observations. In addition, the critical densities of the NH3 (1,1) and (2,2) transitions are an order of magnitude lower than that of N2H+ (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 to trace gas also at densities significantly larger than the critical density and its abundances increases toward the core centers (Tafalla et al., 2002; Crapsi et al., 2007; Caselli et al., 2017); thus, NH3 is expected to follow material traced by N2H+ (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 modelling and varying chemical differentiation. In a detailed study of two starless cores, Tafalla et al. (2004) have found a similar discrepancy in the line widths of NH3 and N2H+, and have attributed the differences in non-thermal line widths to non-LTE effects of the N2H+ lines (cf. Caselli et al., 1995). Because the PdBI data from Henshaw et al. (2014) were fit with an isolated hyperfine component only, modelling it as a simple Gaussian component with $\tau\ll 1$, the relative broadening of N2H+ 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 N2H+ and NH3 across the IRDC, preferentially tracing regions of lower turbulence in ammonia emission. Large differences in abundance profiles of different molecules are commonly found towards low-mass starless cores (e.g., Tafalla et al., 2002; Spezzano et al., 2017) and infrared dark clouds (Feng et al., 2016). For the particular case of NH3 and N2H+, 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, b). 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 NH3 and N2H+. However, both tracers appear to follow the dense gas closely, showing strong linear correlations with mass surface density from Kainulainen & Tan (2013): for NH3 (1,1) (Pearson’s $r=0.73$), NH3 (2,2) ($r=0.71$), and N2H+ (1-0) ($r\sim 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 towards the mass condensations forming the well, and, therefore, the observed line-of-sight velocity gradient magnitudes would be enhanced towards such clumps. However, we find no correlation between the local velocity gradient magnitude and the mass surface density from Kainulainen & Tan (2013) and, furthermore, we find no significant enhancement in the velocity gradient magnitudes towards the locations of mass surface density peaks (Fig. 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 22) from Nguyen Luong et al. (2011), a smooth velocity gradient pattern is seen across the cores. Observationally, a pattern like this can be modelled 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 ${\sim}1$ $\mathrm{km~{}s^{-1}~{}pc^{-1}}$ within three VLA beams diameter: the mean value of \|\nabla\mbox{V_{\mathrm{lsr}}}\| is 1.05 $\mathrm{km~{}s^{-1}~{}pc^{-1}}$ for MDC 22 (1.01 $\mathrm{km~{}s^{-1}~{}pc^{-1}}$ for MDC 7), while the 1-sigma spread in the \|\nabla\mbox{V_{\mathrm{lsr}}}\| values is 0.15 $\mathrm{km~{}s^{-1}~{}pc^{-1}}$ (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 $\mathrm{km~{}s^{-1}~{}pc^{-1}}$ within three VLA beam diameters of MDCs 22 and 7, respectively. Taking the deconvolved FWHM sizes of the 70 $\mathrm{\mu m}$ sources from Nguyen Luong et al. (2011) for both sources (0.12 ad 0.13 pc for MDC 7 and 22, respectively), as values for $R$, we arrive at the specific angular momenta of $J/M\sim 2\times 10^{21}~{}\mathrm{cm^{2}~{}s^{-1}}$, 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 is sensitive not only to ${\sim}0.1$ pc core motions but also to gas motions on larger scales. The seemingly ordered velocity field around MDC 7 and 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 turbulent component at 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 pressurised and turbulent environment. If the origin of angular momentum of the star-forming material is related to the degree of non-thermal motions, the specific angular momentum inherited by the IRDC cores could be boosted (cf. 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 smaller degree of non-thermal 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.

In the first interferometric mapping of the entire IRDC G035.39 cloud, we have mapped the NH3 (1,1) and (2,2) inversion transitions, combining the VLA and GBT observations. 2. 2.

Our data reveal the convoluted kinematics in the cloud, bridging the cloud and core scales. We identify seven principal structures in position-position-velocity 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. 3.

Contrary to expectations from theoretical observational perspectives, we find a systematic line width difference between NH3 and N2H+, with ammonia line width being about 30% smaller, a discrepancy also found towards low-mass dense cores. 4. 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. 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 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\times 10^{21}~{}\mathrm{cm^{2}~{}s^{-1}}$, similar to that of the low-mass dense cores.

VS, JEP, and PC acknowledge the support from the European Research Council (ERC; project PALs 320620).

JCT acknowledges NASA grant 14-ADAP14-0135. This publication made use of ACORNS, a Python package used to link PPV data.

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 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 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 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 multi-scale generalization of the CLEAN algorithm (Cornwell, 2008) implemented within CASA (McMullin et al., 2007) as 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\sigma$, $3\sigma$, and, finally, down to $1\sigma$) for both NH3 (1,1) and (2,2) spectral cubes. For the first run of the algorithm, the CLEAN mask is calculated444The CLEAN mask is calculated as follows. For every channel of the spectral cube, we select an intensity threshold, and create a mask based on it. The mask then undergoes a binary closing operation with a structural element that resembles a circularized VLA beam. Additionally, all the elements of the mask that have fewer than three beams worth of pixels are flagged out to disallow noise features from propagating into the CLEAN mask. For the final, tclean run, the binary closing operation uses a three-channels-wide structural element. This helps to extend the mask from neighboring channels, making it more reliable and continuous - e.g., the emission not deemed significant enough would still be included in the CLEAN mask if its neighboring channels have emission above the threshold level. for each spectral channel based on a dirty image (with a $3\sigma$ 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\sigma$ noise level threshold.

Secondly, we feather the cleaned VLA images with the GBT data, using FEATHER task within the CASA package. We found the default settings for the dish diameter (${\sim}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 $\mathrm{km~{}s^{-1}}$ of the VLA data, and interpolated onto the spatial map used in the VLA imaging (1″ pixels in 240 by 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 (${\sim}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 multi-scale 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 two.

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 signal-to-noise ratio ($\mathrm{SNR}=3$) contours, is shown on 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 $\mathrm{SNR}>3$ closed contours off the main filament of G035.39 shown for comparison on Fig. 10.

To quantify the degree of improvement, we use the fact the (1,1) and (2,2) lines of NH3 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 due to 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 SNR cuts (i.e. on the y-axis we check for correlation of all pixels with $\mathrm{SNR}>\mathrm{x}$). While significant ($r>0.5$) linear correlation is present between (1,1) and (2,2) integrated intensities for all the SNR ranges, the correlation is consistently higher for the imaging setup we use in this work.

Appendix B Channel maps of the NH3 (2,2) line

Channel maps for the NH3 (2,2) line, presented in the same manner as the NH3 (1,1) channel maps in the main text (Fig. 1), are presented in Figure 12. Overall both the integrated intensity structure and the signal to noise levels agree with that of the NH3 (1,1) channel maps.

Appendix C Relaxed censoring PPV-structures

The strict censoring approach (§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 NH3 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 §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 Fig. 3. The velocity gradient magnitude distribution, if derived in the same manner as in §3.2, has a mean value of 0.72 $\mathrm{km~{}s^{-1}~{}pc^{-1}}$, close to 0.75 obtained in §3.2. The spread in the distributions of velocity gradient magnitudes is similar as well: with standard deviations of 0.42 and 0.45 $\mathrm{km~{}s^{-1}~{}pc^{-1}}$ for relaxed and strict censoring strategies, respectively.

Given the similarity of the derived kinematic properties, we argue that the subsequent discussion on IRDC kinematics in §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 Non-thermal line widths and propagation of uncertainties in Mach numbers

To derive the degree of non-thermal gas motions and gas turbulence in G035.39 from the ammonia data, we use the following relations. Firstly, accounting for the channel width broadening effects, we subtract the channel width from the observed full width at half-max maximum of the line:

[TABLE]

where \Delta\mbox{V_{\mathrm{obs}}} is the fitted line width from the observed line profile and \Delta\mbox{V_{\mathrm{chan}}} is taken to be the at 0.2 $\mathrm{km~{}s^{-1}}$. To find the degree of non-thermal gas motions, we then remove the component caused by the thermal broadening of the line, obtaining the non-thermal velocity dispersion from the line width (Myers et al., 1991):

[TABLE]

where \mbox{m_{\rm NH_{3}}}=17~{}m_{\rm H} is the molecular weight of the ammonia molecule.

A useful way to quantify the degree of the non-thermal motions in the gas is to compare those against the local sonic gas speed \mbox{c_{\rm s}}\equiv\sqrt{\frac{\mbox{k_{\rm B}}\mbox{T_{\mathrm{kin}}}}{\mbox{m_{\rm gas}}}}, where $m_{\rm gas}$, the average mass of the particle in the gas medium, is taken to be $2.33$ a.m.u. Following this, a sonic Mach number is computed as the ratio of the two quantities: \mathcal{M}\equiv\mbox{\sigma_{\rm nt}}/\mbox{c_{\rm s}}.

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 pyspeckit package, which provides one-sigma uncertainties on $T_{\mathrm{kin}}$ and $\sigma$. Assuming these uncertainties to be uncorrelated, we propagate the errors as follows:

[TABLE]

where the notations $\sigma_{T_{\rm kin}}$ and $\sigma_{\sigma_{\rm obs}}$ are the uncertainties on the gas kinetic temperatures and observed velocity dispersion, respectively.

Substituting appropriate terms into Equation D1, we arrive at the following expression:

[TABLE]

where the terms \mbox{\sigma_{\rm chan}}=\Delta\mbox{V_{\mathrm{chan}}}/\sqrt{8\ln{2}} and \mbox{\sigma_{\rm obs}}=\Delta\mbox{V_{\mathrm{obs}}}/\sqrt{8\ln{2}} are introduced for clarity.

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

Appendix E Relative orientation of the velocity components