Ubiquitous $\rm NH_3$ supersonic component in L1688 coherent cores

Context : Star formation takes place in cold dense cores in molecular clouds. Earlier observations have found that dense cores exhibit subsonic non-thermal velocity dispersions. In contrast, CO observations show that the ambient large-scale cloud is warmer and has supersonic velocity dispersions. Aims : We aim to study the ammonia ($\rm NH_3$) molecular line profiles with exquisite sensitivity towards the coherent cores in L1688 in order to study their kinematical properties in unprecedented detail. Methods : We used $\rm NH_3$ (1,1) and (2,2) data from the first data release (DR1) in the Green Bank Ammonia Survey (GAS). We first smoothed the data to a larger beam of 1' to obtain substantially more extended maps of velocity dispersion and kinetic temperature, compared to the DR1 maps. We then identified the coherent cores in the cloud and analysed the averaged line profiles towards the cores. Results : For the first time, we detected a faint (mean $\rm NH_3$(1,1) peak brightness $<$0.25 K in $T_{MB}$), supersonic component towards all the coherent cores in L1688. We fitted two components, one broad and one narrow, and derived the kinetic temperature and velocity dispersion of each component. The broad components towards all cores have supersonic linewidths ($\mathcal{M}_S \ge 1$). This component biases the estimate of the narrow dense core component's velocity dispersion by $\approx$28% and the kinetic temperature by $\approx$10%, on average, as compared to the results from single-component fits. Conclusions : Neglecting this ubiquitous presence of a broad component towards all coherent cores causes the typical single-component fit to overestimate the temperature and velocity dispersion. This affects the derived detailed physical structure and stability of the cores estimated from $\rm NH_3$ observations.


Introduction
Star formation takes place in dense cores, which are embedded in molecular clouds. These cores in various molecular clouds are therefore studied in detail in relation to their physical and chemical properties. Across different molecular clouds, the starforming cores are characterised by higher density and lower temperatures, as compared to the parental cloud (Myers 1983;Myers & Benson 1983;Caselli et al. 2002). When studied using molecular transitions tracing higher densities (n(H 2 ) > 10 4 cm −3 ), the cores are revealed to show subsonic levels of turbulence Kirk et al. 2007;Rosolowsky et al. A&A proofs: manuscript no. aanda also noted eight 'quiescent' cores, which they defined as local minima in non-thermal dispersion, in the observation of Serpens in N 2 H + , using BIMA (Berkeley Illinois Maryland Association). They showed that the dispersion increased outwards from the centres of these cores. Using NH 3 (1,1) observations with the Green Bank Telescope, Pineda et al. (2010) studied the transition from the surrounding gas to inside a core, for the first time in the same tracer in the B5 region in Perseus, and reported a sharp transition to coherence (a decrease in the velocity dispersion by a factor of 2 within 0.04 pc). This suggests that we can define the boundaries of coherent cores in a systematic fashion, as regions with subsonic non-thermal linewidths (see also Chen et al. 2019). It should be noted that this definition of cores does not necessarily define cores identically to those defined using continuum emission. Therefore, not all of the 'coherent cores' have continuum counterparts.
In the Green Bank Ammonia Survey (GAS; Friesen et al. 2017), star forming regions in the Gould Belt were observed using NH 3 hyperfine transitions. The first data release (DR1) includs the following four regions in nearby molecular clouds: B18 in Taurus, NGC 1333 in Perseus, L1688 in Ophiuchus, and Orion A North. Using the results from this survey and H 2 column densities derived from Herschel data, Chen et al. (2019) identified 12 coherent structures (termed 'droplet' in the paper) in L1688 and six in B18. They observed that these droplets show high gas density ( n H 2 ≈ 5 × 10 4 cm −3 , derived from their masses and radii) and near-constant, almost-thermal velocity dispersions, with a sharp velocity dispersion increase across the boundary.
This letter presents a new analysis of the GAS DR1 observations to reach unprecedented noise levels towards coherent cores and reveal a broad supersonic component that has been previously unidentified. We focus on L1688, which is part of the Ophiuchus molecular cloud (distance = 138.4±2.6 pc, Ortiz-León et al. 2018, and mass ∼980 M , Ladjelate et al. 2020) because it is one of the nearest star-forming regions and with the most extended NH 3 emission beyond cores among the GAS DR1 regions.

