IN-SYNC. II. VIRIAL STARS FROM SUBVIRIAL CORES—THE VELOCITY DISPERSION OF EMBEDDED PRE-MAIN-SEQUENCE STARS IN NGC 1333

Spectra were obtained using the APOGEE multiobject spectrograph (Wilson et al. 2012). The details of the data reduction and spectral fitting are presented in Cottaar et al. (2014). In brief, model spectra from BT-Settl (Allard et al. 2012) spanning a range of properties are convolved to the resolution of the reduced IN-SYNC spectra. The main parameters that are allowed to vary in fitting the models to the observed spectra are the effective temperature (T_eff), the veiling (r_H), the rotational velocity (vsin i), the surface gravity (log (g)), and the radial velocity (v_r) of the star. Initial parameter uncertainties are estimated from the Markov Chain Monte Carlo fitting process, although these uncertainties are then inflated to match the actual epoch-to-epoch variability seen in the parameters.

IN-SYNC targets in NGC 1333 were chosen from the Cores to Disks (c2d) Spitzer survey of Perseus (Jørgensen et al. 2006; Rebull et al. 2007), supplemented with focused surveys of NGC 1333 (Getman et al. 2002; Gutermuth et al. 2008; Winston et al. 2009, 2010) and candidate members selected based on their mid-IR variability from a preliminary analysis of the light curves obtained in NGC 1333 by the Spitzer YSOVAR program (Rebull et al. 2014).

Stars from all input surveys were given equal weight when assigning targets, and the primary selection criterion was H magnitude. Since the APOGEE spectrograph has a fiber collision limit of 71'', we used three distinct plates in order to achieve nearly complete coverage of bright stars in the center of the cluster. We prioritized assigning fibers to stars with 7.5 mag < H < 13.5 mag; this was limited at the faint end by signal-to-noise ratio (S/N) considerations and at the bright end by the potential for flux bleeding. Sources with H < 12.5 mag were considered highest-priority targets, and then fainter sources were used to fill up a plate. A faint source that was selected for one plate was prioritized on subsequent plates so as to build up the S/N of faint stars.

This fiber assignment scheme ensured that only five NGC 1333 candidate members with H < 13.5 mag were not assigned a fiber on any of the NGC 1333 fiber plug plates. Otherwise, 141 likely NGC 1333 members with H < 13.5 mag were assigned fibers on at least one NGC 1333 plate; of these, 107 were assigned fibers on 2 or more plates, with 79, 61, and 58 members assigned at least 3, 4, or 5 fibers across both observing seasons. The observations are therefore close to complete for the bright stars in NGC 1333 but are significantly incomplete at fainter magnitudes and biased against faint stars in the densest portions of the cluster.

These magnitude limits do not correspond to simple limits on stellar mass, as the intrinsic H-band luminosity of a star in NGC 1333 may be significantly affected by extinction, either local (the envelope around a very young star) or global (from substructure within the dust and gas in the cluster as a whole). The former constraint limits the IN-SYNC sample to relatively older stars. Thus, the vast majority of IN-SYNC stars with a Gutermuth et al. (2008) classification from Spitzer are Class II rather than Class I stars simply because local extinction around Class I stars renders them very faint in H. For the typical 1 Myr old star in the sample, the H-band magnitude limit of 13.5 corresponds to roughly 0.1 M_☉ (see Section 2.4).

IN-SYNC observations of NGC 1333 include observations of stars in the "West-End" of Perseus well outside of NGC 1333. We use a minimal-spanning tree (MST; Gutermuth et al. 2009) on the positions of candidate protostars in the region to define a boundary for the cluster. This method shows a break in the cumulative distribution of span length at 015, or 0.65 pc. Cutting the MST at this span length produces a cluster boundary that corresponds well to what a by-eye identification of the cluster would provide and corresponds to all stars within the boundary of the box 518 < R.A. < 525 and 3115 < decl. <316. Table 1 presents the full list of candidate members in NGC 1333 used in this survey, not all of which were observed.

To compare with our stellar velocities, we measure the velocity and velocity dispersion of the local gas with a combination of different tracers. These include the ¹²CO (1–0) and ¹³CO (1–0) transitions, which trace relatively low volume-density gas; we use the data from the COMPLETE Survey (Ridge et al. 2006). In addition, we use maps from higher critical density transitions in the central region of the cluster. These maps were obtained from the James Clerk Maxwell telescope (JCMT) archive and include ¹²CO, ¹³CO, and C¹⁸O (3–2) from Curtis & Richer (2011). The velocities of dense cores within NGC 1333 are drawn from the N₂H⁺(1–0) observations of continuum sources from Kirk et al. (2007).

