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

NH₃ (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σ on the basis of the rms noise of each channel map. Figure 5 shows the integrated intensity map for the main NH₃ (1, 1) component alone. We have roughly divided the map into 11 regions according to their location in the map or the main embedded protostellar source(s) in the region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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. Toward 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 toward the southwest from the IRAS 2 region. The cloud core near IRAS 4 extends to a wide filament in the southeast, 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 NH₃ velocity structure of NGC 1333. Because 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 NH₃ (1, 1) components to vary. The magnetic dipole-based hyperfine splitting within each of the five groups also has spreads ranging between 11 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 five 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 σ_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:

$$\begin{eqnarray}&&S(v)=\displaystyle \sum _{n=1}^{5}{a}_{n}\displaystyle \sum _{i}{w}_{i}\exp \left(-\displaystyle \frac{v-{v}_{\mathrm{los}}-c{\rm{\Delta }}{f}_{i}/{f}_{0}}{2{{\sigma }_{v}}^{2}}\right),\end{eqnarray} \tag{ 1 }$$

where c is the speed of light, and Δf_i is the frequency offset of the hyperfine component relative to the line center frequency f₀ = 23.6944955 GHz. We solve for a₁, a₂, a₃, a₄, a₅, v_los, and σ_v. We mask the pixels in which none of the channels has more than 4σ detections (noise rms σ = 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–9. The intensity map distributions can be considered to be representative of the brightness temperature distributions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 NH₃ (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₁, a₂, a₄, and a₅ 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 southwest 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 do not 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 subfilaments 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 F₁ = 0–1 satellite component has greater mean intensities over most of the regions than the other satellites, while the F₁ = 2–1 and the F₁ = 1–0 satellites have relatively lesser mean intensities of all the satellites. The F₁ = 1–0 emission has lesser contrast between the strong emission regions and the surrounding regions. 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)–(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(F₁ = 1–0)/T_B(F₁ = 0–1) and T_B(F₁ = 2–1)/T_B(F₁ = 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(F₁ = 1–0)/T_B(F₁ = 0–1) and I = T_B(F₁ = 2–1)/T_B(F₁ = 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 NH₃ (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 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 that 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 that 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.

For the integrated intensity map of NH₃ (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, which would be useful for the analyses in the following section. The NH₃ (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 NH₃ (1, 1) integrated intensity emission. An outflow (from IRAS 2A) between the IRAS 2 region and the West region has NH₃ (2, 2) emission but is absent in the NH₃ (1, 1) map.

Zoom In Zoom Out Reset image size

The satellite components of NH₃ (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 NH₃ (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 NH₃ (1, 1) map to ensure fitting for the same velocity component as in NH₃ (1, 1). In single component regions, solving for the line-of-sight velocity gives the same result in both NH₃ 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 NH₃ (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 NH₃ (2, 2) integrated intensity have lesser NH₃ (2, 2) velocity dispersions and vice versa.

The optical depths can be estimated using the ratios of the main and satellite components of NH₃ (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 (α). The optical depth at any of the NH₃ inversion transitions' line centers τ_m can be calculated by solving for the equation (Mangum & Shirley 2015):

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{R}(m)}{{T}_{R}(s)}=\displaystyle \frac{1-{e}^{-{\tau }_{m}}}{1-{e}^{-\alpha {\tau }_{m}}},\end{eqnarray} \tag{ 2 }$$

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 NH₃ (1, 1), the theoretical value of α is 2/9. For the inner satellites, α is 5/18. For the NGC 1333 region, the LTE assumption is not accurate according to the anomalies presented in Section 4.1. Thus, we investigated a modified LTE assumption using two satellite components having the same theoretical α value. In this method, we calculated the optical depth for the main component separately using the mean outer satellite intensities ${\tau }_{m,\mathrm{outer}}$ and using the mean inner satellite intensities ${\tau }_{m,\mathrm{inner}}$. The mean was taken to reduce the effect of the intensity anomalies and to get a better optical depth estimate. ${\tau }_{m,\mathrm{outer}}$ was found to be systematically higher than ${\tau }_{m,\mathrm{inner}}$ on the basis of the relative intensities. The variations in the intensity distribution get magnified in the τ_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 nonphysical values of τ_m, particularly for ${\tau }_{m,\mathrm{inner}}$. Because of these effects, over the entire map area, the ${\tau }_{m,\mathrm{inner}}$ and ${\tau }_{m,\mathrm{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 }_{(\mathrm{1,1},m)}$ for NH₃ (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 α. 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 }_{(\mathrm{1,1},m)}$ for the next analysis steps.

Zoom In Zoom Out Reset image size

The level populations of the two quantum states (n₊ and n₋) involved in the NH₃ (1, 1) inversion transition are related by an excitation temperature T_ex in the equation

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{+}}{{n}_{-}}=\displaystyle \frac{{g}_{+}}{{g}_{-}}{e}^{-h\nu /{{kT}}_{\mathrm{ex}}},\end{eqnarray} \tag{ 3 }$$

where g₊ and g₋ are the statistical weights of the two levels. The statistical weight ratio is 1 for the inversion transitions of NH₃. The transition frequency is ν. 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 }_{(\mathrm{1,1},m)}$ and the brightness temperature ${T}_{{\rm{b}}(\mathrm{1,1},m)}$ of the main component of NH₃ (1, 1) to determine the NH₃ (1, 1) excitation temperatures ${T}_{\mathrm{ex}(\mathrm{1,1})}$ by solving for the equation

$$\begin{eqnarray}&&{T}_{{\rm{b}}(\mathrm{1,1},m)}=\left(\displaystyle \frac{1}{{{ \mathcal J }}_{\nu }({T}_{\mathrm{ex}})-{{ \mathcal J }}_{\nu }({T}_{\mathrm{bg}})}\right)[1-{e}^{-{\tau }_{(\mathrm{1,1},m)}}],\end{eqnarray} \tag{ 4 }$$

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

Zoom In Zoom Out Reset image size

Similar to T_ex, a rotational temperature T_rot can be used to define the ratio of the NH₃ level populations for (J, K) = (1, 1) and (2, 2) (column densities ${N}_{(\mathrm{1,1})}$ and ${N}_{(\mathrm{2,2})}$, respectively):

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{(\mathrm{2,2})}}{{N}_{(\mathrm{1,1})}}=\displaystyle \frac{5}{3}{e}^{-41.0/{T}_{\mathrm{rot}}}.\end{eqnarray} \tag{ 5 }$$

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 NH₃ (1, 1) optical depth ${\tau }_{(\mathrm{1,1},m)}$, the ratio of the NH₃ (1, 1) and NH₃ (2, 2) main component intensities (${T}_{R(\mathrm{2,2},m)}/{T}_{R(\mathrm{1,1},m)}$), and the ratio of the NH₃ (1, 1) and NH₃ (2, 2) velocity dispersions (${\sigma }_{v(\mathrm{2,2})}/{\sigma }_{v(\mathrm{1,1})}$), as follows (Mangum et al. 1992):

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{rot}} & = & -41.0\left\{\mathrm{ln}\left[-\displaystyle \frac{0.283{\sigma }_{v(\mathrm{2,2})}}{{\tau }_{(\mathrm{1,1},m)}{\sigma }_{v(\mathrm{1,1})}}\right.\right.\\ & & \times \,\left.\left.\mathrm{ln}\left(1-\displaystyle \frac{{T}_{R(\mathrm{2,2},m)}}{{T}_{R(\mathrm{1,1},m)}}(1-{e}^{-{\tau }_{(\mathrm{1,1},m)}})\right)\right]\right\}.\end{array}\end{eqnarray} \tag{ 6 }$$

The rotational temperature can be used to solve for the kinetic temperature T_K in the following equation (Swift et al. 2005) on the basis of statistical equilibrium and detailed balance

$$\begin{eqnarray}&&{T}_{\mathrm{rot}}={T}_{K}\left\{1+\displaystyle \frac{{T}_{K}}{41.5}\mathrm{ln}\left(1+0.6{e}^{-15.7/{T}_{K}}\right)\right\}.\end{eqnarray} \tag{ 7 }$$

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 and 17 K. At the southwest 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.

Zoom In Zoom Out Reset image size

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 and 13 K, less than the temperatures (∼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 NH₃, indicating that the dust is thermodynamically coupled to the dense gas (Forbrich et al. 2014). The dust temperatures, however, are greater than the NH₃ 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 NH₃ 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 NH₃ map than in the dust map.

Zoom In Zoom Out Reset image size

To derive the column density of NH₃, we first express the column density of the upper level of the (1, 1) inversion transition ${N}_{(\mathrm{1,1}),+}$ in terms of the optical depth ${\tau }_{(\mathrm{1,1})}$ and the excitation temperature ${T}_{\mathrm{ex}(\mathrm{1,1})}$ using the following equation

$$\begin{eqnarray}&&{N}_{(\mathrm{1,1}),+}=\displaystyle \frac{3{hJ}(J+1)}{8{\pi }^{3}{\mu }^{2}{K}^{2}}\displaystyle \frac{1}{({e}^{h\nu /{{kT}}_{\mathrm{ex}(\mathrm{1,1})}}-1)}\int {\tau }_{(\mathrm{1,1})}{dv},\end{eqnarray} \tag{ 8 }$$

where the dipole moment of the NH₃ molecule μ 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

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{(\mathrm{1,1}),+}}{{N}_{(\mathrm{1,1}),-}}={e}^{-h\nu /{{kT}}_{\mathrm{ex}(\mathrm{1,1})}},\end{eqnarray} \tag{ 9 }$$

where ${N}_{(\mathrm{1,1}),-}$ is the column density of the lower energy level of NH₃ (1, 1). Hence, the total column density of the NH₃ (1, 1) is

$$\begin{eqnarray}&&{N}_{(\mathrm{1,1})}={N}_{(\mathrm{1,1}),-}+{N}_{(\mathrm{1,1}),+}={N}_{(\mathrm{1,1}),+}({e}^{h\nu /{{kT}}_{\mathrm{ex}(\mathrm{1,1})}}+1).\end{eqnarray} \tag{ 10 }$$

Assuming a single rotational temperature T_rot to define all the (J, K) level populations, the total NH₃ column density N_tot can be estimated from ${N}_{(\mathrm{1,1})}$ and the rotational partition function Q_rot using the following equation (Mangum et al. 1992):

$$\begin{eqnarray}&&{N}_{\mathrm{tot}}=\displaystyle \frac{{N}_{(\mathrm{1,1})}{Q}_{\mathrm{rot}}}{{g}_{J}{g}_{K}{g}_{I}}{e}^{{E}_{(\mathrm{1,1}),+}/{{kT}}_{\mathrm{rot}}},\end{eqnarray} \tag{ 11 }$$

where g_J and g_K are the rotational degeneracies and g_I is the spin degeneracy. For NH₃ (1, 1), their values are 3, 2 and 0.25, respectively. ${E}_{(\mathrm{1,1}),+}$ is the energy of the upper level in the NH₃ (1, 1) inversion transition. The value of ${E}_{(\mathrm{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 NH₃.

By combining Equations (8), (10), and (11), we get

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{tot}} & = & \displaystyle \frac{3{hJ}(J+1)}{8{\pi }^{3}{\mu }^{2}{K}^{2}}\displaystyle \frac{{Q}_{\mathrm{rot}}}{{g}_{J}{g}_{K}{g}_{I}}{e}^{{E}_{(\mathrm{1,1}),+}/{{kT}}_{\mathrm{rot}}}\\ & & \times \,\displaystyle \frac{{e}^{h\nu /{{kT}}_{\mathrm{ex}(\mathrm{1,1})}}+1}{{e}^{h\nu /{{kT}}_{\mathrm{ex}(\mathrm{1,1})}}-1}\int {\tau }_{(\mathrm{1,1})}{dv}.\end{array}\end{eqnarray} \tag{ 12 }$$

Because the optical depths for the different components of NH₃ (1, 1) are different, we express the total optical depth in terms of the main component optical depth of NH₃ (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 NH₃ (1, 1) ($\int {T}_{R(\mathrm{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:

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{tot}} & = & \displaystyle \frac{3{hJ}(J+1)}{8{\pi }^{3}{\mu }^{2}{K}^{2}}\displaystyle \frac{{Q}_{\mathrm{rot}}}{{g}_{J}{g}_{K}{g}_{I}}{e}^{{E}_{(\mathrm{1,1}),+}/{{kT}}_{\mathrm{rot}}}\left[\displaystyle \frac{{e}^{h\nu /{{kT}}_{\mathrm{ex}(\mathrm{1,1})}}+1}{{e}^{h\nu /{{kT}}_{\mathrm{ex}(\mathrm{1,1})}}-1}\right]\displaystyle \frac{1}{{R}_{i}}\\ & & \times \,\left[\displaystyle \frac{\int {T}_{R(\mathrm{1,1},m)}{dv}}{{{ \mathcal J }}_{\nu }({T}_{\mathrm{ex}})-{{ \mathcal J }}_{\nu }({T}_{\mathrm{bg}})}\right]\displaystyle \frac{{\tau }_{(\mathrm{1,1},m)}}{1-{e}^{-{\tau }_{(\mathrm{1,1},m)}}}.\end{array}\end{eqnarray} \tag{ 13 }$$

The expression can be simplified to:

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{tot}} & = & 4.51\times {10}^{12}\ {Q}_{\mathrm{rot}}{e}^{24.35/{T}_{\mathrm{rot}}}\displaystyle \frac{1+{e}^{1.137/{T}_{\mathrm{ex}(\mathrm{1,1})}}}{1.517-{e}^{1.137/{T}_{\mathrm{ex}(\mathrm{1,1})}}}\\ & & \times \,\displaystyle \frac{{\tau }_{(\mathrm{1,1},m)}}{1-{e}^{-{\tau }_{(\mathrm{1,1},m)}}}\int {T}_{R(\mathrm{1,1},m)}{dv}\,{\mathrm{cm}}^{-2},\end{array}\end{eqnarray} \tag{ 14 }$$

where $\int {T}_{R(\mathrm{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:

$$\begin{eqnarray}&&{Q}_{\mathrm{rot}}=\displaystyle \sum _{J}\displaystyle \sum _{K}(2J+1){g}_{K}{g}_{I}{e}^{-{E}_{{JK}}/{{kT}}_{\mathrm{rot}}}.\end{eqnarray} \tag{ 15 }$$

The K-degeneracy g_K is 1 for K = 0 and 2 for $K\ne 0$. The nuclear spin degeneracy value g_I is 0.25 for K = 3n (ortho species) and 0.5 for $K\ne 3n$ (para species), where "n" is a nonnegative 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 NH₃ is shown in Figure 15. The column density values vary over two orders of magnitude, with a median value of 9.9 × 10¹⁴ cm⁻². In the eastern part of the SVS 13 region, the NH₃ column density is as high as 5 × 10¹⁵ cm⁻². This region coincides with the minimum dust temperatures (∼10.5 K) as well as gas kinetic temperatures (∼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 × 10¹⁴ cm⁻²). We overplot the dust-based H₂ column density map contours (from Herschel) for comparison. The dust map traces the same structures as the NH₃ map (at the 36'' Herschel beam scale), including the narrow filament in the SK 1 region. The typical values are 3 × 10²² cm⁻² in the filaments and go up to 9 × 10²² cm⁻² 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 NH₃ column density in that region.

Zoom In Zoom Out Reset image size

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 are 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.

The CLASSy-I N₂H⁺ (J = 1–0) maps at 93.173 GHz have a synthesized beam size of 906 × 758, 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 N₂H⁺ (J = 1–0) integrated intensity map (color), with the NH₃ (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 N₂H⁺ maps have slightly better S/N, which results in the emission having a greater extent, especially near the filaments. The only region of noticeable difference between the N₂H⁺ and NH₃ (1, 1) maps is at the northeastern edge of the NW filament, which has relatively greater emission in N₂H⁺. This location is the same place where higher temperatures were calculated from the NH₃ (1, 1) and NH₃ (2, 2) maps (see Section 5.2).

Zoom In Zoom Out Reset image size

The JCMT Gould Belt survey SCUBA observations included the NGC 1333 region (COMPLETE team 2011), with a map beam size of 15'' at 850 μm. A comparison with the NH₃ integrated intensity map (Figure 16, right) also shows good correlation except in the regions associated with embedded YSOs, where the dust emission strongly increases because of 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 southwest 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 NH₃ and N₂H⁺ are poor tracers of individual YSOs. This conclusion is also clear from the NH₃ column density map in Figure 15.

The high values in the NH₃ temperature map have good correspondence with the outflows of the region (Plunkett et al. 2013) observed using CO maps and Hα 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 are 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.

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

Zoom In Zoom Out Reset image size

Importantly, the velocity difference between an ionic species (N₂H⁺) and a neutral species (NH₃) 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 (≳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., Klessen et al. 2005; Nakamura & Li 2005; Vázquez-Semadeni et al. 2011; Chen & Ostriker 2012, 2014). Considering that 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 ≳1 km s⁻¹ velocity difference between N₂H⁺ and NH₃ (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 F₁ = 1–0 isolated hyperfine component of N₂H⁺ (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 ≳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 N₂H⁺ and another by NH₃ and hence the significantly larger velocity differences.

Because 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 N₂H⁺ and NH₃, 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 H¹³CO⁺. 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.

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 NH₃ (1, 1) map of the SE filament indicates a presence of two parallel substructures, as seen also in H¹³CO⁺ and N₂H⁺ observations. The analysis of the H¹³CO⁺ map of the same region showed that one of the two substructures has a significant velocity gradient across it (Dhabal et al. 2018). With our large-area, high-resolution NH₃ (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 (∼0.5 pc) turbulent cell moving toward 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.

Zoom In Zoom Out Reset image size

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 H¹³CO⁺ (Dhabal et al. 2018), which is also seen in N₂H⁺ (Figure 17) and NH₃ (Figure 8). This feature is consistent with the filament formation model proposed by Chen & Ostriker (2015; see also C.-Y. Chen et al. 2019, in preparation). 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 {\rm{\Delta }}{{v}_{h}}^{2}/({GM}(r)/L)$, is typically ∼0.1 from numerical simulations (C.-Y. Chen et al. 2019, in preparation). Here, Δ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 Δv_h measurement. For the SE filament, we estimated the mass per unit length on the basis of the Herschel column density map obtained from the Gould Belt Archive (Sadavoy et al. 2014), and derived M/L ≈ 20 ${M}_{\odot }$ pc⁻¹ (Dhabal et al. 2018). Combining with the measurement Δv_h ≈ 0.35 km s⁻¹, we get C_v ∼ 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 is additional observational evidence of the turbulent cell collision scenario. First, previous studies on polarimetry around this area using optical/near-IR extinction (Alves et al. 2011) and dust thermal emission (JCMT BISTRO; Y. Doi et al. 2019, in preparation) 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). Second, the locations along the high-velocity gradient ridge also have relatively high-velocity dispersions (>0.5 km s⁻¹) in the NH₃ (1, 1) map that are typical for shocked layers. Third, 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 on the basis of 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 aforementioned 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 substructures in the SE and NE filaments, and all along the proposed compressed gas layer. Similar observations of parallel substructures in the Serpens filaments and in observations by other authors (Hacar et al. 2013; Beuther et al. 2015) 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 (C.-Y. Chen et al. 2019, in preparation). 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 substructures. This hypothesis may be tested using optically thin tracers of very dense gas such as H¹³CN 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, on the basis of hydrodynamic simulations of turbulent clouds. In this model, the subfilaments 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 subfilaments (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 redshifted and blueshifted 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

6.3.3. Star Formation