Connecting the Scales: Large Area High-resolution Ammonia Mapping of NGC 1333

TL;DR

This study maps dense gas in NGC 1333 using ammonia observations, revealing a large-scale velocity gradient and potential shock-induced star formation, combining interferometric and single-dish data for high-resolution insights.

Contribution

It provides the first large-area, high-resolution ammonia maps of NGC 1333, integrating VLA and GBT data to analyze gas dynamics and star formation processes.

Findings

Identified a large-scale velocity gradient in NGC 1333.

Detected high velocity dispersion indicating shock layers.

Correlated ammonia features with star-forming regions.

Abstract

We use NH3 inversion transitions to trace the dense gas in the NGC 1333 region of the Perseus molecular cloud. NH3(1,1) and NH3(2,2) maps covering an area of 102 square arcminutes at an angular resolution of ~3.7" are produced by combining VLA interferometric observations with GBT single dish maps. The combined maps have a spectral resolution of 0.14 km/s and a sensitivity of 4 mJy/beam. We produce integrated intensity maps, peak intensity maps and dispersion maps of NH3(1,1) and NH3(2,2) and a line-of-sight velocity map of NH3(1,1). These are used to derive the optical depth for the NH3(1,1) main component, the excitation temperature of NH3(1,1), and the rotational temperature, kinetic temperature and column density of NH3 over the mapped area. We compare these observations with the CARMA J=1-0 observations of N2H+ and H13CO+ and conclude that they all trace the same material in these dense star forming regions. From the NH3(1,1) velocity map, we find that a velocity gradient ridge extends in an arc across the entire southern part of NGC 1333. We propose that a large scale turbulent cell is colliding with the cloud, which could result in the formation of a layer of compressed gas. This region along the velocity gradient ridge is dotted with Class 0/I YSOs, that could have formed from local overdensities in the compressed gas leading to gravitational instabilities. The NH3(1,1) velocity dispersion map also has relatively high values along this region, thereby substantiating the shock layer argument.

Figures

Tables

Table 1. Table 1: Observed Molecular Lines. The transtion frequencies are obtained from the CDMS (Müller et al., 2005 ) .

Molecular Species	Transition	Frequency (MHz)	Channel Width (kHz)	Frequency Range (MHz)
NH3	(J,K)±=(1,1)+-(1,1)-	23694.4955	3.9	5.574
NH3	(J,K)±=(2,2)+-(2,2)-	23722.6333	3.9	5.574
NH3	(J,K)±=(3,3)+-(3,3)-	23870.1292	3.9	1.997
NH3	(J,K)±=(4,4)+-(4,4)-	24139.4163	3.9	1.997
15NH3	(J,K)±=(1,1)+-(1,1)-	22624.9295	7.8	1.997
HNCO	J(Ka,Kc)=1(0,1)-0(0,0)	21981.5726	7.8	1.997
H2O	J(Ka,Kc)=6(1,6)-5(2,3)	22235.0798	7.8	1.997
CCS	N=1-2, J=2-1	22344.0308	7.8	1.997
CC34S	N=1-2, J=2-1	21930.4860	7.8	1.997

Table 2. Table 2: Properties of the original and combined NGC 1333 images

Images	u-v range	Channel Widtha	Synthesized Beam	Sensitivity
VLA	25 - 1000 m	3.9 kHz	3.76′′$\times$ 3.34′′	0.007 Jy/beam
GBT	0 - 100 m	5.7 kHz	31.8′′$\times$ 31.8′′	0.083 Jy/beam
Combined	0 - 1000 m	11.4 kHz	3.93′′$\times$ 3.41′′	0.004 Jy/beam

Equations

$$S(v)=\sum_{n=1}^{5}a_{n}\medop\sum_{i}w_{i}\exp{\bigg{(}-\frac{v-v_{los}-c\Delta f_{i}/f_{0}}{2{\sigma_{v}}^{2}}\bigg{)}}$$

$$\frac{T_R(m)}{T_R(s)}=\frac{1-e^{-{\tau }_m}}{1-e^{-\alpha {\tau }_m}}$$

$$\frac{n_{+}}{n_{-}}=\frac{g_{+}}{g_{-}}{e}^{-h\nu /k{T}_{ex}}$$

$$T_{b(1,1,m)}=\Big{(}\frac{1}{\mathcal{J}_{\nu}(T_{ex})-\mathcal{J}_{\nu}(T_{bg})}\Big{)}[1-e^{-\tau_{(1,1,m)}}]$$

$$\frac{N_{(2,2)}}{N_{(1,1)}}=\frac{5}{3}{e}^{-41.0/T_{rot}}$$

$$T_{rot}=-41.0\bigg{\{}\ln{\Big{[}-\frac{0.283\sigma_{v(2,2)}}{\tau_{(1,1,m)}\sigma_{v(1,1)}}\ln{\Big{(}1-\frac{T_{R(2,2,m)}}{T_{R(1,1,m)}}(1-e^{-\tau_{(1,1,m)}})\Big{)}}\Big{]}}\bigg{\}}$$

$$T_{rot}=T_{K}\Big{\{}1+\frac{T_{K}}{41.5}\ln\big{(}1+0.6e^{-15.7/T_{K}}\big{)}\Big{\}}$$

$$N_{(1,1),+}=\frac{3hJ(J+1)}{8{\pi }^3{\mu }^2{K}^2}\frac{1}{(e^{h\nu /k{T}_{ex(1,1)}}-1)}\int {\tau }_{(1,1)}dv$$

$$\frac{N_{(1,1),+}}{N_{(1,1),-}}=e^{-h\nu /k{T}_{ex(1,1)}}$$

$$N_{(1,1)}=N_{(1,1),-}+N_{(1,1),+}=N_{(1,1),+}(e^{h\nu /k{T}_{ex(1,1)}}+1)$$

$$N_{tot}=\frac{N_{(1,1)}{Q}_{rot}}{g_J{g}_K{g}_I}{e}^{E_{(1,1),+}/k{T}_{rot}}$$

$$N_{tot}=\frac{3hJ(J+1)}{8{\pi }^3{\mu }^2{K}^2}\frac{Q_{rot}}{g_J{g}_K{g}_I}{e}^{E_{(1,1),+}/k{T}_{rot}}\frac{e^{h\nu /k{T}_{ex(1,1)}}+1}{e^{h\nu /k{T}_{ex(1,1)}}-1}\int {\tau }_{(1,1)}dv$$

$$N_{tot}=\frac{3hJ(J+1)}{8\pi^{3}\mu^{2}K^{2}}\frac{Q_{rot}}{g_{J}g_{K}g_{I}}e^{E_{(1,1),+}/kT_{rot}}\bigg{[}\frac{e^{h\nu/kT_{ex(1,1)}}+1}{e^{h\nu/kT_{ex(1,1)}}-1}\bigg{]}\frac{1}{R_{i}}\bigg{[}\frac{\int T_{R(1,1,m)}dv}{\mathcal{J}_{\nu}(T_{ex})-\mathcal{J}_{\nu}(T_{bg})}\bigg{]}\frac{\tau_{(1,1,m)}}{1-e^{-\tau_{(1,1,m)}}}$$

$$N_{tot}=4.51\times 10^{12}{Q}_{rot}{e}^{24.35/T_{rot}}\frac{1+e^{1.137/T_{ex(1,1)}}}{1.517-e^{1.137/T_{ex(1,1)}}}\frac{{\tau }_{(1,1,m)}}{1-e^{-{\tau }_{(1,1,m)}}}\int T_{R(1,1,m)}dvc{m}^{-2}$$

$$Q_{rot}=\sum _J\sum _K(2J+1)g_K{g}_I{e}^{-E_{JK}/k{T}_{rot}}$$

Full text

ISM: clouds, ISM: kinematics and dynamics, ISM: molecules, ISM: structure, stars: formation

††journal: ApJ

1 Introduction

Maps spanning large to small scales are fundamentally important to understand star formation in molecular clouds. The relative importance of processes such as turbulence, magnetic fields, gravity and chemical evolution vary at the different scales, all contributing to the birth of a protostar. In addition, feedback from star formation also affects the molecular cloud environment over a range of distances. Large area gas and dust surveys of molecular clouds ranging from parsecs to about 1000 AU are required to get the complete picture of the formation and evolution of dense structures that eventually form stars.