We have used the publicly available (André et al. 2010) Herschel data for NGC 1333 to construct a map of the dust column density over NGC 1333. This map was created from fitting the Herschel 160–500 μm data with a single-temperature modified blackbody where the dust opacities at each wavelength are given by the opacities in Ossenkopf & Henning (1994) for dust grains with thin ice mantles, rather than assuming a simple power-law modification of the blackbody spectrum (i.e., taking a single value of β). Choosing instead to use β = 2, normalized at 230 GHz (a typical assumption; see Schnee et al. 2010), produces a dust mass that is 10% greater. For the Ossenkopf & Henning (1994) model the opacity at 500 μm  κ₅₀₀, is 0.05 cm² g⁻¹; we assume a gas-to-dust ratio 100:1.

To account for large-scale gradients present in the Herschel data, the zero point of this map was set by matching the column densities obtained around the edges of NGC 1333 with the COMPLETE (Ridge et al. 2006) extinction map based on Two Micron All Sky Survey (2MASS) photometry. The COMPLETE extinction map is lower resolution than our new Herschel column density map (25 versus 36''), missing structure evident in the Herschel column density map. Additionally, the extinction map is also significantly biased at the position of the cluster by the presence of many embedded red stars. The Herschel dust column density map, anchored by the reliable (i.e., noncluster) portions of the extinction map, therefore provides the best available tracer of the dust (and therefore gas) mass in NGC 1333. For comparison, the total mass of the cluster gas is 20% lower when estimated from the COMPLETE extinction map rather than the Herschel-derived column density map.

The column density map derived from Herschel is shown in Figure 1, along with all the stars in the IN-SYNC target catalog; the stars observed by IN-SYNC are a subset of these objects.

Zoom In Zoom Out Reset image size

Table 2 shows the best-fit stellar parameters for all APOGEE spectra obtained in NGC 1333. For this analysis we have applied quality criteria to these fits. Specifically, we exclude all spectra with an S/N < 20 and stars for which the reduced χ² of the best model fit was >10, as these spectra lead to unreliable parameter determinations.

We then trim stars with very low or high rotational velocity (v sin i < 5 km s⁻¹ or v sin i > 75 km s⁻¹). The former cut removes only one star, 2M03291184+3121557, which is near the edge of the cluster and has a velocity far from the cluster mean; this is likely a contaminating field star. The latter cut may remove genuine cluster members, typically very early type stars dominated by hydrogen lines. The large v sin i for these stars means that they have broad absorption lines; these broad lines mean that the precision of the radial velocity fit is low, and therefore the stars exhibit a much broader spread in radial velocity than stars with more secure fits (see the abrupt increase in velocity spread in Figure 2).

Zoom In Zoom Out Reset image size

These cuts on v sin i are similar to the cuts applied in Cottaar et al. (2015) except that we adopt an upper limit of v sin i = 75 km s⁻¹ rather than 150 km s⁻¹, as this corresponds to the observed v sin i threshold for dramatically increased radial velocity scatter in NGC 1333. The poor radial velocity fits for these hot early-type stars mean that we are unable to make any conclusive statements about the velocity dispersion at the high mass/effective temperature end of the distribution. Using the Dotter et al. (2008) pre-main-sequence isochrones and an age of 1 Myr (see Figure 3), this cut effectively corresponds to removing all stars with M ≳ 3.5 M_☉.

Zoom In Zoom Out Reset image size

We remove all stars from the sample with a radial velocity uncertainty greater than 1.1 km s⁻¹, as these stars contribute little information about the velocity dispersion. We also remove stars with strong radial velocity variability. We estimate this in the same manner as Cottaar et al. (2015), by calculating the χ² statistic for the model of constant radial velocity as

$$\begin{equation} \chi ^2 = \sum _i{\frac{(v_i - \mu)^2}{\sigma _i^2}}, \end{equation} \tag{ 1 }$$

where v_i and σ_i are the best-fit radial velocity and associated uncertainty in each epoch, respectively, and μ is the uncertainty-weighted mean radial velocity. If the probability of obtaining at least this large a χ² is less than 10⁻⁴, we flag the star as a radial velocity variable and exclude it from the following velocity analysis. A total of 6% of stars are flagged as having variable radial velocities in this fashion.

Finally, after applying these cuts, there are three stars that lie closer to the main-sequence isochrone than the 1–2 Myr isochrone where most of the target clusters (see Figure 3(b), where these three stars are shown with star symbols). These three stars (2M03290289+3116010, 2M03293476+3129081, and 2M03295048+3118305) are all candidate members, rather than confirmed young stars, they are radial velocity outliers, and two of them lie on the outskirts of the cluster (Figure 3(a)); we therefore conclude that these three stars are most likely contaminating field stars, rather than genuine members of the cluster, and we exclude them from further analysis. This leaves 70 stars for the analysis of the cluster's velocity dispersion.

Table 3 shows the uncertainty-weighted mean parameters for all stars observed by IN-SYNC, along with how many spectra were used in the determination and a flag to indicate stars that were identified as having a variable radial velocity.