Ammonia maps
We use the NH 3 (1,1) and (2,2) maps from GAS DR1 (Friesen et al. 2017). The observations were carried out using the Green Bank Telescope (GBT) to map NH 3 emission in the star-forming regions in the Gould Belt with A V >7 mag using the seven-beam K-Band Focal Plane Array (KFPA) at the GBT. Observations were performed in frequency switching mode with a frequency throw of 4.11 MHz (≈ 52 km s −1 at 23.7 GHz). The spectral resolution of the data is 5.7 kHz, which corresponds to ≈ 72.1 m s −1 at 23.7 GHz (i,e, the approximate frequency of observations). The extents of the maps were selected using continuum data from Herschel or JCMT, or extinction maps derived from 2MASS (Two Micron All Sky Survey). To convert the spectra from frequency to velocity space, the rest frequencies for NH 3 (1,1) and (2,2) lines were considered to be 23.6944955 GHz and 23.7226336 GHz, respectively.
The parameter maps of L1688, which were released in DR1, are very patchy, particularly for the kinetic temperature. Therefore, to obtain a good estimate of the Mach number throughout the cloud, which requires determining the temperature, and thus, a good detection of the (2,2) line, the data were smoothed by convolving them to a beam of 1 (GBT native beam at 23 GHz is ≈31 ). The data cube was then re-gridded to avoid oversampling. The relative pixel size was kept the same as in the original GAS maps, at one-third the beam width. The median noise level achieved as a result is 39 mK in the regions of interest (see Sections 3.2 and 3.3). In comparison, the median noise in the same regions was 131 mK in the original GAS DR1 maps. 1

Line fitting
We fitted NH 3 line profiles to the data with the pyspeckit package (Ginsburg & Mirocha 2011), which uses a forward modelling approach, following the process described in Friesen et al. (2017). The range in velocity to fit is determined from the average spectrum of the entire region ( Figure F.1). Since the orthoto-para ratio in the region is not known, we only used the p-NH 3 column density in the fitting process, and did not attempt to convert it to total NH 3 column density 2 .
The model produces synthetic spectra based on initial guesses provided for the input parameters: excitation temperature (T ex ), kinetic temperature (T K ), para-ammonia column density (N(p-NH 3 )), velocity dispersion (σ v ), and line-of-sight central velocity (v LSR ) of the gas (Section 3.1 in Friesen et al. 2017). A non-linear gradient descent algorithm, MPFIT (Markwardt 2009) is then used to determine the best-fit model, and the corresponding values of the parameters. A good set of initial conditions is necessary to ensure that the non-linear least-squares fitting does not get stuck in a local minimum. The value of v LSR is critical, and, therefore, we used the first-order moments of the (1,1) line as initial guesses. The second-order moments were used as guesses for velocity dispersion. For the other parameters, we used guesses based on the GAS DR1 map of the respective parameter. The initial guesses used in our work are: log 10 (N p−NH 3 /cm −2 ) = 14, T K = 20 K, and T ex = 5 K. These numbers are within the range of the values reported in the GAS DR1 maps, and therefore, they are reasonable guesses. As a test, we checked the fits with varying guesses and determined that the exact values of the initial parameters did not affect the final results, as long as they were within a reasonable range.
In this work, we use the cold_ammonia model in the pyspeckit library, which makes the assumption that only the (1,1), (2,2) and (2,1) levels are occupied. It is also assumed that the excitation temperature, T ex , is the same for the (1,1) and (2,2) transitions, as well as their hyperfine components.

Identification of coherent cores
To study the spectra towards cores, we first define 'coherent cores' as the region with a sonic Mach number < 1 and larger than a beam (similar to Pineda et al. 2010;Chen et al. 2019) 3 . The sonic Mach number, M S , is defined as the ratio of the nonthermal velocity dispersion, σ NT (see Appendix A) to the sound 1 The GAS data use a tapered Bessel function as the gridding kernel (advocated for by Mangum et al. 2007), to achieve maximum resolution. The kernel is not strictly positive, which leads to adjacent pixels having anti-correlated noise. When smoothing over anti-correlated values, the noise level drops faster than what is expected from independent data. Therefore, we see the noise level drop faster than the inverse of the beam radius. 2 To compare N(p-NH3) to the total NH 3 column densities reported in other works, an easy conversion is to multiply it by two, with the ortho-to-para ratio assumed to be unity. speed in the medium, where c S is the one-dimensional sound speed in the gas, Here k B is the Boltzmann's constant, T K is the kinetic temperature in the region, and µ gas = 2.37 amu is the average molecular mass (Kauffmann et al. 2008). Based on this definition, we identify 12 regions in L1688 as coherent cores. They are shown in Figure 1. See Appendix B for a brief description of the cores.