In going from the less dense regions of the clouds to denser filaments and on to cores, the transitions in the morphological and kinematic properties provide important clues about the underlying processes. The scaling of velocity dispersions with cloud sizes indicates that compressible turbulence plays a key role in molecular clouds (Larson, 1981; Boldyrev et al., 2002; Hennebelle & Falgarone, 2012). Turbulence can affect star formation both by supporting a cloud against gravitational infall on large scales and by forming local density fluctuations that may give rise to self-gravitating cores at small scales (Ballesteros-Paredes et al., 1999; Padoan & Nordlund, 2002; Mac Low & Klessen, 2004). Filament-like structures, that are prevalent in molecular clouds (Mizuno et al., 1995; André et al., 2010) harbor many star forming cores. Near protostellar cores, there is increasing evidence of a sharp transition from turbulent to thermal line widths (Pineda et al., 2010; Seo et al., 2015). This behavior suggests that dense cores may form as pressure-confined structures within filaments and evolve to gravitationally bound cores before undergoing collapse to form a protostar (Pattle et al., 2015; Kirk et al., 2017). Magnetic fields within the cloud can also affect core fragmentation and collapse by providing pressure support (Chen & Ostriker, 2012). By comparing the kinematics of neutral and ionic species in the gas (Flower, 2000), it is possible to determine the presence and significance of magnetic fields.

Radiation feedback from the embedded protostars also plays an important role in suppressing fragmentation of cores (Krumholz, 2006). Numerical models show that protostellar outflows and wind from intermediate-mass stars can drive turbulence in the molecular clouds (Carroll et al., 2009; Offner & Arce, 2015). This behavior is suggested by observations of cavities and shells in Perseus as well (Quillen et al., 2005; Arce et al., 2011). Gas temperature maps that are sensitive from large to small scales can establish the importance of the feedback processes involved in star forming regions (Offner et al., 2009; Foster et al., 2009).

In this paper, we study the NGC 1333 region in Perseus using high angular resolution Very Large Array (VLA) observations of NH3 inversion transitions. NH3 is an abundant molecule that traces gas of density greater than $10^{4}$ cm*-3* (Shirley, 2015) and is a late-depleter (Bergin & Langer, 1997). Nitrogen containing molecules like NH3 and N2H*+* remain in the gas phase longer than carbon and oxygen-containing molecules such as HCO*+, HCN, and CO. Hence, it is particularly suitable to study cores (Rosolowsky et al., 2008) and dense regions of the molecular clouds that are expected to participate in star formation in the future. Most previous observations using NH3* involved small target areas around cores. The Green Bank Ammonia Survey (GAS) (Friesen et al., 2017) is the first large-scale NH3 survey of all the major clouds in the Gould Belt. The 32*′′* beam of the GBT at 23 GHz, however, is insufficient to resolve the inner structure of cores. The VLA observations reported here make it possible to probe the regions close to the cores at a resolution of 1025 AU and study how the emission is connected to the large scale emission over the molecular cloud.

The NGC 1333 region in Perseus is a reflection nebula at a distance of 299 $\pm$ 17 pc (Zucker et al., 2018). It is known for its active low-to-intermediate mass clustered star formation (Hatchell et al., 2007; Walawender et al., 2008). There exists a wealth of information for this region, which is rich in sub-mm cores (Dodds et al., 2015), YSOs (Foster et al., 2009; Johnstone et al., 2010), outflows (Plunkett et al., 2013), Herbig-Haro objects (Bally et al., 1996) and masers (Lyo et al., 2014). The Spitzer cores to disks (c2d) (Jørgensen et al., 2006; Young et al., 2015) legacy survey, the Herschel PACS and SPIRE images (André et al., 2010), the James Clerk Maxwell Telescope (JCMT) Gould Belt survey (Chen et al., 2016) and the CARMA CLASSy observations (Storm et al., 2014; Dhabal et al., 2018) all provide complementary data to track the YSOs and put the VLA data in larger context.

The layout of the paper is as follows. The VLA observational set-up is discussed in Section 2. In Section 3, we discuss how the interferometric visibility data is combined with single-dish image data from the Green Bank Telescope (GBT) to produce maps that are sensitive to large-scale structure. Using the VLA NH3 (1,1) and NH3 (2,2) observations, in Section 4 we obtain integrated intensities, peak intensities, line-of-sight velocities and velocity dispersions. In Section 5, the primary results are used to derive optical depths, excitation temperatures, rotational temperatures, kinetic temperatures and column densities. On comparing with the CARMA N2H*+* maps and the dust maps, we establish the morphology and kinematics of the clouds in Section 6. We also discuss some of these results in the context of filament formation and their relation to star formation. The main results are summarized in Section 7. Appendix A has a brief analysis of the sensitivity of a weighting parameter involved while combining single-dish and interferometric data.

2 VLA Observations of NGC 1333

The K-band (18-26 GHz) observations were carried out using the VLA array of 25-m antennas in D configuration providing a resolution of about 4*′′*. The NGC 1333 region is covered by 87 pointings to cover an area of 102 square arcminutes. They were observed in multiple sessions, with 8-12 minutes integration time on each pointing. The observations were carried out from March 2013 to May 2013.

The correlator was configured to observe the NH3 (1,1), NH3 (2,2), NH3 (3,3) and NH3 (4,4) lines with a channel width of 3.9 kHz (corresponding to about 0.049 km/s). The maximum sub-band bandwidth of 4 MHz was not sufficient to cover all the satellites of the ammonia inversion transitions. So two partially overlapping sub-bands were used for NH3 (1,1) and NH3 (2,2). The relative intensities of the satellites are negligible for NH3 (3,3) and NH3 (4,4). For these lines, we used a single sub-band covering only the main transition line. The CC34S N=1-2, J=2-1, HNCO 1(0,1)-0(0,0), H2O maser 6(1,6)-5(2,3), CCS N=1-2, J=2-1 and 15NH3 (1,1) lines were also available to be observed in the same configuration (see Table 1). A channel width of 7.8 kHz was used for these observations to improve the signal-to-noise ratio (SNR), anticipating their relatively low brightness temperatures. Two continuum bands of width 128 MHz were also observed at 22.3 GHz and 23.9 GHz.

Flagging and calibration were carried out using the CASA EVLA pipeline. J0336+3218, 3C84 and 3C147 were used as phase, bandpass and flux calibrators, respectively. The calibrated visibility data in the measurement sets were cleaned using the CASA tclean routine. The synthesized beam over the channels has a median size of 3.76*′′* $\times$ 3.34*′′. For the NH3* lines, the mean value of the spectral channel RMS noise is 6.8 mJy/beam at the native channel width. Channel binning was carried out to improve the sensitivity and look for detections in the rarer species, at the cost of velocity resolution. NH3 (1,1) and NH3 (2,2) were detected over large areas in the observed regions. NH3 (3,3) was only detected at the outflows in the NGC 1333 region. Seven 22 GHz H2O masers were also visually identified from the corresponding position-position-velocity data cube. In the remaining five molecular species transitions, no emission was detected above the noise.

3 Combination of Interferometric and Single Dish Data

The maps produced from the VLA data alone have strong negative sidelobes, which is typical of maps produced from interferometric visibility data due to the lack of small u-v spacings (Kogan, 2000). In our VLA data, the minimum u-v spacing is 25 m. Consequently, in many areas of the map the negative sidelobes corresponding to strong emission elsewhere suppress the actual flux levels and spatial variations even though the interferometer is sensitive to small scale structure. This issue can be remedied by using complementary single dish observations (Vogel et al., 1984). The NGC 1333 region has been observed in NH3 (1,1), NH3 (2,2) and NH3 (3,3) with the Green Bank Telescope (GBT) as part of the Green Bank Ammonia survey (GAS) (Friesen et al., 2017) and the data are publicly available (Friesen, 2017). These data were combined with the high resolution VLA maps to make the maps sensitive to large scale structure as well. In the rest of the paper, we only discuss the NH3 inversion transition observations in the NGC 1333 region.