For comparison with the gas velocities, we have converted all stellar heliocentric radial velocities into the velocity frame of the gas data. This is the kinematical local standard of rest (LSRK), defined as a solar motion of 20 km s⁻¹ in the direction of α(J2000), δ(J2000) = (18:03:50.29, +30:00:16.8). No other velocity correction, such as correction for the gravitational redshift (Pasquini et al. 2011) or convective blueshift (Shporer & Brown 2011), has been applied. These corrections, which account for the difference between the velocity of the photosphere and the velocity of the center of mass of the star, are typically on the order of a few hundred meters per second (although the convective blueshift is poorly constrained for young stars). From a comparison with literature results, Cottaar et al. (2014) estimates that the systematic uncertainty in the absolute zero point of the IN-SYNC radial velocity system is on the order of 0.5 km s⁻¹.

From the Dotter et al. (2008) isochrones and the age of the cluster (1–2 Myr) from Figure 3(b), we can assign a mass to each star based on the effective temperature. The lowest mass calculated for the Dotter et al. (2008) isochrones is 0.1 M_☉ (corresponding to T_eff = 3000 K), and we do not see many stars with temperatures far below this. The mass function (number of stars per unit mass) increases down to the 0.1 M_☉ limit, suggesting that we are reasonably sampling the masses down to 0.1 M_☉. We are therefore sampling masses between 3.5 M_☉ and down to close to the hydrogen-burning limit. We are not sensitive to the velocity distribution of brown dwarfs or early-type stars.

The determination of the intrinsic velocity dispersion of the young stellar population in NGC 1333 requires two main corrections. First, the radial velocities of the young stars have significant and nonuniform uncertainties. Second, the velocity distribution will be inflated by the presence of binaries, particularly close binaries. Multiple epochs of radial velocity data allow for the identification of some binaries, but not all, particularly since 40% of the stars in our NGC 1333 sample only have one observation meeting our signal-to-noise ratio (S/N) and goodness-of-fit criteria.

We take three distinct approaches to infer the intrinsic velocity dispersion: (1) use the velbin package (Cottaar & Hénault-Brunet 2014) to model the velocities at all epochs and the influence of binaries simultaneously, (2) use an outlier-resistant analytic estimate for the velocity dispersion in the case of nonuniform uncertainties, and (3) trim all stars with radial velocity uncertainty greater than 0.5 km s⁻¹ and use outlier-resistant estimates of the width of the distribution. Each of these approaches makes different assumptions, so it is important to check that they produce consistent estimates for the intrinsic velocity dispersion. For comparison, using just the sample standard deviation to estimate the width of the velocity distribution (i.e., ignoring errors and binary contamination) provides an estimate of 1.1 ± 0.1 km s⁻¹.

We first infer the intrinsic velocity dispersion of NGC 1333 using the velbin package introduced in Cottaar & Hénault-Brunet (2014). This package generates a large sample of binary stars, with the mass ratio and orbital properties drawn from literature values. Specifically, we use the lognormal period distribution from Raghavan et al. (2010), the nearly flat mass ratio distribution from Reggiani & Meyer (2013), and the flat eccentricity distribution from Duchêne & Kraus (2013). This distribution is then sampled and compared against the observed velocities and velocity errors from the real data set to infer the fraction of sources that are binaries and thus deduce the intrinsic velocity width of NGC 1333's stars after accounting for the influence of binaries. The underlying velocity distribution is assumed to be Gaussian, and so the posterior probability distribution function is inferred for the center of the velocity distribution, the intrinsic velocity dispersion, and the binary fraction.

Figure 4 shows the posterior probability distribution functions inferred for NGC 1333's velocity distribution, marginalized over the fraction of binaries (which is poorly constrained with a fairly flat posterior distribution between 20% and 70%) and over both the binary fraction and the central velocity. The latter posterior probability distribution for the intrinsic velocity dispersion, σ₀, is only slightly asymmetric, with a most likely value of 0.92 ± 0.12 km s⁻¹ and a 95% credible interval of [0.72, 1.13] km s⁻¹. The central velocity, v_cen, is 8.02 ± 0.31 km s⁻¹. For comparison, the centroid velocity of the ¹³CO (1–0) gas in this region ranges from 7 to 9 km s⁻¹ (Quillen et al. 2005).

Zoom In Zoom Out Reset image size

The correction for binarity is relatively small because we have data over a three-year baseline and have already removed 6% of the stars as radial velocity variable (and thus likely binaries; see Section 2.4). The remaining stars are known not to have a large radial velocity signature over the three years they were observed and are therefore likely to either (1) not be in short-period binaries or (2) not have an edge-on inclination that produces a large radial velocity signature. Since we marginalize over all binary fractions in Figure 4, this estimate provides a conservative estimate for the intrinsic velocity dispersion, σ₀. If we fix the binary fraction at extreme values, we get the following: for a binary fraction of 80%, σ₀ = 0.89 ± 0.09 km s⁻¹; for a binary fraction of 20%, σ₀ = 0.98 ± 0.10 km s⁻¹.