Averaging spectra towards the coherent cores
We first averaged the spectra towards each coherent core defined in Section 3.2, to get a higher signal-to-noise ratio. To avoid any possible line broadening due to averaging in a region with velocity gradients, we aligned the spectra at each pixel within a core before averaging. We took the velocity at a pixel, which was determined from the single-component fit at that pixel, and using the channelShift function from module gridregion in the GAS pipeline 4 , we shifted the spectra by the corresponding number of channels. Then, we averaged the resultant spectra from all pixels inside a core, which are now essentially aligned at v=0. This averaging allowed us to reach a typical noise level of 18 mK in the average spectra, which is almost seven times better than the median noise level of 131 mK in the individual spectra towards the cores in the original GAS maps.
We checked if smoothing the data (Section 2) before averaging had any effect on the results and found that the change is not significant. See Appendix C for details.

Presence of a second component in core spectra
The top panels in Figure 2 show the average spectra in the Oph-F core, with a single-component fit (see Appendix D for the guesses used in the fit process), and corresponding residual. Oph-F is shown here as an example, the average spectra for the other coherent cores are shown in Appendix G (Figures G.1 to G.11). From the residuals, it is evident that the fit does not recover all of the flux at the positions of the hyperfine lines. In particular, the single-component fit misses the flux present in a broad component, which is ubiquitously present towards all cores.
Therefore, we added a second component to the fit (see Appendix D for the fit procedure). The resultant two-component fits and the individual components for Oph-F are shown in the bottom panel of Figure 2, and in Figures G.1 to G.11 for the other cores. The two-component fit is clearly an improvement (better χ 2 ), as seen from the respective residuals. We employed the Akaike Information Criterion (AIC, see Appendix E) to verify that for each core, considering a second component corresponds to a significant improvement in the fit to the spectra. The difference in AIC values from single-component fits to twocomponent fits, ∆ AIC = AIC 1−comp. − AIC 2−comp. , are shown in Table G.1. 4 https://github.com/GBTAmmoniaSurvey/GAS/tree/master/GAS In model comparison, the model with the lowest AIC value is the preferred one, and, therefore all regions are better fit by a two-component fit. However, for coherent cores L1688-d12 and L1688-SR3, the value of ∆ AIC is small (see Table G.1) and the average spectra fits (Figures G.8 and G.11) only show marginal improvements with the two-component fit. Therefore, we note that these two regions are better fit with a two-component model, but their fit-derived results are not as well-constrained as in the other ten coherent cores.
The kinetic temperatures, velocity dispersions, LSR velocities, and p-NH 3 column densities for the single-component fits and for each individual component in the two-component fits are shown in Table G.1. We see that one of the two components is subsonic (M S < 1) across all cores, whereas the other component has M S ≥1. We refer to the subsonic and supersonic components as 'narrow' and 'broad' components, respectively.
We note that the residuals from the two-component fit in some cores, in particular Oph-E, might suggest the presence of a third component or a non-Gaussian component. However, we do not see this towards all of the cores. Even when this occurs, the indication is not very clear : residuals comparable to the noise or not more than one or two channels wide. Therefore, we did not fit more than two components to the spectra in the cores.