The strategy described in Koda et al. (2011) is adopted to combine the interferometric and single dish data. This process involves converting the single dish data to visibilities to fill the inner part of the u-v plane, followed by weighting these visibilities with respect to the interferometric visibilities. This combined visibility dataset is then inverse Fourier transformed and cleaned to produce the final image. There are multiple ways to combine the data from single dish telescopes and interferometers. Other methods involve combining the data in the image domain (Helfer et al., 2003), or further Fourier transforming the image domain data and then combining them (Weiß et al., 2001) (as is done in the CASA feather routine). Since inverse Fourier transforming the visibilities to produce the interferometric image is a non-linear process, combining them after this step leads to loss of information from the interferometric component. Thus we combine the data in the visibility domain.

The GBT single dish maps have a spectral resolution of 5.7 kHz ($\sim$ 0.07 km/s at 23.7 GHz). We binned the channels by a factor of 2 to improve the SNR, while retaining 5 channels to cover the typical FWHM linewidth of 0.7 km/s. The resulting channel width is about 3 times the VLA channel width. A total of 400 channels cover all the NH3 (1,1) and NH3 (2,2) hyperfine components. We converted the map from K to Jy/beam using a factor of 0.69 taking into account the GBT beam size of 31.8*′′* $\times$ 31.8*′′*. The mean RMS of the resulting SD map is 0.059 Jy/beam.

After extracting the relevant region of the map, we generate fake visibility data from the single dish data cubes using the tp2vis program (Koda et al., 2011). This involves sampling the GBT single dish map at the VLA pointings with a 25 m primary beam (same as the VLA beam) and with u-v distances up to 100 m (equivalent to the GBT diameter). The single dish visibilities are inverse Fourier transformed to check how the resulting map compares with the input single dish map. The flux levels are within 10% of original map, and there are no noticeable differences in the intensity distributions.

The weight densities (weights/GHz/pointing/arcmin2) of the single dish visibilities are compared with those of the interferometric visibilities. Accordingly, the single dish weights are scaled to roughly match them (Figure 1). This ensures that the sensitivities transition smoothly from the larger scales to the smaller scales. Further analysis of the effect of weight scaling on the resulting maps is carried out in Appendix A. The combined set of visibilities is now used to generate the dirty map and cleaned using the qac_clean routine that uses the CASA tclean routine (Teuben et al., 2018). We use the multiscale deconvolver with scales of 0*′′, 4′′* and 12*′′*. The cleaning proceeds by a set of major cycles transforming between visibility and image domains while progressively updating the model and residual images. Within each major cycle, an inner loop of minor cycles cleans the map in the image domain. We needed to induce more major cycles without cleaning too deep in the corresponding minor cycles to avoid picking up interferometric artifacts. So a ‘cycleniter’ parameter of 1000 was used. Based on the noise RMS of 0.004 Jy/beam, we also use a ‘threshold’ parameter value of 0.02 Jy/beam, which decides the minimum target peak flux at the end of each minor cycle for each channel. Typically, the channels containing signal are cleaned in 10-15 major cycles. In regions having strong emission, the flux remaining in the residuals is less than 10% of the total flux.

The final maps thus obtained do not have negative sidelobes and show the high resolution features, while being sensitive to large scale structure. The synthesized beam is 3.93*′′* $\times$ 3.41*′′. The flux values are compared to the interferometric maps obtained from the VLA data alone in Figure 2. Typically, the average flux in regions of strong emission in the VLA maps is about 20% lesser than the corresponding flux in the combined maps, but the local spatial variations of the flux are very similar. The maps are also smoothed to the GBT beam size and compared to the original GBT map (Figure 3). On binning the combined map to the GBT-only map pixel size, the total flux matches to about 1.5%, while the variability in scales of the GBT beam is about 5%. Since the NH3* (2,2) emission traces warmer gas, it is less extended, and consequently the flux levels in the combined map are better matched with the interferometric map fluxes than for NH3 (1,1).

4 Results

The NH3 inversion transitions involve quantum tunneling of the N nuclei through the plane of the H nuclei in the same rotational quantum state (J,K). The inversion transitions however have hyperfine components due to electric quadrupole coupling (quantum number denoted by F1) and magnetic dipole coupling (quantum number denoted by F). In the case of NH3 (1,1), there are 18 hyperfine components that are bunched into 5 groups of 2, 3, 8, 3 and 2 lines based on the F1 transition. Their theoretical relative intensities are 0.111, 0.139, 0.5, 0.139 and 0.111, respectively. Similarly, the NH3 (2,2) transition has 24 hyperfine components also divided into 5 groups of 3, 3, 12, 3 and 3 lines, having relative intensities 0.05, 0.052, 0.796, 0.052 and 0.05, respectively (Ho & Townes, 1983; Mangum & Shirley, 2015). In both cases, the central line is the main component, and the other four components are called satellite components.

4.1 Ammonia (1,1) Results

NH3 (1,1) main and satellite lines are detected over large areas in the NGC 1333 region. The relative intensities of the hyperfine components, however, vary in different regions from their theoretical relative intensities (see top two panels of Figure 4). In many regions, we also note the presence of two velocity components in the same line-of-sight. The velocity peaks are separated by as much as 1.5 km/s in some locations, indicating that the two velocities belong to independent components (for example, see third panel of Figure 4). Integrated intensity maps are generated from the position-position-velocity data cubes, by using channels near the line centers containing signal, and clipping them at 3 $\sigma$ based on the RMS noise of each channel map. Figure 5 shows the integrated intensity map for the main NH3 (1,1) component alone. We have roughly divided the map into 11 regions based on their location in the map or the main embedded protostellar source(s) in the region.

IRAS 4, IRAS 2 and SVS 13 are well-known YSO systems, each of which has been resolved into multiple sources (Walawender et al., 2008), and are all located in the central cloud core. There is a known Herbig-Haro object HH 12 in the northern region, and it is located at the apex of the two filaments (NW and NE filaments). To the west of the HH 12 region, there is a region of modest emission (NW region), but no protostars are detected there. There is a similar relatively quiescent area (West region), about 1*′* west of the IRAS 2 region. The IRAS 7 protostellar cluster region is located to the east of the NE filament. Towards the south, there is a narrow ‘u’ shaped filament connecting the IRAS 2 and the IRAS 4 regions. At the bottom of the ‘u’ is a known Class 0 YSO source – SK 1 (Dionatos & Güdel, 2017). There is another narrow short filament in this area extending towards the south-west from the IRAS 2 region. The cloud core near IRAS 4 extends to a wide filament in the south-east, which has been studied in Dhabal et al. (2018) as part of the CLASSy-II regions.

Spectral line fitting is carried out to obtain the NH3 velocity structure of NGC 1333. Since we noted that the relative intensities of the satellites do not match theoretical predictions in many areas, we allow the peak intensities for each of the five NH3 (1,1) components to vary. The magnetic dipole based hyperfine splitting within each of the 5 groups also have spreads ranging between 11 kHz and 56 kHz. As some of them are resolvable by our channel width of 11.4 kHz, we fit the spectra to a model having all 18 hyperfine components. Even though we vary the relative intensities ($a_{n}$) of the five groups of hyperfine components ($n$ = 1,2,3,4,5), we assume that the relative intensities ($w_{i}$) of the components (denoted by $i$) within each of the 5 groups to be based on the theoretical values (as given in Table 15 of Mangum & Shirley, 2015). We assume a single line-of-sight velocity $v_{los}$ and a Gaussian line shape with a single dispersion value $\sigma_{v}$ for all the 18 components. Thus we fit the spectrum at each pixel as a function of velocity $S(v)$ using the following equation:

[TABLE]

where $c$ is the speed of light, $\Delta f_{i}$ is the frequency offset of the hyperfine component relative to the line center frequency $f_{0}$ = 23.6944955 GHz. We solve for $a_{1}$, $a_{2}$, $a_{3}$, $a_{4}$, $a_{5}$, $v_{los}$ and $\sigma_{v}$. We mask the pixels in which none of the channels have more than 4 $\sigma$ detections (noise RMS $\sigma$ = 0.004 Jy/beam). Least-squares fits are obtained after providing reasonable limits to the parameter search spaces. Pixels for which the fit solution involves any of the parameters reaching the limits are considered to be bad fits and are excluded. The pixels containing emission from outflows are identified by their unusually wide Gaussian fits and are also omitted. The results are shown in Figures 6, 7, 8 and 9. The intensity map distributions can be considered to be representative of the brightness temperature distributions.