We also consider the effect of a period cutoff on the lognormal period distribution from Raghavan et al. (2010). This has a small influence since the inferred binary fraction increases/decreases in order to match the observational constraints. If the binary fraction is held fixed at 50%, then there is some sensitivity to adopting period cutoffs. We consider imposing a cutoff on the maximum period. A semimajor cutoff of 1 AU provides an estimate of σ₀ = 0.98 ± 0.09 km s⁻¹, a semimajor cutoff of 10 AU gives σ₀ = 0.89 ± 0.10 km s⁻¹, and a semimajor cutoff of 100 AU (or more) gives σ₀ = 0.92 ± 0.10 km s⁻¹. This behavior is because the binaries that can most increase our observed velocity dispersion are binaries with intermediate periods—short-period binaries produce a large radial velocity signature that is easily ruled out by our multi-epoch observations, while long-period binaries produce small radial velocity variations. All these variations are well within our uncertainty given for σ₀. Our estimate of σ₀ is therefore relatively robust against changes in the assumed binary population.

The second estimate for the velocity dispersion considers binaries as velocity outliers and infers the parameters of the velocity distribution using robust estimators. For this purpose we use the median and the interquartile range to infer the center and width of the distribution. A single velocity for each star is calculated as the weighted mean over all observed epochs, and a single velocity uncertainty is calculated as the median error across the epochs. This approach therefore uses less information than the velbin method.

To account for the nonuniform errors in the epoch-averaged radial velocity measurements, we adopt the outlier-resistant estimator given by Ivezić et al. (2014)

$$\begin{equation} \sigma _0 = \sqrt{\zeta ^2 \sigma _G^2 - e_{50}^2}, \end{equation} \tag{ 2 }$$

where σ_G is the unbiased estimator of σ for a Gaussian based on the interquartile range:

$$\begin{equation} \sigma _G = 0.741 (q_{75} - q_{25}) \end{equation} \tag{ 3 }$$

and

$$\begin{eqnarray} e_{50} &=& \mathrm{median}(e_i),\\ \zeta &=& \frac{\mathrm{median}(\sigma _i)}{\mathrm{mean}(\sigma _i)}, \end{eqnarray} \tag{ 4 }$$

and

$$\begin{equation} \sigma _i = \sqrt{\sigma _G^2+e_i^2-e_{50}^2}, \end{equation} \tag{ 6 }$$

where the errors on individual stellar radial velocities are denoted as e_i. Uncertainties on this estimator are calculated from bootstrapping (Efron 1979). For the defined sample of stars in NGC 1333, this estimate of σ₀ is 1.04 ± 0.18 km s⁻¹. This broader confidence interval is indicative of the relative instability of this estimator for small samples (Ivezić et al. 2014).

Finally, based on the previous estimates of the intrinsic velocity dispersion, we trim all stars with radial velocity uncertainty >0.5 km s⁻¹. This allows us to approximate the errors as roughly uniform (and small, compared to the intrinsic velocity dispersion) and simply calculate

$$\begin{equation} \sigma _0 = \sqrt{\sigma _{r}^2 - e_{50}^2}, \end{equation} \tag{ 7 }$$

where σ_r is some robust estimator of the dispersion. Using either the median absolute deviation (Muller 2000) or σ_G for σ_r produces comparable results with σ₀ = 1.12 ± 0.18 km s⁻¹. This estimate for the dispersion is actually slightly greater than the estimate just from the sample standard deviation of the untrimmed data, as many of the stars with radial velocity uncertainty between 0.5 and 1.1 km s⁻¹ lie near average velocity. The uncertainty on this estimate is again fairly large, because we have significantly reduced the amount of data used in this estimate. Note that Cottaar et al. (2015) also trim all stars with radial velocity uncertainty >0.5 km s⁻¹ in their analysis of IC 348.

The estimates of the intrinsic velocity dispersion of the stars in NGC 1333 from velbin, the outlier-resistant analytic estimate, and the velocity dispersion of the error-trimmed subset are therefore all consistent with one another, within their relatively large uncertainties. We proceed with the velbin estimate of σ₀ = 0.92 ± 0.12 km s⁻¹, as this estimate uses the most information.

Figure 5 shows the positions and velocities of the dense cores and IN-SYNC stars compared with the low-density gas tracer, ¹³CO (1–0), which shows the general cloud gas. The cores and stars have similar, although not identical, spatial distributions. The cores have a small velocity dispersion (the one-dimensional velocity dispersion is 0.51 ± 0.05 km s⁻¹) and are strongly correlated with the highest-intensity regions of the diffuse gas; this confirms what other studies have found—dense cores are not moving ballistically with respect to their surrounding diffuse gas (Walsh et al. 2004; Kirk et al. 2007, 2010).

Zoom In Zoom Out Reset image size