Detection of a broad component in cores
A narrow and a broad component of ammonia emission has been detected in all of the coherent cores in L1688 ( Figure 1). The broad component in the cores is fainter, with a mean (1,1) peak brightness temperature (T p ) ≈0.25 K and T p,(2,2) < 0.13 K, except in Oph-A (T p,(2,2) ≈0.27 K) and Oph-C (T p,(2,2) ≈0.16 K). Because of such low intensities, this component cannot be detected from the previous ammonia maps. As a comparison, the original GAS maps had noise levels of ≈131 mK in the coherent core regions, which is not low enough for a 3σ detection of the (1,1), and in most cases, this is comparable to the T MB of the broad component in the (2,2). Therefore, using the original GAS data, one would not be able to detect the broad component, and the temperature and linewidth of the narrow core component would be affected. As can be seen in Table G.1, the mean noise level in the average spectra in the cores is ≈18 mK in this work, enabling us to comfortably fit even the faint broad component in these spectra for both (1,1) and (2,2).

Spectra in the ambient cloud outside cores
To check if the two components, which are seen towards the coherent cores, are present throughout L1688, we looked at the spectra outside of the cores. Since the spectra outside of the cores are even fainter than those inside of them, we have to consider the average spectra in a small region, away from the cores. Therefore, we defined a rectangular box (shown by the blackdashed box in Figure 1) that is representative of the ambient cloud in L1688, as : being at least two beams away from any coherent cores, to avoid contamination from the core spectra; not including, or being very close( one beam) to known protostars; and being away from the western edge of the cloud, which is affected by a strong external illumination (Habart et al. 2003).
The top panels of Figure  Top panels : Average NH 3 (1,1) and (2,2) spectra of Oph-F core, with a single-component fit. The blue-dotted line shows the median noise level in Oph-F in the original GAS data, which clearly shows that the broad component cannot be detected in the individual spectra in GAS DR1 data. Bottom panels : Same spectra, with two-component fit (green). The narrow (red) and broad (blue) components are also shown separately.
two-component fit was attempted, we did not see a narrow component (bottom panel of Figure 3). AIC indicates that the twocomponent fit is a better model, with both components being supersonic (Mach number >1). Therefore, we conclude that these are two broad components in the ambient cloud, which are separated in centroid velocity by ≈0.6 km s −1 . The two components in the spectra have equal dispersions of 0.35 km s −1 . The kinetic temperatures of these two components are also approximately equal, 16.9 ± 0.2 K and 17.0 ± 0.5 K. These values are comparable to those from the broad components in the cores, and higher than the typical temperature of the narrow components.
The Mach numbers of these two cloud components are comparable to that of the broad components seen in the cores (Table  G.1).
The residuals from the single-component fit to the ambient cloud spectra are not as suggestive as in the cores. Therefore, the single-component model can also be considered as a reasonable fit. Even then, the comparison between the ambient cloud and the broad component towards the cores, as mentioned above, remains the same. The temperature, velocity dispersion, and the Mach number in the ambient cloud, from a single-component fit, The difference is as high as ≈0.1 km s −1 in the case of SR1. On average, the single-component fit overestimates the kinetic temperature by 10% and the velocity dispersion by 28%, as compared to the narrow component. Correspondingly, the Mach number is overestimated, on average, by 32%, and the non-thermal dispersion is typically overestimated by 39%. The average kinetic temperatures and velocity dispersions inside the cores in the GAS DR1 parameter maps, which were obtained using single-component fits, are very similar to the single-component results we report. Therefore, the aforementioned comparison with narrow component results still holds, if the DR1 results (final two columns in Table G.1) are considered. On average, the GAS DR1 results show a 14% overestimate in T K and a 34% overestimate in σ v , as compared to our narrow component results.