For all these maps, we calculate a typical uncertainty, which is reported in the respective figures. For the peak intensity maps and velocity maps, while least-squares fitting the spectra for each pixel, we get the variance of each parameter estimate by diagonalizing the corresponding covariance matrices. The median value of the variances of all the included pixels is reported as the uncertainty. For the integrated intensity maps, the uncertainty is derived from the combination of the peak intensity and velocity dispersion values and uncertainties.

The peak intensity map for the main NH3 (1,1) component has a median value of 30 mJy/beam. This value corresponds to a brightness temperature of 4.9 K. The maximum brightness temperature goes up to as much as 20 K near the IRAS 2 sources. The four satellite components $a_{1}$, $a_{2}$, $a_{4}$ and $a_{5}$ have median brightness temperatures of 2.4 K, 2.0 K, 2.1 K and 2.1 K, respectively. The ratios between the different components are not consistent throughout the map, but typically indicate an optical depth of greater than 1 for the main component, while the satellite components can be considered to be optically thin.

On comparing the main component intensity map (Figure 6) to the integrated intensity map (Figure 5), we find many differences in the maximum intensity regions. The integrated intensity map has greater intensities near the Class 0/I YSOs, but the peaks of the brightness temperatures are well distributed throughout, including the filaments and regions having low integrated emission. This difference is particularly evident in the south west quadrant of the map. The SK 1 region, the West region and the southern part of the IRAS 2 region all have greater brightness temperatures compared to the northern part of the IRAS 2 region. Yet, the IRAS 2 northern region has greater integrated intensities.

Correspondingly, we see in Figure 9 that the velocity dispersions are greater in the regions containing the protostars (although they don’t always peak at the YSO locations, but around them). They are particularly high near the boundaries of the SVS 13 region going up to as much as 0.8 km/s compared to the nominal values of about 0.3 km/s in the filaments. Some of the high dispersions correspond to known locations of outflows as well.

In some regions, we see sharp boundaries in the dispersion map where areas having different estimates of dispersion are adjacent to each other. The sharp boundaries are caused by the presence of multiple velocity components of comparable intensity. In such cases, depending on the relative intensities, the spectral line fitting, which fits a single component, may select one of the components or fit both of them with a Gaussian having a greater width. These multiple velocity components also tend to be located near the protostellar sources. In addition, the SE filament has a distinct narrow region of higher dispersion as well. This feature is caused by the presence of the two parallel velocity coherent sub-filaments partially overlapping in the line-of-sight (Dhabal et al., 2018). In the narrow region, the two components, which have comparable intensity get fit by a Gaussian with a greater dispersion value. These multiple velocity component spectra also affect the velocity map and the intensity maps, but the effects are less drastic in those maps.

The intensities of the satellite components also show variations, even though they have a good positive correlation of about $0.8-0.9$ to each other and to the main component. In the ideal case of no anomalies, the correlation is expected to be 1. The F1=0-1 satellite component has greater mean intensities over most of the regions than the other satellites, while the F1=2-1 and the F1=1-0 satellites have relatively lesser mean intensities of all the satellites. The F1=1-0 emission has lesser contrast between the strong emission regions and the surroundings. A mechanism to explain the deviations from the theoretical relative intensities of the satellites was suggested by Stutzki & Winnewisser (1985). In their paper, the hyperfine anomalies are caused by selective trapping of infrared photons in the (J,K) = (2,1) to (1,1) transition. They are generally associated with maser activity and are postulated to be associated to star forming clumps of high density. A suggested indicator of this effect involves the ratio of the intensities of the outer satellites or the inner satellites. Typically $T_{B}$(F1=1-0)/$T_{B}$(F1=0-1) and $T_{B}$(F1=2-1)/$T_{B}$(F1=1-2) are both expected to be less than one in the selective trapping based anomalous regions (Stutzki et al., 1984).

We generated maps (not shown) of these two intensity ratios $O$ = $T_{B}$(F1=1-0)/$T_{B}$(F1=0-1) and $I$ = $T_{B}$(F1=2-1)/$T_{B}$(F1=1-2). The anomalous regions were defined to be ones where $O$ is less than 0.8 or greater than 1.2. By this definition, 55% of the unmasked pixels are found to be anomalous. We found that the anomalous regions are spread throughout the map area, but are in lesser concentration in regions of high NH3 (1,1) integrated intensity. The lower ratio ($O<0.8$) regions comprise 67% of the anomalous pixels, while the higher ratio regions ($O>1.2$) are concentrated mainly near the IRAS 4 region. Also we found no correlation between $O$ and $I$ (Spearman correlation coefficients between -0.1 and +0.1), even after selecting only the lower ratio anomalous pixels. Thus, the selective trapping mechanism fails to account for the relative intensity anomalies in the NGC 1333 region.

The line-of-sight velocity map shows rich variation throughout the map with values ranging from 6.25 km/s to 9.25 km/s. The velocity gradient in the SE filament (Dhabal et al., 2018), extends to the IRAS 4 region and further north to the SVS 13 region. The west part of the SVS 13 region has velocities varying monotonically by 2.0 km/s in the north-south direction. Similar large gradients are present in the northern part of the HH 12 region (in the east-west direction), and immediately to the west of the IRAS 4 sources. The IRAS 2 region also has a large magnitude velocity gradient in the SE-NW direction, but it is not monotonic. Most of these large gradients are caused by a combination of (a) multiple velocity-coherent components with varying intensities across the field of view and (b) a gradient across a single velocity-coherent component.

The NE filament has a single velocity coherent component in most of the northern half and it has a gradient across it which changes direction once. Going from east to west, the gradient changes from positive to negative near the ridge. In the southern half, again there is a gradient across the filament which changes from negative in the east to positive in the west, but this seems to be caused by multiple velocity components in the line of sight. The velocity dispersion map of this region shows a sharp contrast region along the filament ridge. The NW filament, SK 1 region, the West region and the IRAS 2 have relatively lesser variations of about 0.5 km/s. The NW region has the least line-of-sight velocity variations among all fields.

4.2 Ammonia (2,2) Results

For the integrated intensity map of NH3 (2,2) (Figure 10, left), we used a lower clip value of 2 times the RMS noise to cover a greater area of the map, that would be useful for the analyses in the following section. The NH3 (2,2) detections are spread over a relatively lesser area than (1,1) and they peak closer to the Class 0/I YSOs compared to the NH3 (1,1) integrated intensity emission. An outflow (from IRAS 2A) between the IRAS 2 region and the West region has NH3 (2,2) emission but is absent in the NH3 (1,1) map.

The satellite components of NH3 (2,2) have very low relative intensities and are indistinguishable from the noise for most regions (see Figure 4, bottom panel). We fit the central component of the NH3 (2,2) spectra with 12 hyperfine components, using their theoretical weight ratios. We fixed the line-of-sight velocities from the corresponding values of the NH3 (1,1) map to ensure fitting for the same velocity component as in NH3 (1,1). In single component regions, solving for the line-of-sight velocity gives the same result in both NH3 transitions, but this is not always followed in regions having multiple components in the line of sight. Hence, we solved for the velocity dispersion and a single peak intensity for the main NH3 (2,2) component. The peak intensity map is shown in Figure 10, right. The median value is 9 mJy/beam corresponding to a brightness temperature of 1.5 K. The maximum brightness temperature is 7.5 K near the protostellar sources in the SVS 13 region. For the (2,2) transition, there is a better match between the integrated intensity map and the peak intensity map. Only the IRAS 2 region shows a reversal in the areas of maximum intensity for the two maps, i.e., in this region the areas having greater NH3 (2,2) integrated intensity have lesser NH3 (2,2) velocity dispersions and vice-versa.

5 Analysis

5.1 Optical Depth

The optical depths can be estimated using the ratios of the main and satellite components of NH3 (1,1) assuming the hyperfine levels to be in local thermodynamic equilibrium (LTE). In this approximation, the optical depth ratios of the main to satellite components are equal to their respective theoretical intensity ratios ($\alpha$). The optical depth at any of the NH3 inversion transitions’ line centers $\tau_{m}$ can be calculated by solving for the equation (Mangum & Shirley, 2015):

[TABLE]