The radial velocity errors on the stars would tend to diminish the appearance of any real correlation between the stellar velocities and that of the diffuse gas. Nonetheless, the structure of stars with low radial velocity errors reveals some cases where the stars are not well correlated with the diffuse gas. Comparing the radial velocity of the IN-SYNC stars with the first moment (i.e., the intensity-weighted mean) of the emission profile from the cloud gas (both ¹³CO (1–0) and C¹⁸O (3–2)) along the line of sight toward each of the stars shows no correlation between the diffuse gas velocity and the stellar velocity. This lack of a correlation arises because the first moment of the cloud gas is essentially the same at all locations in the cluster, while the stars have a broad spread in radial velocities.

Figure 6 displays the difference between stellar radial velocity and the centroid (first-moment) velocity of the ¹³CO (1–0) and C¹⁸O (3–2) gas. In both cases, the mean offset is close to zero, and the standard deviation of the offsets is roughly 1 km s⁻¹, which turns out to be comparable with the line width of the ¹³CO (1–0) gas.

Zoom In Zoom Out Reset image size

This measurement of the velocity dispersion in NGC 1333 represents the earliest measurement of this quantity.¹⁵ NGC 1333 is still actively forming stars, with dense gas cores that are either starless or host to early (Class 0/I) protostars. These cores lie in roughly the same spatial area as the young stars measured by IN-SYNC but have a much tighter velocity dispersion. Taking only the cores from Kirk et al. (2007) that lie within our defined spatial boundary, the one-dimensional velocity dispersion is 0.51 ± 0.05 km s⁻¹.

Note that the determination of the dense core velocity dispersion is not corrected for the effects of binarity. One might imagine, for instance, that a binary system with separate accretion envelopes could produce a skewed line profile in N₂H⁺. In the case of two envelopes with high relative velocities, the N₂H⁺ line profile could separate into two velocity components, although this situation is indistinguishable from multiple widely separated cores within the beam along the line of sight.

We can directly compare the velocity dispersion of the population of dense N₂H⁺ cores and the population of stars observed with IN-SYNC. Because the two populations lie in similar positions within the cluster, this comparison is relating the velocity dispersion of stars just before they form (presumably stars inherit the systemic velocity of the dense cores out of which they form) and that about 1–2 Myr after their birth.

Figure 7 shows the histograms of velocities for dense cores and stars. The distributions have consistent central velocities within their uncertainties (v_cen(cores) = 7.59 ± 0.09 km s⁻¹, v_cen(stars) = 8.02 ± 0.31 km s⁻¹). Because the stellar distribution appears broader owing to the significant and nonuniform uncertainties on the radial velocity determinations, we plot the Gaussian distributions inferred in Section 3.1, as well as the simple fit to the velocity distribution. The dense core velocity dispersion is significantly less broad than the stellar velocity dispersion inferred with velbin. The dense cores have a similar velocity dispersion to the line width seen in the dense gas tracer (C¹⁸O (3–2)), but significantly smaller dispersion than the diffuse gas (C¹³O (1–0)).

Zoom In Zoom Out Reset image size

The full virial equation for a molecular cloud is

$$\begin{equation} \frac{1}{2}\ddot{I} = 2(\mathcal {T} - \mathcal {T_S}) + \mathcal {M} + \mathcal {W}, \end{equation} \tag{ 8 }$$

where $\ddot{I}$ denotes the acceleration of the expansion/contraction of the cloud, $\mathcal {T}$ is the kinetic energy of particles and gas within the cloud, $\mathcal {T_S}$ is the surface term (surface pressure), $\mathcal {M}$ is the energy in the magnetic field, and $\mathcal {W}$ is the gravitational potential energy. A simple virial analysis normally assumes that one can neglect the surface term and the magnetic field, and that a cloud's expansion or contraction is not accelerating ($\ddot{I} = 0$) so that

$$\begin{equation} - \mathcal {W} = 2\mathcal {T}. \end{equation} \tag{ 9 }$$

Solving this equation for a spherical distribution with a power-law density distribution, ρ(r)∝r^−k, gives (see Bertoldi & McKee 1992) a virial velocity dispersion, σ_vir, via

$$\begin{equation} \sigma _{\mathrm{vir}} = \sqrt{\frac{a M G}{5 R}}, \end{equation} \tag{ 10 }$$

where R is the radius, M is the mass of the region under consideration, and a is a geometric factor of the order of unity:

$$\begin{equation} a = \frac{1-k/3}{1-2k/5}. \end{equation} \tag{ 11 }$$

Evaluating the virial state NGC 1333 therefore requires knowledge of the mass of the cluster at a given radius. This mass includes both the gas mass and the stellar mass. The projected gas distribution (Figure 1) shows significant substructure and nonazimuthal symmetry, and so estimates of the virial velocity assuming a smooth distribution will necessarily be only a rough approximation.