Discussion and conclusions
For the first time, a faint, broad component has been detected towards coherent dense cores. When we consider a twocomponent fit to the spectra, it enables us to measure the velocity dispersion and kinetic temperature of the coherent core more accurately. We find that the typically-used single-component fit towards dense cores overestimates their temperatures and dispersions, on average, by 10% and 28%, respectively. 5 Weighted averages of ∆ T K and ∆σ v . The 1.2 K systematic offset in temperature in the cores derived from the single-component fit is important, as it is comparable to the temperature drop observed towards the centre of the cores when compared to the outer 1 or 2 (Crapsi et al. 2007;Harju et al. 2017). This decreased temperature of the cores has an effect on the chemistry inside the cores (Caselli & Ceccarelli 2012;Bergin & Tafalla 2007), since, at the volume densities of L1688 dense cores (∼ 10 5 cm −3 ), gas and dust are thermally coupled (Goldsmith 2001) and surface chemistry rates exponentially depend on the dust temperature (see Hasegawa et al. 1992).
Without the detection of the faint broad component, both the temperature and velocity dispersion in the cores are overestimated. This changes the estimates of the dynamical stability for the cores.
In the ambient cloud surrounding the coherent cores, we do not see the narrow component. Instead, at the position of our 'ambient cloud box' (Figure 1), we see two broad components that are ≈0.6 km s −1 apart in centroid velocity, with equal temperatures and velocity dispersions. This might suggest the presence of two molecular clouds at slightly different velocities, which are probably merging at the location of L1688. The collision between these two clouds might result in a local density increase, where the narrow features are produced following a corresponding dissipation of turbulence, thus creating the observed coherent cores with subsonic linewidths.
On the other hand, if the faint broad features seen towards the coherent core position trace the less dense and more turbulent cloud along the line of sight of the coherent core, then we can measure the relative motions between the cores and the cloud in a direct way. Accordingly, we measured the relative velocity between the narrow and broad component, δv (n−b) . Twothirds of the coherent cores show subsonic |δv (n−b) | values, which suggests that most cores are dynamically coupled to their natal cloud (Kirk et al. 2010). The standard deviation of δv (n−b) is 0.35 km s −1 , which is lower than the broad component's typical velocity dispersion (0.48 km s −1 ). Therefore, the motions of the cores are lower than the typical motions in the surrounding gas (see also Walsh et al. 2007;Kirk et al. 2007Kirk et al. , 2010. Bailey et al. (2015) have also reported similar findings in simulated observations.
The cloud velocity dispersions that we derived from NH 3 are lower than those obtained using lower density tracers (CO) in the Article number, page 5 of 15 A&A proofs: manuscript no. aanda earlier studies. Therefore, the results from lower density tracers might overestimate the local degree of turbulence in the cloud, or NH 3 based measurements provide a strict lower bound to it.
Without the required sensitivity in the observations, the faint broad component is not detected, and, therefore cannot be considered in the fit. This causes both the kinetic temperature and the velocity dispersion in the cores to be overestimated. Therefore, our results suggest that with deeper observations, we obtain better estimates of the core properties. The non-thermal component of the velocity dispersion is calculated by removing the thermal dispersion for the observed molecule (σ T ) from the total observed velocity dispersion (σ v ) : where σ chan is the contribution due to the width of the channel and the thermal component of the dispersion observed in NH 3 is with µ NH 3 = 17 amu, being the mass of the ammonia molecule.
The typical values of T K and σ v in this study are 14 K and 0.2 km s −1 , respectively (see Section 4.1). At T K =14 K, the thermal component, σ T,NH 3 is 0.082 km s −1 . When a Hann-like kernel is used in the spectrometer, then the correction factor for the measured velocity dispersion is given by (see Leroy et al. 2016;Koch et al. 2018) where ∆ is the spectral resolution and k is dependent on the Hann-like function applied. In the case of the VEGAS spectrometer used in the GAS observations, k is 0.11, which corresponds to a correction term of After applying this small correction, we obtain a typical nonthermal component of the velocity dispersion, σ NT = 0.182 km s −1 .
Apart from these previously identified cores, we also identified three more subsonic regions : L1688-SR1 (south of Oph-B1 in Motte et al. 1998), L1688-SR2 (three islands east of Oph-A), and L1688-SR3 (near the western edge of the map). Ladjelate et al. (2020) identifies a prestellar core at the position of SR1 and an unbound starless core at the position of SR3. SR1 and SR3 are associated with a local peak in NH 3 integrated intensity. SR2 does not contain any bound core in the Herschel data, nor does it show a local NH 3 peak.
Oph-B3, Oph-D and H-MM1 are also identified as L1688-c4, L1688-d11 and L1688-d9, respectively, in Chen et al. (2019). The two structures that we identify as Oph-A are likely associated with Oph A-N2 and Oph A-N6, as identified in N 2 H + by Di . The continuum cores Oph-B1 and Oph-B2 from Motte et al. (1998) are not subsonic (Friesen et al. 2009) and, therefore, they are not considered as 'coherent cores' by our definition.