where $T_{R}(m)$ and $T_{R}(s)$ are the main and satellite component source radiation temperatures. Their ratio is equal to the ratio of the respective relative intensities $a_{m}$ and $a_{s}$ obtained from the spectral line fitting. For either of the outer satellites of NH3 (1,1), the theoretical value of $\alpha$ is 2/9. For the inner satellites, $\alpha$ is 5/18. For the NGC 1333 region, the LTE assumption is not accurate based on the anomalies presented in Section 4.1. So we investigated a modified LTE assumption using two satellite components having the same theoretical $\alpha$ value. In this method, we calculated the optical depth for the main component separately using the mean outer satellite intensities $\tau_{m,outer}$ and using the mean inner satellite intensities $\tau_{m,inner}$. The mean was taken to reduce the effect of the intensity anomalies, and to get a better optical depth estimate. $\tau_{m,outer}$ was found to be systematically higher than $\tau_{m,inner}$ based on the relative intensities. The variations in the intensity distribution get magnified in the $\tau_{m}$ maps and using the mean of the two inner or the two outer satellites does not nullify this effect. In some regions, the intensity ratios are so skewed, that we get non-physical values of $\tau_{m}$, particularly for $\tau_{m,inner}$. Due to these effects, over the entire map area, the $\tau_{m,inner}$ and $\tau_{m,outer}$ only have a modest positive correlation of about 0.31.

Thus, the assumption involving two satellite components is only partially successful in deriving a consistent optical depth. Consequently, we add another modification to the LTE approximation by using all the four satellite intensities. To obtain a representative $\tau_{(1,1,m)}$ for NH3 (1,1), we (i) scale up the outer satellite intensities by a factor of 5/4 (ratio of theoretical intensities), (ii) take the mean of all four intensities and (iii) solve equation 2 using a value of 5/18 for $\alpha$. The resulting optical depth map is shown in Figure 11. All the regions are fairly optically thick (median value 2.1), with no particularly identifiable distribution pattern. We use this value of $\tau_{(1,1,m)}$ for the next analysis steps.

5.2 Temperatures

The level populations of the two quantum states ($n_{+}$ and $n_{-}$) involved in the NH3 (1,1) inversion transition are related by an excitation temperature $T_{ex}$ in the equation

[TABLE]

where $g_{+}$ and $g_{-}$ are the statistical weights of the two levels. The statistical weight ratio is 1 for the inversion transitions of NH3. The transition frequency is $\nu$. The Boltzmann constant and the Planck constant are denoted by $k$ and $h$, respectively. The sum of $n_{+}$ and $n_{-}$ gives the total level population of the corresponding J=K quantum state.

If we assume that a single value of $T_{ex}$ is applicable to all the hyperfine components of the inversion transition, we can use the optical depth $\tau_{(1,1,m)}$ and the brightness temperature $T_{b(1,1,m)}$ of the main component of NH3 (1,1) to determine the NH3 (1,1) excitation temperatures $T_{ex(1,1)}$ by solving for the equation

[TABLE]

where we have assumed a beam filling factor of 1. The (1,1) inversion transition frequency is $\nu$ and $\mathcal{J}_{\nu}(T)$ is the Rayleigh-Jeans equivalent temperature at a frequency $\nu$ of a black body at temperature $T$. It is given by $\mathcal{J}_{\nu}(T)$ = $\frac{h\nu}{k}/[\exp(\frac{h\nu}{kT})-1]$. The background radiation temperature $T_{bg}$ is 2.73 K. The resulting excitation temperature map is shown in Figure 12. The median $T_{ex(1,1)}$ value is 9.4 K and it increases to a maximum of 20 K in the IRAS 2 region.

Similar to $T_{ex}$, a rotational temperature $T_{rot}$ can be used to define the ratio of the NH3 level populations for (J,K) = (1,1) and (2,2) (column densities $N_{(1,1)}$ and $N_{(2,2)}$, respectively):

[TABLE]

The factor 5/3 comes from the J-degeneracies and 41.0 K corresponds to the energy difference between the (1,1) and the (2,2) levels. From this equation, it is possible to derive $T_{rot}$ in terms of the NH3 (1,1) optical depth $\tau_{(1,1,m)}$, the ratio of the NH3 (1,1) and NH3 (2,2) main component intensities ($T_{R(2,2,m)}/T_{R(1,1,m)}$) and the ratio of the NH3 (1,1) and NH3 (2,2) velocity dispersions ($\sigma_{v(2,2)}/\sigma_{v(1,1)}$) as follows (Mangum et al., 1992):

[TABLE]

The rotational temperature can be used to solve for the kinetic temperature $T_{K}$ in the following equation (Swift et al., 2005) based on statistical equilibrium and detailed balance

[TABLE]

The rotational and kinetic temperature maps are shown in Figures 13 and 14, respectively. Both maps have very similar distributions. The median values of $T_{rot}$ and $T_{K}$ are 13.0 K and 12.4 K, respectively. In central NGC 1333, there are many regions of high temperature near the three main Class 0/I YSO groups. In addition, the west part of the SVS 13 region has higher temperatures between 14 K and 17 K. At the southern-west tip of this region, the temperature reaches 20 K. The northern part of the NW filament near the HH 12 region also has higher temperatures of $16-18$ K compared to the other filaments.

Figure 14 also has the Herschel dust temperature map contours (Pezzuto et al., 2017) overlaid for comparison. Most of the filaments have dust temperatures between 11 K and 13 K, less than the temperatures ($\sim$ 14 K) in the less dense areas of the cloud. The dust temperature values in the filaments and quiescent regions are in fair agreement to the kinetic temperatures obtained for NH3, indicating that the dust is thermodynamically coupled to the dense gas (Forbrich et al., 2014). The dust temperatures, however, are greater than the NH3 temperatures near the YSOs, with values greater than 20 K near IRAS 4, IRAS 2 and SVS 13. The dust temperatures are also greater than 16 K near IRAS 7 and HH 12. In these regions, the effective dust temperatures are elevated because the line-of-sight dust emission is dominated by contributions from the warm regions near the protostars, unlike the NH3 which traces the entire column of dense gas (Schnee et al., 2006). We also note that in regions closer to the YSOs, there are greater temperature variations in the high resolution NH3 map than in the dust map.

5.3 Column Density

To derive the column density of NH3 we first express the column density of the upper level of the (1,1) inversion transition $N_{(1,1),+}$ in terms of the optical depth $\tau_{(1,1)}$ and the excitation temperature $T_{ex(1,1)}$ using the following equation

[TABLE]

where the dipole moment of the NH3 molecule $\mu$ has a value of 1.468 Debye. Using equation 3, the ratio of the column densities of the two inversion states in (1,1) can be written as

[TABLE]

where $N_{(1,1),-}$ is the column density of the lower energy level of NH3 (1,1). Hence, the total column density of the NH3 (1,1) is

[TABLE]

Assuming a single rotational temperature $T_{rot}$ to define all the (J,K) level populations, the total NH3 column density $N_{tot}$ can be estimated from $N_{(1,1)}$ and the rotational partition function $Q_{rot}$ using the following equation (Mangum et al., 1992):

[TABLE]

where $g_{J}$ and $g_{K}$ are the rotational degeneracies and $g_{I}$ is the spin degeneracy. For NH3 (1,1) their values are 3, 2 and 0.25, respectively. $E_{(1,1),+}$ is the energy of the upper level in the NH3 (1,1) inversion transition. The value of $E_{(1,1),+}/k$ is 24.35 K where $k$ is the Boltzmann constant. This derivation uses a partition function that includes both ortho and para species of NH3.

By combining equations 8, 10 and 11, we get

[TABLE]

Since the optical depths for the different components of NH3 (1,1) are different, we express the total optical depth in terms of the main component optical depth of NH3 (1,1) using a relative intensity factor $R_{i}$ = 0.5. By using equation 4, we can express this optical depth in terms of the integrated intensity of the main component of NH3 (1,1) ($\int T_{R(1,1,m)}dv$). We use a correction factor of $\tau/(1-e^{-\tau})$ for the optically thick case (Goldsmith & Langer, 1999) which is applicable for most of the regions being analyzed. The resulting column density equation is thus:

[TABLE]

The expression can be simplified to:

[TABLE]

where $\int T_{R(1,1,m)}dv$ is in K km/s. The optical depths, excitation temperatures and rotational temperatures were calculated in the previous sections. To calculate the partition function $Q_{rot}$, we use the formula:

[TABLE]