We compare two estimates of the gas mass within NGC 1333. The first is the mass estimate derived from the Herschel column density map calculated in Section 2.3. The second is using the COMPLETE ¹²CO and ¹³CO maps to derive the excitation temperature of the CO and then using the X-factors calculated for this region by Pineda et al. (2008) to convert the CO intensity to a total mass. The comparison between these two methods is shown in Figure 8. The conversion between CO and total mass depends on relations that were not calibrated in the central region of NGC 1333, since the extinction map used in the conversion was unreliable there (Pineda et al. 2008). For this reason, we use the mass estimated from the Herschel dust column density map, although the difference between the two profiles is not large.

Zoom In Zoom Out Reset image size

Figure 8 also shows the enclosed mass for ρ(r)∝r⁻¹ and ρ(r)∝r⁻² profiles. The observed mass profile is intermediate between these two simple profiles, matching the k = 1 profile at radii >0.6 pc, and close to the k = 2 profile at radii <0.6 pc. This change in profile occurs at the edge of the region of high column density (i.e., ∼15 mag of A_V) seen in the Herschel column density map (see Figure 1). The majority of the stellar population is contained within 800'' (0.97 pc at a distance of 250 pc), and we use this as the fiducial cluster radius.

IN-SYNC does not provide a complete census of all the stars in NGC 1333, so an estimate of the total stellar mass from these data requires significant extrapolation. Our input catalog contains 205 objects in the region considered. Low-mass regions such as this one have a typical mean stellar mass of 0.5 M_☉. In order to show how this mass is distributed in Figure 8, we simply assign 0.5 M_☉ to every star in the IN-SYNC catalog, giving a total mass of 102 M_☉. For comparison, Lada et al. (1996) estimate a total stellar cluster mass of 45 M_☉ over a similar region to that which is considered here. This is obviously a fairly rough approximation, but the stellar mass is much smaller than the gas mass and therefore has relatively little influence on the virial velocity of the cluster.

The virial velocity given by Equation (10) can be evaluated as a function of radius. Since the projected mass profile implies that the true density profile is between ρ(r)∝r⁻² and ρ(r)∝r⁻¹, we use a single value of ρ(r)∝r^−1.5 to provide a continuous value for the virial velocity. This leads to a value of a = 1.25. We show the result of this calculation in gray in Figure 9, along with an estimate of the uncertainty, which is dominated by the systematic uncertainty on the mass estimate, which is a combination of an uncertain zero point in our column density map, the uncertain dust emissivity, and our assumption of single temperature component along the line of sight. We adopt a factor of two uncertainty on the mass to account for these sources of uncertainty. The virial velocity takes on a roughly constant value beyond 0.6 pc (as would be expected for the ρ(r)∝r⁻¹ profile, which is observed at large radii) of σ_vir = 0.79 ± 0.20 km s⁻¹.

Zoom In Zoom Out Reset image size

Another estimate of the virial velocity comes from assuming that the diffuse gas is roughly in virial equilibrium, and therefore that the virial velocity dispersion of NGC 1333 can be estimated from the line width of the diffuse gas. As shown in Figure 7, the line profile for ¹³CO (1–0) integrated over the cluster region is non-Gaussian. We adopt the common approach of measuring the FWHM of the emission and then converting to an (effective) σ for a Gaussian distribution. The result of measurement as a function of radius is shown in green in Figure 9. There are several important caveats with this measurement: the presence of significant outflow energy and momentum in NGC 1333 (e.g., Lefloch et al. 1998; Knee & Sandell 2000; Plunkett et al. 2013) means that the CO in NGC 1333 may be supervirial; on the other hand, ¹³CO (1–0) may not be optically thin, so optical depth effects may increase the observed line width.

Finally, we can estimate the intrinsic velocity dispersion for both stars and dense cores, binned within a given radius. For this purpose we use velbin for the stars and the sample standard deviation for the cores (as we did in Section 3.1 and Section 3.3 for the full cluster). These measurements are displayed in Figure 9.

Figure 9 shows that the various estimates of the velocity dispersion approach constant values beyond about 0.4 pc. Inside this radius there is a suggestive dip in all measurements of the velocity dispersion, although for any individual tracer this is not a statistically significant decrease. At small radii, neglecting the surface term, $\mathcal {T_S}$, in Equation (8) is obviously incorrect since the surrounding mass of the cluster gas makes a significant contribution. Furthermore, it is difficult to interpret the measured radial velocity of stars at radii smaller than the orbital radii. Therefore, we take the velocity dispersions measured at large radii as indicative of the true state of the cluster.

The stellar velocity dispersion of 0.92 ± 0.12 km s⁻¹ is therefore consistent with our estimate of the cluster's virial velocity (0.79 ± 0.20 km s⁻¹) and the line width of the surrounding diffuse gas (1.1 ± 0.1 km s⁻¹); the dense cores have a velocity dispersion of 0.51 ± 0.05 km s⁻¹, less than the velocity dispersion of the diffuse gas, and consistent with subvirial.