Appendix C: Checking effect of smoothing the data on the final results
To check if there is any effect due to smoothing the data (see Section 2), we performed the same averaging on the original GAS data, using the velocity map published in DR1 for aligning the spectra. Since the DR1 velocity map was not as extensive, this test could not be performed for all the cores, as we need the velocity information to align the spectra. In the cores where this check was possible, the results from the new fits were not significantly different (≈1.4% change in T K and ≈4.7% change in σ v ).
In comparison, the difference in the estimates from the narrow component and from single-component fit results were 10% in T K and 28% in σ v (Section 4.4).

Appendix D: Fitting the average spectra in the cores
To fit the average spectra in cores, we followed the same process as described in Section 3.1. We used the module specfit in pyspeckit to fit single spectra. Since we were fitting the spectra towards the cores, we used, as our initial guess, slightly lower values of 12 K and 4 K for T K and T ex , respectively. As the spectra are aligned at v=0, we used v LSR =0.0 km s −1 as our initial guess. For σ v , we used the mean value inside all cores, which is σ v ≈ 0.19 km s −1 , from the extended velocity dispersion map obtained from the fit process described in Section 3.1. While fitting two components to the spectra, we added a twocomponent cold_ammonia model as a new model in the fitter. We then used this model to fit the spectra. As the residuals from single-component fits suggest the presence of a narrow and a broad component, we used two values, one from each side of the guess used in the single-component fit (0.19 km s −1 ) as our initial σ v guesses for the two components. After checking the fit with varying initial guesses, we found that small differences in T K and v LSR in the initial guesses also aid in arriving at a stable fit. Therefore, we used the following initial guesses in the twocomponent fits:

Appendix E: Model selection using AIC estimator
We fitted the averaged spectra considering models with a different number of components in ammonia, with each component being modelled as a cold_ammonia spectrum in pyspeckit. The number of parameters used in the models is five per component. To select the best model to fit a spectra, we checked the Akaike Information Criterion (AIC), which determines if the quality of the model significantly improves considering the increase in the number of parameters used. Assuming that each channel in the spectra has a constant Gaussian error, which is equal to the noise in the spectra, AIC is related to the χ 2 of the data as : where k is the number of parameters used to model the data and C is a constant, which depends on the noise in the spectra, and is given by : where N is the number of channels in the spectra. The model with the lowest AIC value is considered to be the best model. In our case, the channels are not completely independent and, therefore the actual number of independent channels is different. This effect, however, is systematic and it does not change for models with different number of components. The change in the AIC value from one model to another is the important term, and the actual AIC values are less significant.
Using Equation E.1 we can rewrite the change in AIC value from a single-component fit to a two-component fit, ∆ AIC , as where ∆k is the change in number of parameters and ∆χ 2 is the change in χ 2 , from a single-component fit to a two-component fit. Each model component uses five independent parameters and, therefore ∆k = −5. Also, ∆χ 2 can be divided into two groups: channels with and without emission, or ∆χ 2 = (∆χ 2 ) with emission + (∆χ 2 ) without emission .

(E.4)
Outside the emission region the models are identical and, therefore, where, n line is the number of channels with emission, σ is the noise in the spectra, and rms emission is the residual root mean square (rms) of the model in the emission region. The comparison of the σ emission and the noise, σ, is taken as a measure of determining the goodness of the fit in line-fitting software, such as CLASS in GILDAS 6 . Then, in our comparison of single-and two-component fits, we have ∆ AIC = −10 + n line σ 2 rms 2 emission, 1−comp − rms 2 emission, 2−comp . (E.6) We note thatn line depends on the width of the spectra and, therefore, it is different for each coherent core. From Equation E.6, we see that for a low variation in the AIC, the improvement in the residuals of the fit is also small. We find that for the two coherent cores with low ∆ AIC (<30), d12 and SR3, the improvement in the rms emission is less than 10% of the noise in the respective spectra. For the ambient cloud spectra, with ∆ AIC =94, rms emission improves by 10% of the noise level, from a singlecomponent to a two-component fit. For every other region, we see an improvement in rms emission by >15%.
Owing to our limited number of data points, it is not possible to definitively determine what change in AIC represents a significant improvement in the fit upon the addition of the second component. From the corresponding improvement in rms emission , we note that in our current work, only the regions with ∆ AIC > 150 show a significant improvement in the fit.