The K-degeneracy $g_{K}$ is 1 for K = 0 and 2 for K $\neq$ 0. The nuclear spin degeneracy value $g_{I}$ is 0.25 for K = 3n (ortho species) and 0.5 for K $\neq$ 3n (para species) , where ‘n’ is a non-negative integer. The energy levels $E_{JK}$ up to J=5 are available in Table 9 of Mangum & Shirley (2015) and are sufficient for our rotational energy temperature ranges.

The resulting column density map of NH3 is shown in Figure 15. The column density values vary over two orders of magnitude, with a median value of 9.9 $\times$ 1014 cm*-2*. In the eastern part of the SVS 13 region, the NH3 column density is as high as 5 $\times$ 1015 cm*-2*. This region coincides with the minimum dust temperatures ($\sim$ 10.5 K) as well as gas kinetic temperatures ($\sim$ 10 K) in NGC 1333. The filaments have lesser column density than the cloud center, with the NW filament having the least column densities (mean value of 7 $\times$ 1014 cm*-2*). We overplot the dust-based H2 column density map contours (from Herschel) for comparison. The dust map traces the same structures as the NH3 map (at the 36*′′* Herschel beam scale), including the narrow filament in the SK 1 region. The typical values are 3 $\times$ $10^{22}$ cm*-2* in the filaments and goes up to 9 $\times$ 1022 cm*-2* at the cloud cores. The only exception is the region close to the IRAS 4 sources, where the dust-based column density is much higher relative to the NH3 column density in that region.

The expressions for the optical depth, the temperatures and the column density depend on the many terms and their error estimation is more involved. The previously determined uncertainties of the terms in the respective formulae are used as standard deviations of zero-centered normal distributions and added to the nominal pixel values of the terms to obtain alternative values. These alternative values are used to obtain alternative estimates of the expression being calculated at each pixel. The difference between the nominal expression and this alternative expression gives an error map. The standard deviation of the distribution of values in these error maps gives the respective uncertainties.

6 Discussion

One important question about NGC 1333 that arises is how well the different NH3 maps show the true characteristics of its dense gas and dust. In the previous section, we saw good spatial correlation with the Herschel temperature and column density maps, but that comparison was limited by the 10 times larger beam size of Herschel-derived maps compared to our VLA and GBT combined maps. In this section, the regions are compared with the CLASSy-I N2H*+* maps covering the entire NGC 1333 region (Storm et. al. in prep.) and with the JCMT 850 $\upmu$m dust emission map (Chen et al., 2016) from the JCMT Gould Belt Survey.

6.1 Cloud Morphology

The CLASSy-I N2H*+* (J=1-0) maps at 93.173 GHz have a synthesized beam size of 9.06*′′* $\times$ 7.58*′′, a channel width of 195.2 kHz (0.63 km/s) and a sensitivity of 0.08 Jy/beam (Storm et al., 2014). Figure 16 (left) shows the N2H+* (J=1-0) integrated intensity map (color), with the NH3 (1,1) integrated intensity contours overlaid on it. Both integrated intensity maps include all the respective hyperfine components. The two maps show very good agreement on all scales. The N2H*+* maps have slightly better SNR which results in the emission having a greater extent, especially near the filaments. The only region of noticeable difference between the N2H*+* and NH3 (1,1) maps is at the north-eastern edge of the NW filament, which has relatively greater emission in N2H*+. This location is the same place where higher temperatures were calculated from the NH3* (1,1) and NH3 (2,2) maps (see Section 5.2).

The JCMT Gould Belt survey SCUBA observations included the NGC 1333 region (COMPLETE team, 2011), with a map beam size of 15*′′* at 850 $\upmu$m. A comparison with the NH3 integrated intensity map (Figure 16, right) also shows good correlation except in the regions associated with embedded YSOs, where the dust emission strongly increases due to increasing column density and temperature on small scales. On close comparison, we suspect that the JCMT map is offset by about half the beam size in the south-west direction. This offset is particularly evident in the NE filament, and the right arm of the narrow ‘u’ shaped filament in the SK 1 region. One clear fact from the JCMT comparison is that the integrated intensity maps of NH3 and N2H*+* are poor tracers of individual YSOs. This conclusion is also clear from the NH3 column density map in Figure 15.

The high values in the NH3 temperature map have good correspondence with the outflows of the region (Plunkett et al., 2013) observed using CO maps and $H\alpha$ maps. The Herbig Haro objects HH 7-11 are associated with the high temperatures in the northern part of the SVS 13 region. The high temperature region in the HH 12 region and the connected NE filament is identified in the literature as a bow shock structure (Knee & Sandell, 2000). The western part of the SVS 13 region and the parts of the IRAS 2 region are heated by the outflow from IRAS 2A.

6.2 Cloud Kinematics

To analyze the kinematics of N2H*+* relative to those of NH3, spectral fitting was carried out for all the seven hyperfine components of N2H*+, with variable amplitudes for the three F1 quantum number based groups of lines. The resulting line-of-sight velocity map and its difference with the convolved and regridded NH3* velocities are shown in Figure 17. The two maps are very well correlated in most regions to within the error in the N2H*+* spectral fitting (0.18 km/s). The results indicate that NH3 is tracing the same material as N2H*+* and the dust, and that most features in the NH3 map can be expected in high resolution maps of other similar dense gas tracers.

Importantly, the velocity difference between an ionic species (N2H*+) and a neutral species (NH3*) can be used to infer ion-neutral drift, or ambipolar diffusion (Mestel & Spitzer, 1956). The general agreement in the velocity maps presented in Figure 17 therefore indicates that there is good ion-neutral coupling over most of the NGC 1333 region (Flower, 2000). This behavior is consistent with theoretical predictions that ambipolar diffusion is less efficient in quiescent environments (see, e.g., Spitzer, 1968; Mouschovias, 1979, for the typical timescale of ambipolar drift).

Interestingly, the IRAS 4 and the SVS 13 regions show relatively large differences between the two velocities. These differences could be related to multiple velocity components as discussed in the following paragraph, but could not be ascertained using our data. Alternatively, the differences could also be related to large ion-neutral drift ($\gtrsim 0.2$ km/s), which could be explained by shock accelerations, as proposed by Li & Nakamura (2004) and examined and confirmed in various numerical simulations (e.g., Nakamura & Li, 2005; Klessen et al., 2005; Vázquez-Semadeni et al., 2011; Chen & Ostriker, 2012, 2014). Considering these two regions potentially represent the shock front of an expanding bubble that is responsible for the formation of the SE filament (see discussions below in Section 6.3.1), our results could be the first observational validation of shock-accelerated ambipolar diffusion.

We note that there are two small regions with $\gtrsim 1$ km/s velocity difference between N2H*+* and NH3 (the dark green area on the eastern side of SVS 13 and the dark red area on the western side of IRAS 4). In fact, on investigating our data with the F1=1-0 isolated hyperfine component of N2H*+* (J=1-0), we found that multiple overlapping velocity-coherent components are commonly seen in those locations. In most cases, these multiple components have similar relative intensities for both molecular species, and thus the spectral fitting, which does a single velocity fit, gives comparable line-of-sight velocities. In the two small regions with $\gtrsim 1$ km/s velocity differences, however, the relative intensities of the two velocity-coherent components in the two species are reversed, resulting in one getting selected by N2H*+* and another by NH3 and hence the significantly larger velocity differences.

Since a lot of these regions have multiple velocity components, it is challenging to get a good estimate of the velocities of each of these components using complex molecular tracers like N2H*+* and NH3, given their many hyperfine components. Attempts to fit these spectra with two velocity components produced multiple solutions, that are very different from each other but have almost equal fitting errors. This issue is further compounded by the variation of the relative intensities of the hyperfine components from region to region. To get a good representation of all the multiple velocity components, some of these regions (especially near the Class 0/I YSOs) need to be observed with dense gas tracers like H13CO*+*. It is optically thin in most regions and has no hyperfine splitting, thereby allowing easier identification of the different velocity components. This tracer would also allow better estimation of the velocity dispersion of each component. A single component fit sometimes overestimates the dispersion in the dynamic environment near protostars having multiple velocity components.

6.3 Filaments in NGC 1333 and their relation to Star Formation

We identified three filaments in the NGC 1333 region (the SE, NE and NW filaments) and some narrower filaments in the SK 1 region. The FWHM values of the filaments range from 0.015 pc in the SK 1 region to 0.05 pc in the NW filament. The corresponding deconvolved widths from the JCMT map after taking into account the larger beam size, show that the filament widths in dust are greater than those in dense gas by about 50%. These results match those found for the CLASSy-II filaments (Dhabal et al., 2018).