Lada et al. (1996) first noted the bimodal spatial distribution of young stars in NGC 1333, which exhibits distinct northern and southern clusters in the stars observable in the near-infrared. The relatively unbiased survey of Gutermuth et al. (2008) confirmed the existence of these two subclusters, which can also be seen in Figure 3(a) as the clusters of stars around declination 3137and 3127, respectively. This clustering is less prominent in Figure 5 since the northern cluster contains fainter stars on average; fainter stars tend to have larger radial velocity uncertainties and are therefore suppressed in this figure.

Given the presence of subclusters, does it make sense to consider the velocity dispersion and virial velocity of the cluster as a whole? Using the Gutermuth et al. (2008) division of NGC 1333 into two clusters separated by a declination of 313, we use velbin to calculate the v_c and σ₀ for the northern and southern clusters. For the northern cluster, v_c = 7.9 ± 0.42 km s⁻¹ and σ₀ = 0.92 ± 0.18 km s⁻¹. For the southern cluster, v_c = 8.1 ± 0.40 km s⁻¹ and σ₀ = 0.94 ± 0.15 km s⁻¹. There is thus no detectable kinematic difference between the two subclusters, and considering them separately produces roughly the same result for the stellar velocity dispersion. We therefore proceed with considering the velocity dispersion of the full cluster.

Quillen et al. (2005) report a velocity gradient in ¹³CO (1–0) of 1 km s⁻¹ in the north-south direction across NGC 1333. However, this gradient is fit in a larger region than we consider here, extending 15' further south, and the magnitude of the gradient is strongly influenced by the substantial blueshift in emission in the southern portion of their field. The diffuse gas shown in Figure 5 certainly displays some large-scale structure but does not appear to exhibit a single consistent gradient. This is confirmed with higher-excitation CO lines; Bieging et al. (2014) present CO (2–1) and ¹³CO (2–1) maps of NGC 1333 and speculate that the patchy variations in centroid velocity seen in these tracers could be ascribed to a cloud–cloud collision.

We find no statistically significant velocity gradient in the stellar velocities across NGC 1333. For instance, considering just a north-south gradient, the best-fit velocity gradient across the region is −2.9 ± 2.2 km s⁻¹ deg⁻¹ (12.6 ± 9.6 km s⁻¹ pc⁻¹). Since the area under consideration is one-third of a degree, the uncertainty on the magnitude of this gradient is on the order of 0.7 km s⁻¹, leaving open the possibility of a gradient that is large enough to contribute to the observed spread in radial velocities seen in the region. Ultimately, we rely on the fact that we cover the same spatial region with our three tracers (diffuse gas, stars, and dense cores) and assume that any large-scale gradients or patchy structures affect the three tracers to a similar degree.

Proszkow et al. (2009) show that clusters that start with a subvirial velocity distribution become slightly supervirial as the cluster collapses, while an initially virial velocity distribution remains roughly constant. Nonspherical elongation of the cluster, particularly along the line of sight, makes it difficult to distinguish between the case of subvirial velocity dispersion and projection effects. Nonetheless, our results are broadly consistent with these simulations that start subvirial and then quickly become virial.

Offner et al. (2009) examined the velocity dispersion of young stars in a simulation of turbulent star formation. They find that star clusters forming in turbulent virialized clouds naturally begin with subvirial velocity dispersions; however, these subvirial velocity dispersions persist for one free-fall time. The dynamical time, t_dyn = R/σ_vir, of NGC 1333 is 1.1 Myr, so our result is in conflict with this prediction. Furthermore, the stars in the Offner et al. (2009) simulation retain a strong correlation with the centroid velocity of the gas in which they are embedded; our stars show no such correlation. This simulation seems to predict the behavior of the starless and protostellar cores traced by N₂H⁺, but not the behavior of the 1–2 Myr old stars that are observable in the near-infrared.

Both Kruijssen et al. (2012) and Girichidis et al. (2012) have recently studied the dynamical state of young stars in simulations. Kruijssen et al. (2012), in simulations with rather more mass than NGC 1333 (10⁴ M_☉ versus the 10³ M_☉ in NGC 1333), find that after one free-fall time most of the stars are in subclusters with relatively little gas left, and the stars are in viral equilibrium. Our study suggests that gas expulsion does not play a critical role in the virialization of young stars; the mass of NGC 1333 is still dominated by gas, and yet its young stars are virialized. Girichidis et al. (2012) study a slightly earlier stage, when only 20% of the mass is in protostars. In this study, only the central portions of the cluster are virialized; the external stars contribute to a subviral velocity distribution for the cluster as a whole. We see no evidence of this trend; the central regions and the outskirts of the cluster appear to have a similar viral state.