6.3.1 Filament formation by colliding turbulent cells

The NH3 (1,1) map of the SE filament indicate a presence of two parallel sub-structures, as seen also in H13CO*+* and N2H*+* observations. The analysis of the H13CO*+* map of the same region showed that one of the two sub-structures has a significant velocity gradient across it (Dhabal et al., 2018). With our large-area, high-resolution NH3 (1,1) maps, we can put this in perspective of the kinematics of the entire cloud. We note that in the southern part of NGC 1333 there is a large region having line-of-sight velocities around $6-7$ km/s, which are at least 1 km/s less than those in adjacent regions. We see a large velocity gradient all along the southern edges of the SE filament, the IRAS 4 region, the SVS 13 region and the IRAS 2 region. The gradient is as much as 2 km/s across 0.09 pc in the SVS 13 region. It continues into the IRAS 4 region and the SE filament, where the gradients are comparable ($1-1.5$ km/s across 0.04 pc). The IRAS 2 region also shows a ridge of high velocity resulting in a sharp gradient in the SE-NW direction. We propose that this behavior could indicate a large-scale ($\sim$ 0.5 pc) turbulent cell moving towards the cloud center from the south (in the projected sky plane). The collision between this cell with the cloud could result in the formation of a layer of compressed gas, from which dense star-forming clusters (IRAS 4 and IRAS 2) emerge, as well as an elongated structure (the SVS 13$-$IRAS 4$-$SE filament system) at the shock front. This scenario is illustrated in Figure 18.

The importance of supersonic turbulence in structure formation is well established by numerical simulations with or without magnetic fields (Padoan et al., 1999; Banerjee et al., 2006; Pudritz & Kevlahan, 2013). Generally driven by stellar winds and supernovae (Mac Low & Klessen, 2004), these large-scale supersonic flows generate locally planar shock layers of compressed gas and help in transferring energy to smaller scales. In such shock layers, dense structures can form, either seeded by turbulent velocity perturbations (Chen & Ostriker, 2014, 2015) or at the intersection of shock fronts (Pudritz & Kevlahan, 2013). With increasing density, gravity starts becoming dominant, amplifying the initial anisotropies and resulting in the formation of protostellar systems.

In the SE filament region, a velocity gradient was identified across the filament with H13CO*+* (Dhabal et al., 2018), which is also seen in N2H*+* (Figure 17) and NH3 (Figure 8). This feature is consistent with the filament formation model proposed by Chen & Ostriker (2015) (see also Chen et al. in prep.). In this model, the velocity gradient is a projection effect of gravity-induced accretion within the shocked layer created by colliding turbulent cells. We note, however, that in this model the ratio between gas kinetic energy and gravitational energy, $C_{v}\equiv\Delta{v_{h}}^{2}/(GM(r)/L)$, is typically $\sim~{}0.1$ from numerical simulations (Chen et al. in prep.). Here, $\Delta v_{h}$ is half of the velocity difference across the filament out to a transverse distance r from the filament spine, and $M(r)/L$ is the mass per unit length of the filament measured at the same distance from the spine as the $\Delta v_{h}$ measurement. For the SE filament, we estimated the mass per unit length based on the Herschel column density map obtained from the Gould Belt Archive (Sadavoy et al., 2014), and derived $M/L\approx 20$ M*☉*/pc (Dhabal et al., 2018). Combining with the measurement $\Delta v_{h}\approx 0.35$ km/s, we get $C_{v}~{}\sim~{}1.6$, much larger than those expected from gravity-induced velocity differences. This disparity suggests that the SE filament is at the front-end of the expanding bubble, and the velocity gradient observed across the filament is due to the shock, because gravity from the filament itself is not sufficient to generate that velocity difference.

There are additional observational evidences of the turbulent cell collision scenario. Firstly, previous studies on polarimetry around this area using optical/near-IR extinction (Alves et al., 2011) and dust thermal emission (JCMT BISTRO; Doi et al. in prep.) have both found that the magnetic field within the SE filament is likely parallel to the filament, which is a feature of shock-compressed gas because shocks only amplify the magnetic field component that is parallel to the shock front (see, e.g., Shu, 1992). Secondly, the locations along the high velocity gradient ridge also have relatively high velocity dispersions ($>$ 0.5 km/s) in the NH3 (1,1) map that are typical for shocked layers. Thirdly, we see many protostars along the boundary of this proposed compressed gas layer. This region has the highest number of Class 0/I YSOs in NGC 1333. They are divided into the multiple systems IRAS 2, IRAS 4 and SVS 13. The star formation in this region could be triggered by gravitational instabilities in the dense gas layer. Similar interpretations have been made for cloud-cloud collisions based on observations (Duarte-Cabral et al., 2011; Nakamura et al., 2014; Dewangan & Ojha, 2017), and such collisions are often reported in simulations (Vázquez-Semadeni et al., 2003; Inoue & Fukui, 2013; Takahira et al., 2014). Finally, water masers are identified near all three afore-mentioned protostellar systems along the velocity gradient ridge (see Figure 5). The relative kinetic energies of shocked and unshocked gas are known to power masers; hence their presence is also an indication of shocks propagating in dense regions (Elitzur et al., 1989; Hollenbach et al., 2013).

We detect multiple parallel sub-structures in the SE and NE filaments, and all along the proposed compressed gas layer. Similar observations of parallel sub-structures in the Serpens filaments and in observations by other authors (Beuther et al., 2015; Hacar et al., 2013) indicate that such behavior is an important feature related to filament formation. It may be argued that this phenomenon is a density selection effect, wherein the filament spine has a density much greater than the critical density of the used tracers (Chen et al. in prep.). Thus, for these tracers, the gas in the filament spine emits less than the surrounding, less-dense part of the filament. Hence, a single physical filament could appear to be composed of multiple parallel sub-structures. This hypothesis may be tested using optically thin tracers of very dense gas such as H13CN emission.

Considering the alternative that the multiple parallel structures are physically separate, a common formation mechanism like the ‘fray and gather’ model (Smith et al., 2016) has been proposed, based on hydrodynamic simulations of turbulent clouds. In this model, the sub-filaments are formed first and are then gathered together by large scale motions of the cloud. Alternatively, it could be explained as multiple filaments forming in the dense gas layer created by colliding turbulent cells, which on projection in the sky plane can seem parallel to each other. Another explanation involves formation of a wide filament that fragments into multiple sub-filaments (Tafalla & Hacar, 2015).

The previously available GBT maps of NGC 1333 (Friesen et al., 2017), by themselves, were unable to resolve the rich structure within the filaments, all of which have narrower widths than the GBT beam size at 22 GHz. The red-shifted and blue-shifted parts of the cloud from which we propose the turbulent cell theory is also discernible in the GBT-only velocity map. The narrow shock zone with the high velocity gradient and high velocity dispersion, however, is resolved only in the combined VLA and GBT maps.

6.3.2 Effect of Outflows on Filaments

Outflows are known to clear the gas surrounding forming stars, and eventually lead to disruption of the larger cloud (Arce, 2002). They could also gravitationally unbind the gas in the dense core and stop the stellar mass accretion phase (Tafalla & Myers, 1997). In the NGC 1333 region, there is more energy in outflows than is observed in bulk, turbulent motions. Thus despite having a tiny fraction of the total mass, outflows have a key role in shaping the region. It is estimated that at the current rate of energy ejection from the outflows in NGC 1333, the entire cloud will be dispersed in 10 million years (Curtis et al., 2010).

The NE filament and the NW filament are in locations known to be impacted by outflows. Knee & Sandell (2000) showed that the cavity between these filaments is filled with high-velocity gas from several outflows. They identified HH 12 as the leading bow shock of a Class 0 source in the SVS 13 multiple system. The region inbetween these filaments is clear of dust and gas, and many evolved YSOs and main sequence stars are present here. Past outflows from these sources could have played an important role in the formation of the NE and NW filament by clearing out the cloud. Quillen et al. (2005) interpreted the two filaments as the walls of ancient outflow cavities which are no longer actively driven.