The difference between the 0.92 km s⁻¹ stellar velocity dispersion and the 0.5 km s⁻¹ dense core velocity dispersion can be explained in a number of ways. First, consider the additional terms in the virial equation that influence the gas and stars differently. After their formation, stars will cease to feel both the external pressure and the magnetic field. Since the dynamical time of NGC 1333 is 1.1 Myr, it is possible for the velocity dispersion of 1–2 Myr old stars to evolve subsequent to their formation.

We first consider the influence of external pressure on the region. We can calculate the external pressure using the formula from Lada et al. (2008):

$$\begin{equation} P_{\mathrm{S}} = 4.5 \times 10^3 \phi _G k A^2_V, \end{equation} \tag{ 12 }$$

where ϕ_G is a geometric factor of the order of 1 and A_V is the average column density in the region outside the cloud, which we can estimate as 2 mag of A_V for NGC 1333 from the COMPLETE extinction map. The energy due to surface pressure is simply

$$\begin{equation} \mathcal {T_S} = 4 \pi R^3 P_S. \end{equation} \tag{ 13 }$$

Our estimate for $\mathcal {T_S}$ is therefore 1.3 × 10⁴⁵ erg, which is less than 10% of the potential energy in the cluster ($|\mathcal {W}| = 3.3\times 10^{46}$ erg). The external pressure does not seem to be significant.

We cannot estimate the strength of the magnetic field in NGC 1333 directly, but if we assume that the difference between the dense core kinetic energy and the energy of the potential well in which they are embedded comes purely from the magnetic field, then we can calculate the strength of that field as

$$\begin{equation} \mathcal {M} = \frac{B^2 R^3}{6}, \end{equation} \tag{ 14 }$$

which gives B = 66 μG. This is a reasonable strength for the magnetic field in a region as dense as NGC 1333. Crutcher (2012) define an empirical relation between the strength of the magnetic field and the density of a region. For regions denser than the threshold particle density of n₀ = 300 cm⁻³,

$$\begin{equation} B = B_0 \left(\frac{n}{n_0}\right)^{0.65}, \end{equation} \tag{ 15 }$$

where B₀ = 10 μG. Assuming spherical symmetry, the particle density, n, of NGC 1333 is roughly 6 × 10³cm⁻³, and thus B = 70 μG.

Whether a magnetic field of this strength would actually provide support will depend at least partly on the field's morphology. In general, an ordered magnetic field will not provide support against collapse along magnetic field lines, but any realistic magnetic field configuration will involve turbulence translating support against collapse into all directions (see McKee & Ostriker 2007, and references therein). The exact details depend on how tangled the magnetic field is and how much energy density is in large-scale field components; further investigation with simulations is required.

The diffuse gas has a broader velocity dispersion, and it might be expected to be more strongly influenced by the magnetic field. However, the diffuse gas could well be out of virial equilibrium owing to the injection of energy from outflows from the existing young stars in NGC 1333.

It is therefore quite possible that NGC 1333 has a magnetic field strong enough to explain the dense cores' subvirial motion. In this picture, the dense cores' velocities are constrained by the magnetic field; upon their formation the stars are freed from the magnetic field and evolve to the velocity dispersion dictated by gravity.

The most difficult problem with this explanation is that the young stellar population in NGC 1333 does not seem to be in the relaxed, centrally condensed configuration expected for a virialized stellar population. As discussed in Section 4.1, the stellar population exhibits significant structure in the form of two subclusters. Furthermore, an analysis of the radial density profile of the young stars in NGC 1333 by Gutermuth et al. (2008) finds a flat central surface density profile, which is interpreted as evidence of very little dynamical evolution.

We consider a second explanation: that the cluster NGC 1333 might be in a state of global collapse. The IN-SYNC stars, which are older than the current dense cores, therefore formed when the cluster was more extended than it is today, presumably with an initially subvirial dispersion (like that possessed by the current dense cores). As the cluster globally collapsed, the potential energy of the stellar configuration was transformed into additional kinetic energy, rendering the stars dynamically hotter than the current population of dense cores. Cottaar et al. (2015) have recently proposed that IC 348 is in a state of global collapse based on an independent line of reasoning. This explanation would therefore lend support to the idea that young clusters are generally in a state of global collapse.

In this picture the subvirial velocity dispersion of dense cores can be explained if dense cores form at the convergent point of large-sale turbulent flows (Elmegreen 2007; Gong & Ostriker 2011) or if the cores form from a small number of velocity-coherent structures (perhaps filaments; Hacar & Tafalla 2011). The dense cores in NGC 1333 do exhibit coherence in position–velocity space (see top right panel of Figure 5), so this is a reasonable explanation. If the IN-SYNC stars come from similarly substructured initial conditions, where those structures are now undergoing global collapse, then perhaps the stellar velocity dispersion could be inflated within the substructures before the substructure is erased. In this model it is unclear why the stellar subclusters in NGC 1333 have consistent central velocities; the naïve expectation in the case of global collapse is that they would have different bulk motions. Additional modeling of this point is required.