The NE filament shows major velocity variations along its length ($\sim$ 0.7 km/s) and across its width ($\sim$ 0.5 km/s) with a ridge of maxima in the north and a ridge of minima in the south. In addition to the evolved sources on the west of the NE filament, the outflow from one of the IRAS 7 sources (Bally et al., 1996) can also contribute to the morphology of the NE filament from the eastern side. Based on the CO outflows identified by Plunkett et al. (2013), we find that one of their outflow candidates ‘C3’ is aligned perfectly along part of the eastern wall of the NE filament.

The narrow short filament extending towards the south-west from the IRAS 2 region is also found to be parallel to the southern lobes of the outflows from IRAS 2A and IRAS 2B (Plunkett et al., 2013). In addition, the filament is coincident with the overlap region of these two outflows lobes. Although these evidences suggest that the short filament may have been formed in the shock front produced by the swept up gas between these two parallel outflows, this filament has relatively low NH3 (1,1) velocity dispersion values ($\sim$ 0.16 km/s) in Figure 9. This is not expected for shocked regions. The narrow ‘u’ shaped filament in the SK 1 region has also been proposed to be a dust shell by Lefloch et al. (1998). However, as with the other filament in this region, this filament has low velocity dispersion values, which is not expected in bow shocks. The cavity to its north could still have been cleared out by outflows from one or more of the sources in the cloud core. Similarly, the IRAS 2 outflows may be responsible for clearing the material on the two sides of the narrow south-west filament.

6.3.3 Star Formation

We find that many of the Class 0/I YSOs are embedded in the cloud, but not necessarily confined to the filaments. Of the isolated protostars, there is one Class 0/I source at the southern tip of the NW filament, and another similar source SK 1 in the ‘u’ shaped filament. Most of the remaining Class 0/I sources are present in protostellar multiple systems like IRAS 2, IRAS 4 and IRAS 7. Of these, seven of the sources are present along the proposed high-velocity-gradient ridge identified in Section 6.3.1. The compression of the gas to very high densities can make that gas gravitationally supercritical and induce cores and protostars to form within them (Balfour et al., 2015; Smith et al., 2014). Belloche et al. (2006) came to similar conclusions about the origin of IRAS 4A and IRAS 4B based on the high mass infall rate measured in the IRAS 4A envelope, which is only possible in a model involving a strong external compression wave.

Volgenau et al. (2006) used high resolution maps of the IRAS 2 and IRAS 4 protostellar systems to show that the cores harboring these sources are not quiescent homologous structures. Even on the scales of the protostellar envelopes, turbulence is the major contributor to the observed line-widths. Many of the sources such as IRAS 4A, IRAS 4B and IRAS 2A are close binary protostars. (Of these, IRAS 2A is not resolved in Spitzer data; we mark it as a single protostar although it has been confirmed to be a binary using radio data.) Simulations have also shown that fast compression is often associated with core fragmentation (Hennebelle et al., 2003) which could explain the recent results reported by Sadavoy & Stahler (2017) that most stars form initially in wide binaries. In addition, based on the uncorrelated direction of the outflows from the binary system of IRAS 2A, Tobin et al. (2015) concluded that it is evidence of core/envelope fragmentation by turbulence and not disk fragmentation due to gravitational instability. These observations indicate that turbulence at different scales plays an important role in the star formation process, particularly in the southern half of NGC 1333.

There is a string of Flat YSOs in the less dense area north of the SVS 13 region. Along the NW filament, there is a sequence of five Class II YSOs. The relatively evolved YSOs are spread throughout the cloud including the lower density regions. Indeed, outflows from these more-evolved objects may have previously cleared out some regions of the cloud. The formation of the Class 0 YSO sources – VLA 42 in the HH 12 region and SK 1 have been proposed to be triggered by outflow bow shocks by Sandell & Knee (2001). We do not find sufficient evidence for such shocks in the case of SK 1, because of the low velocity dispersion values along the proposed shock layer. In the case of VLA 42, however, there is evidence of high temperatures ($>$ 16 K) and large velocity dispersions ($\sim$ 0.4 km/s) in its environs over a large area (compared to the VLA beam area). Bright $H\alpha$ emission is also detected in this region (Walawender et al., 2005), which further supports a shock-triggered star formation scenario for VLA 42. Thus, outflows could be instrumental in regulating star formation in NGC 1333.

7 Summary

In this paper, NH3 observations are used to study the morphology, kinematics and temperatures of the entire NGC 1333 region at an angular resolution of 4*′′*. We showed the importance of high-resolution, large-area maps that are sensitive to structure at all the spatial scales. To attain coverage, VLA interferometric data were combined with GBT single dish maps in the visibility domain. We have studied the effects of the various parameters that are important in the joint deconvolution of the visibility data.

For further analysis, we used data cubes corresponding to NH3 (1,1) and NH3 (2,2) inversion transitions, which have the brightest emission among all the observed transitions. The final data cubes have a channel width of 11.4 kHz, a synthesized beam of 3.93*′′* $\times$ 3.41*′′, and a sensitivity of 4 mJy/beam. We produced the integrated intensity maps for both NH3* (1,1) and NH3 (2,2). Line fitting was carried out taking into account all the hyperfine components of NH3 (1,1) to obtain the peak intensity maps of the main and satellite components, the line-of-sight velocity map and the dispersion map. For NH3 (2,2), we used a similar method on the main component to obtain its peak intensity map and the dispersion map. We also derived maps of the optical depth for the NH3 (1,1) main component, the rotational temperature, the kinetic temperature and the column density of NH3.

Using these maps, we studied the distribution and properties of the dense gas in NGC 1333. Our main conclusions are summarized below.

• •

In locations near the cloud center and in some filaments, multiple velocity components in the cloud are observed.

• •

In many locations, there are large deviations from theoretical expectations in the relative intensities of the NH3 (1,1) main and satellite groups of hyperfine components. The observed anomalies in the NGC 1333 region do not support the theory involving hyperfine selective trapping (Stutzki & Winnewisser, 1985). Although these anomalies affect the optical depth map, the effects on the temperature maps are negligible.

• •

There is very good correlation of the NH3 (1,1) maps with the morphology and kinematics of the corresponding N2H*+* (J=1-0) maps. In regions away from the protostellar sources, the dust maps from JCMT and Herschel also match well with the NH3 (1,1) integrated intensity map. They all trace the same material in these regions.

• •

The northern part of the region has many active outflows, which possibly cleared out a cavity between two filaments.

• •

In the southern part of the map, based on the continuous large velocity gradient pattern, the presence of Class 0/I YSOs along the region of the gradient and other kinematic signatures, we have argued that this is a region of compressed gas formed by the collision of large-scale turbulent cells.

Overall, the NH3 observations of the entire NGC 1333 region are very rich in information and offer many clues about star formation in molecular clouds. We have postulated theories to explain some of these observations. Further analysis of the data needs to be carried out to study the regions around each of the protostars, by investigating the relation between the various map features at all scales.

Appendix A Weighting of Single Dish and Interferometric Visibilities

The value of the weighting parameter used prior to visibility combination in Section 3 involves matching the weight densities of the single dish visibilities and the interferometric visibilities. Since the single dish weight density falls off rapidly with uv distance, the matching is not very exact. In this section, we investigate the sensitivity of the maps and the total flux with the relative weighting.

Assuming the GBT weighting parameter that is used for the data analysis to be 1, we use GBT weights in the range of 0.25 to 4 on a single channel containing signal to obtain the corresponding dirty maps and study the progression of the cleaning process in each of them. We used the same ‘cycleniter’ parameter as before (1000), and set a threshold of 3 times the RMS noise in the dirty maps. As shown in Figure 19, in the nominal case the total flux remains close to the GBT flux (scaled for the beam area) throughout the cleaning cycles. For higher (lower) than nominal GBT weights, the total flux in the dirty map is very high (low), but with progressive iterations the final cleaned flux reaches closer to that of the nominal weight case. Eventually, they end up within a factor of 2 of the total flux in the nominal case for a weight factor of 4. Most of the contribution to this difference comes from the extent of the emission, which increases with increasing GBT weights. The differences in regions of peak emission are lesser ($<$ 10% for a weight factor of 4). The model flux and residual flux also end up at different ratios with respect to the total flux in the different cases. The beam areas only vary by about 2% for a factor of 4 in the weights. We compare the outcomes of the final cleaned channel maps for GBT weights 0.25 and 4 in Figure 20.