Dynamics in Star-forming Cores (DiSCo): Project Overview and the First Look toward the B1 and NGC 1333 Regions in Perseus

The internal velocity structure within dense gaseous cores plays a crucial role in providing the initial conditions for star formation in molecular clouds. However, the kinematic properties of dense gas at core scales (~0.01 - 0.1 pc) has not been extensively characterized because of instrument limitations until the unique capabilities of GBT-Argus became available. The ongoing GBT-Argus Large Program, Dynamics in Star-forming Cores (DiSCo) thus aims to investigate the origin and distribution of angular momenta of star-forming cores. DiSCo will survey all starless cores and Class 0 protostellar cores in the Perseus molecular complex down to ~0.01 pc scales with<0.05 km/s velocity resolution using the dense gas tracer N$_2$H$^+$. Here, we present the first datasets from DiSCo toward the B1 and NGC 1333 regions in Perseus. Our results suggest that a dense core's internal velocity structure has little correlation with other core-scale properties, indicating these gas motions may be originated externally from cloud-scale turbulence. These first datasets also reaffirm the ability of GBT-Argus for studying dense core velocity structure and provided an empirical basis for future studies that address the angular momentum problem with a statistically broad sample.


INTRODUCTION
In molecular clouds (MCs), multi-scale supersonic flows compress material to initiate the creation of overdense regions, which may shrink to become prestellar cores and then collapse gravitationally to form protostellar systems (Shu et al. 1987).Turbulence is therefore considered one of the key agents affecting the dynamics of the star forming process in MCs, in combination with magnetic fields and gravity, at all physical scales and throughout different evolutionary stages (McKee & Ostriker 2007;Girichidis et al. 2020).There have been several observational projects aiming to resolve the gas dynamic properties at larger (MCs; e.g., the CARMA Large Area Star Formation Survey (CLASSy), Storm et al. 2014; the Green Bank Ammonia Survey (GAS), Friesen et al. 2017) and smaller (protostellar disks; e.g., the VLA/ALMA Nascent Disk and Multiplicity (VANDAM) ★ E-mail: cheyu.c@gmail.comSurvey, Tobin et al. 2016b; the Mass Assembly of Stellar Systems and their Evolution with the SMA (MASSES) survey, Stephens et al. 2018) scales.These observations revealed features that generally agree with the picture that multi-scale turbulence is dynamically significant during the star formation process, including misalignment between outflow direction and core-scale magnetic field orientation (Hull et al. 2014;Xu et al. 2022), misalignment between disk orientation and filament structure (Stephens et al. 2017), the correlation between dense core location and turbulence structure (Chen et al. 2019b), and the misalignment between core elongation and local magnetic field orientation (Chen et al. 2020a;Pandhi et al. 2023).On the other hand, the velocity information within dense cores, especially on the crucial scales of 0.01 − 0.05 pc that contain the bulk of the material destined to be incorporated into disks and subsequently accreted onto the star or planets, remains less well-investigated.
As the immediate precursors of protostars, dense molecular cores possess critical information about the physical properties at the ear-liest stage of star formation.In particular, rotation within dense cores is important in the evolution leading to the creation of protostellar systems (see review in Li et al. 2014).The angular momentum of star-forming cores is hence a critical quantity in shaping the outcome of core collapse: whether a single star or multiple system is formed (Offner et al. 2023), and whether a large or small disk is produced (Li et al. 2014).It is also important for determining whether the angular momentum of core material is conserved during the collapse into a star-disk system or whether it is significantly reduced; the latter would point to an efficient mechanism of angular momentum removal, such as magnetic braking (see review in Li et al. 2014).
While there have been several proposals on the origin of core-scale angular momentum, including the cloud-scale turbulence (Burkert & Bodenheimer 2000;Chen & Ostriker 2018), the dynamic interaction between local clumps (Kuznetsova et al. 2019), and filament fragmentation (Misugi et al. 2019), it is still not well understood through observations, because this scale (∼ 7 − 60 ′′ in nearby clouds with distance  ∼ 150 − 300 pc) is generally beyond the resolution of single-dish telescopes for low-J molecular transitions in wavelengths ∼ 3 mm regime that trace the cold dense gas, while most of the interferometers cannot cover such a large spatial area with the necessary velocity resolution within a reasonable amount of time.Because of these constraints, linear fitting is generally applied to observed velocity gradients across cores, regardless of the complex nature of the velocity field.It is assumed that rigid-body rotation applies and that the angular speed is roughly the gradient of line-ofsight velocity (Goodman et al. 1993).
It is known that some dense cores show clear gradients in lineof-sight velocities, while others have a relatively random velocity field (see e.g., Goodman et al. 1993;Caselli et al. 2002).Previous observations (e.g., Goodman et al. 1993;Caselli et al. 2002;Pirogov et al. 2003;Chen et al. 2007;Tobin et al. 2011;Chen et al. 2019c;Pandhi et al. 2023) as well as fully-3D MHD simulations (Chen & Ostriker 2015, 2018) have found a power-law relationship,  ∼  1.5 , between the total specific angular momentum , defined as the total angular momentum  divided by mass , and radius  for dense cores/clumps with radii ∼ 0.005 − 10 pc (see e.g., Li et al. 2014;Chen & Ostriker 2018).More recently, higher angular resolution maps enabled the derivation of the specific angular momentum radial profile,  (), which shows that this radial dependence is closer to  () ∝  1.8 (e.g., Pineda et al. 2019;Gaudel et al. 2020).The  −  correlation over such a large range of spatial scales suggests that either gas motion in cores originates at scales much larger than the core size, or the observed rotation-like features arise from sampling of turbulence at a range of scales (Burkert & Bodenheimer 2000).
However, observational statistics are quite poor on the crucial subcore scales of ∼ 0.01−0.05pc due to the aforementioned instrumental constraints.For example, the recent large-scale NH 3 survey GAS (Friesen et al. 2017) was able to reveal dense gas dynamics in all the northern Gould Belt clouds with   > 7 mag, but has very limited resolution at sub-pc scale due to the 30 ′′ beam.The previous interferometer survey CLASSy (Storm et al. 2014) was able to resolve N 2 H + down to ∼ 7 ′′ scale in three nearby star-forming clouds, but it lacks the necessary spectral resolution for the fine velocity structure within cores.More discussions on recent observational surveys can be found in the review by Pineda et al. (2023).
Here, we present the first results from the Dynamics in Starforming Cores (DiSCo) Survey with the recently commissioned Argus focal plane array on the Green Bank Telescope (GBT).Unlike interferometers which have limited spatial frequency coverage, GBT-Argus has full spatial sampling for both compact and extended emission.With the unique combination of angular and velocity resolution, full spatial-scale recovery, and imaging speed, GBT-Argus is ideal to survey a large sample of dense cores in nearby star-forming MCs with high resolution and sensitivity.Built on a successful pilot study (Chen et al. 2019a), the DiSCo survey aims to provide an unprecedented observational view of how the original core material forms and falls into the protostellar system.
The outline of this paper is as follows.We introduce the DiSCo project in Section 2, including detailed descriptions on target selection (Section 2.1), observations (Section 2.2), and our data reduction and spectral line fitting routines (Section 2.3).The cleaned results are presented in Section 3, where we discuss our core identification method (Section 3.1) and the derived core-scale angular momentum (Section 3.2).Further discussions are presented in Section 4.1, and we summarize our results in Section 5.

PROJECT OVERVIEW
DiSCo is a GBT large program (GBT20A-322, PI: C.-Y. Chen) aiming to systematically characterize the internal velocity structures of young dense cores in star-forming regions.The Perseus Molecular Cloud was chosen because of its relatively-close distance ( ≈ 300 pc; Zucker et al. 2019), a wide range of levels of star formation activity (e.g., Storm et al. 2014), and excellent ancillary data (see e.g., Tobin et al. 2016a;Stephens et al. 2018;Chen et al. 2019b;Hsieh et al. 2019).

Target Selection
We construct a target catalog containing a sample of 108 cores in Perseus.Based on previous N 2 H + observations (e.g., Storm et al. 2014) and our pilot studies (Chen et al. 2019a), we decided to focus on Class 0 protostellar cores because cores at later evolutionary stages (Class I and beyond) tend to have less N 2 H + emission (see e.g., Hsieh et al. 2019) .We used the catalog from the VANDAM Survey (Tobin et al. 2016a) for the Class 0 protostars in Perseus.Among the 45 Class 0 protostars in Perseus, 21 of them are single stars, while the remaining 24 are either a compact multiple system itself, or part of a multiple system.These are all included in the DiSCo project.There are 8 protostars in B1 and 21 in NGC 1333.
For starless cores, we cross-compared and combined 3 separate catalogs to achieve a complete census.Most of the cores are based on the N 2 H + dendrogram leaves identified in CLASSy (Storm et al. 2014)) which provided the most direct comparison to the DiSCo project.However, CLASSy only covers the densest areas in three of the six sub-regions (B1, NGC 1333, and L1451) in Perseus (see Fig. 1).Since 1-to-1 correlation is often found in N 2 H + and NH 3 cores (e.g., Johnstone et al. 2010), we also consider NH 3 dendrogram leaves identified in GAS (Friesen et al. 2017)) which covers all   > 7 mag sub-regions in Perseus.To make sure that we include isolated cores as well, dense cores identified in JCMT-SCUBA 850 m observation with bright N 2 H + emission detected by the 30 m IRAM telescope are also included (Kirk et al. 2007;Johnstone et al. 2010).
From the dendrogram analysis, we only consider cores with minor axes ≳ 20 ′′ (i.e. at least two beams across the narrowest region within the core) and aspect ratio ≲ 2.5 (to avoid highly elongated structures like filaments).To limit the observational time to a reasonable amount, we select only bright cores with peak N 2 H + intensities ≳ 2.0 K.We also checked the YSO catalogs from Jørgensen et al. (2007) and the VANDAM survey to make sure these cores are starless.63 starless cores are identified as DiSCo targets, among which there are 13 starless cores in B1 and 24 in NGC 1333. Figure 1 shows an annotated Herschel column density map (resolution 36.3 ′′ ; André et al. 2010;Pezzuto et al. 2021) illustrating the locations of our targets.

Observations
Argus is a 16-pixel focal plane array operating in the 74 − 116 GHz range on the GBT (Sieth et al. 2014), which is designed for efficient and sensitive large-area mapping.As GBT provides a beam size of ∼ 9 ′′ at 90 GHz, the 4 × 4 Argus array with each receiving pixel separated by 30.4 ′′ on the sky is able to significantly improve the mapping speed.Also, Argus has a low system temperature by using advanced Monolithic Millimeter-wave Integrated Circuit (MMIC) technology, which gives receiver noise temperatures of less than 53 K per pixel.
The observations presented here were conducted using the GBT from semester 2020A through 2022A with the on-the-fly (OTF) method (see Mangum et al. 2007 and references therein).The total observing time is about 116 hr.We started each session with observations of a calibrator to adjust the telescope surface for thermal corrections, to determine pointing corrections and receiver focus.We used the Argus receiver to map N 2 H + J=1-0 emission with the VEGAS backend, which was configured to a rest frequency of 93173.704MHz for N 2 H + using mode 6 with 187.5 MHz of bandwidth and 1.43 kHz (∼ 0.0046 km s −1 ) spectral resolution.System temperature calibration for all 16 receiver pixels was done before observing the science targets, and we performed calibration, pointing, and focus scans every 30-50 minutes depending on the weather.Based on our previous pilot study, we set the integration time per beam for the 16 Argus beams to be 25 − 27 seconds in order to achieve our desired sensitivity of ∼ 0.05 K in brightness temperature.We keep the sampling rate ∼ 1.2 − 1.8 ′′ per sample so that it is less than 1/3 of the beam size (∼ 9 ′′ ).This means the integration time (data output rate) is either 1 second with a map scan rate of ∼ 1.2 − 1.8 ′′ per second ( int = 1) or  int = 2 second with a scan rate ∼ 0.6 − 0.9 ′′ per second.Frequency switching was used with offsets of −12.5 and +12.5 MHz.The calibrated beam size is 9.4 ′′ for N 2 H + J = 1-0 emission.

Data reduction and spectral fitting
Argus uses the chopper wheel method for calibration, which is the standard procedure in mm and sub-mm spectral line observation (Kutner & Ulich 1981).The data is thus calibrated in   * scale, which corrects for atmospheric attenuation, resistive losses and rearward spillover and scattering.The GBT weather database returns the atmospheric temperature and opacity information.The main-beam efficiency of Argus is adopted to be  MB = 0.505 for N 2 H + at 93.17 GHz with a typical measurement uncertainty of ∼ 20%, and the flux uncertainties adopted in this work are statistical only.We refer the readers to the GBT-Argus calibration memo, Frayer et al. (2019), for more details.
We performed standard calibration using GBTIDL including baseline subtraction, and gbtgridder was used to make data cubes with pixel size 2 ′′ × 2 ′′ from maps per frequency channel using a Gaussian kernel.When re-gridding the data, we chose to combine 5 channels for N 2 H + to reach velocity resolution ≈ 0.023 km s −1 .The sensitivity of our final N 2 H + maps is ∼ 0.2 − 0.25 K.We used the Python package PySpecKit (Ginsburg & Mirocha 2011;Ginsburg et al. 2022) for spectral fitting, which simultaneously fits the 7 hyperfine lines of N 2 H + J=1-0 and returns the fitted centroid velocity, linewidth, excitation temperature, and optical depth.We adopted a signal-to-noise cut to peak line intensity of / > 5 for the N 2 H + data when performing the fitting, and a / > 2 mask when generating N 2 H + moment 0 maps, which are constructed by integrating emission from all 7 hyperfine structures of N 2 H + (see Figs. 2−3).

RESULTS
The observations toward B1 and NGC 1333 of Perseus are summarized in Fig. 2 and 3. Zoom-in velocity maps of individual cores are   000, 1-17 (2023) shown in Fig. 4-6 (see Sec. 3.1 below).Here, we use the protostar names from previous studies (e.g., Tobin et al. 2016a) to represent the corresponding protostellar cores (mostly Per-emb-(number) with some exceptions), and label starless cores as (shortened region name)-(number) (e.g., B1-2, N1333-5).Tables 1−3 list the measured and derived DiSCo core properties for B1 and NGC 1333 (single and multi-peak cores), respectively.Note that there are some offsets between the Herschel column density peak and the DiSCo N 2 H + integrated emission peak, but these are mostly insignificant.This could be due to the distribution of N 2 H + being not spherically symmetric, or simply caused by resolution effects, because the Herschel column density resolution is ∼ 4× worse than GBT-Argus (36.3 ′′ vs. 9 ′′ ).Some more discussions on peak offset can be found in Chen et al. (2019a, Sec. 4.2). ) and integrated N 2 H + intensity (grey contours; in the level of 1 K•km s −1 ).Green ellipses/arrows represent fitting results using Herschel column density contours, while purple ellipses/arrows showing cores identified using N 2 H + emission.Yellow stars mark the location of Class 0 protostars, while green circles represent dense cores (as N 2 H + peaks) identified in previous literature (see Sec. 2.1 for more details).Local maxima of Herschel column density (gray diamonds) and DiSCo integrated N 2 H + emission (cross) are also marked.The beam is shown in the upper left corner, and a physical scalebar of 0.02 pc is shown in the B1-2/B1-4 panel.The coordinates of protostars and N 2 H + peaks associated with starless cores are given in Table 1.  2.
MNRAS 000, 1-17 ( 2023)     2, but for multi-peak targets in NGC 1333.For each multi-peak core (separated by horizontal lines), we listed the identified N 2 H + peaks within the core (protostellar sources first, then starless sources) and the corresponding fitting results.We prioritize our core identification routine to have as few N 2 H + peaks within a core as possible, but we also provide the measurement using core boundaries identified with Herschel column density.For example, Bolo-45 and N1333-21 each has a N 2 H + -defined core boundary, but they belong to the same dense core when Herschel column density is used to define the core boundary.On the other hand, some sources are too clustered to have separate core boundaries, and we hence provide measurements inside the core boundaries that include multiple peaks (e.g., the SVS 13 region).

Core Identification
In general, gas column density is the preferred data for measuring the dense core boundary, because it is fitted from dust thermal radiation and thus is less biased by excitation conditions (e.g., gas density or temperature) than spectral lines.However, since dense cores commonly consist of cold, dense gas that favors molecular lines like N 2 H + , NH 3 , or C 18 O, these tracers are often considered good core material tracers (see e.g., GAS, CLASSy). 1 Following our pilot study (Chen et al. 2019a), we consider either the Herschel column density or N 2 H + integrated emission (moment 0) contours as the core boundaries.For protostellar cores, Herschel column density is used unless there is no closed contour available (due to resolution or clustered environment); in such cases, we considered N 2 H + integrated emission instead.On the other hand, all starless cores are defined using N 2 H + integrated emission, because very few of them are resolved in Herschel and we would like to keep our analysis as consistent as possible.
Instead of finding a background value by fitting a Gaussian shape to the intensity (as was done in Chen et al. 2019a), we adopted a similar approach as considered in the GRID-core algorithm (Gong & Ostriker 2011;Chen & Ostriker 2014, 2015;Gong & Ostriker 2015), and identify core boundaries as the biggest closed contours that only include one local maximum ("boundary contours"; see column "BC level" in Tables 1−3).These contours were drawn from the peak value of each core (see Tables 1-3) with negative stepsize Δ log( H /cm −2 ) = 0.05 for Herschel contours, and Δ N 2 H + = 0.1 (B1) and 0.2 K• km s −1 (NGC 1333) for N 2 H + contours.We note that the biggest advantage of defining cores using contours is that this method better preserves the natural shape of the core, while Gaussian fitting implicitly presumes axisymmetry in cores.
We calculated most of the core properties (including the mean lineof-sight velocity  los , mean velocity gradient ∇, and the dispersion of the velocity gradient direction  ∇ ) using these core-boundary contours, and fit such contours with ellipses (see Figs. 4-6) to estimate the sizes, aspect ratios, and orientations (column "fit PA") of the cores (Tables 1-3).Note that we already excluded small, compact cores when selecting targets, and thus even the smallest cores identified in DiSCo have diameters ≳ 14 ′′ , or 1.5 times of the ∼ 9 ′′ beam.
Among the total of 66 targets in B1 and NGC 1333, two appear to be not core-like in N 2 H + (see Fig. 7), both in NGC 1333.There is basically no N 2 H + emission associated with protostellar core Per-emb-24 (Fig. 7, left), which could suggest that this protostar is actually at a later evolutionary stage instead of a Class 0 source (see e.g., Tobin et al. 2016a).On the other hand, while being previously identified as a core using GAS data with 30 ′′ resolution, starless "core" N1333-2 is actually a ring-like structure (Fig. 7, right).These two targets were thus removed from the analysis.
In some cases, multiple protostars and/or local maxima of N 2 H + emission are too close to one another, and there is no closed contour (either  H column density or N 2 H + emission) with a single target.This is the case for Per-emb-6/Per-emb-10, B1-bN/B1-bS (see Fig. 4), Per-emb-27/N1333-8, Per-bolo-45/N1333-21, Per-emb-12/Per-emb-13/IRAS4B' (see Fig. 5), and the SVS 13 region (see Fig. 6).Also, 1 While we do not include a direct comparison between column density and N 2 H + emission in this work, some insight can be inferred from Figures 4  and 5 by comparing the elliptic core boundaries fitted from column density and the N 2 H + integrated emission contours.The column density-defined core boundaries and the brightest regions in N 2 H + mostly overlap, suggesting that N 2 H + is a good substitute of column density to trace core material.
when the second local maximum around the core area is weak or small in size, we loosen our criterion to include such substructures inside the identified core boundary.This is the case for N1333-1, N1333-13, N1333-17, N1333-22, and N1333-23 (see Figs. 5 and 6; the larger core area are marked with dark purple ellipses).All these cores are categorized as 'multiples' in our analysis.
Interestingly, if we define the normalized intensity variation across the core to be where  boundary is the intensity of the contour level corresponding to the core boundary, we find that the normalized intensity difference Δ/ peak tends to have a power-law correlation with size, as shown in Fig. 8.While both the log of  H column density and N 2 H + emission were used as  here, 2 we see that they follow a similar power-law correlation to the core size, roughly as Δ/ peak ∝ .
While neither N 2 H + emission intensity nor H 2 column density can be directly linked to gas density without further information of gas properties, we note that such correlation resembles a Bonner-Elbert sphere with a flat central region surrounded by a power-law dropoff in density distribution.This seems consistent with the expected power-law density profile of dense cores,  ∝  −  , from observations (see e.g., Adams 1991;Ward-Thompson et al. 1999;Shirley et al. 2000;Bacmann et al. 2000;Alves et al. 2001;Evans et al. 2001;Motte & André 2001;Kirk et al. 2005), theories (e.g., the singularisothermal sphere in Shu et al. 1987 and the Bonnor-Ebert sphere in Ebert 1955 andBonnor 1956), and simulations (e.g., Gong & Ostriker 2011, 2015;Offner et al. 2022).

Velocity Gradient and The 𝐽 − 𝑅 correlation
Following our pilot study, we calculate the gradient of  lsr locally at each map pixel first by averaging the difference of  lsr among nearest neighbors, then take the vector average of ∇ lsr to derive the mean gradient direction (the arrows in Figures 4-5).The length of this mean vector is then the mean linear velocity gradient at core scale, ∇.Assuming the observed velocity gradient is due to core rotation, the angular velocity of the core is simply  = ∇ (Goodman et al. 1993), which can then be combined with the fitted core radius  to derive the specific angular momentum  ≡ / ∼  2  =  2 ∇.We then compare the  −  relation from our DiSCo data with previous observations (Goodman et al. 1993;Caselli et al. 2002;Chen et al. 2007;Pirogov et al. 2003;Tobin et al. 2011;Yen et al. 2015) and simulations (Chen & Ostriker 2018), which is shown in Fig. 9 (left panel).
Clearly, DiSCo cores blend nicely in with the existing trend that roughly follows a power-law correlation  ∝  1.5 (see e.g., Chen et al. 2019a).The least-square fit of DiSCo data gives  ∝   with  = 1.55 ± 0.15, and  = 1.59 ± 0.04 when combining all data (from both observations and simulations) shown in Figure 9 (left panel).However, we note that while the core assembly roughly follows the same power law, the scatter in  for a given radius is pretty large, more than one order of magnitude.Considering the definition of  ≡  2 ∇ already guarantees sharp dependence on , this suggests that the dependence of ∇ on core size is relatively weak and has large scatter.We demonstrate this in the upper-right panel of Fig. 9, which  plots the radial dependence of ∇.This shows that the mean velocity gradient within individual cores ∇ has a noisy negative correlation with the core radius, which is the main source of uncertainty in the  ∝   power law dependence.

Solid-body Rotation?
The protostellar evolution of angular momentum at core scale is critical in shaping the follow-up formation of disks and/or multiple systems (Offner et al. 2023), and a rotating dense core (plus turbulence as velocity perturbations) is a generally-adopted initial condition in numerical studies (e.g., Seifried et al. 2013;Lam et al. 2019;Tsukamoto et al. 2020, or see references in the recent reviews by Pineda et al. 2023 andZhao et al. 2020).The most direct observable feature of solid-body rotation would be a uniform gradient of line-ofsight velocity across the core, as discussed in Goodman et al. (1993).The derivation of the specific angular momentum  ≡ / =  2 ∇ is also based on such an assumption.However, as previous observations also reported, only very few dense cores show clear, monotonic velocity gradient within them (see e.g., Caselli et al. 2002).Recent observations by Pineda et al. (2019) also revealed features of nonsolid-body rotation within dense cores.Such deviation from the idealized assumption of solid-body rotation could be the result from turbulent motion at the core scale (see e.g., Chen & Ostriker 2018).
To investigate such internal motion, we measured the dispersion of line-of-sight velocity gradient directions,  ∡ (∇,∇) , where ∡(x, y) represents the angle between the two vectors x, y.This quantity is listed as the uncertainty of  ∇ in Tables 1-3, which indicates whether or not the mean velocity gradient within the core is welldefined.By plotting it against core size (Fig. 9, bottom right), we see that larger cores tend to have more chaotic internal motion (i.e., larger dispersion in velocity gradient direction).This could be due to the turbulence cascade from the parent cloud.However, we also note that these larger, more turbulent cores are mostly protostellar cores and/or multiples, and thus their velocity fields could have been distorted by gravitational collapse or the interaction between the other cores.
We further investigate the internal velocity variation for cores from different environments and at different evolutionary stages.Fig. 10 shows the step-style histograms of the dispersion of velocity gradient direction,  ∡ (∇,∇) , within starless (blue) and protostellar (orange) cores, as well as those with multiple cores (green), for both B1 (left panel) and NGC 1333 (middle panel) regions, and for the two regions combined (right panel).For rigid-body rotation, there should be a prominent direction of velocity gradient within the core, and thus the angle dispersion of ∇ should be small.This is clearly not the case for the DiSCo targets presented here, as most of the cores appear to have angle dispersion ∼ 30 − 60 • in their velocity gradient direction.Interestingly, there are clearly more starless cores with small (≲ 30 • ) dispersion in ∡(∇, ∇) than protostellar cores, which makes sense because more severe, global collapse may have started in protostellar  cores, and thus the velocity structure would not be as smooth as in early-stage starless cores.We also note that, while there are not many samples, cores with multiple emission peaks ("multiples" in the plot) do not seem to differ significantly from protostellar cores, which is generally consistent with recent theoretical work by Smullen et al. (2020) that there seems to be very little correlation between core properties and multiplicity.
The large dispersion in velocity gradient direction indicates that solid-body rotation may not be a good assumption for core-scale gas motion.More importantly, we would like to emphasize that even if the cores show prominent, uniform velocity gradient direction, it is not guaranteed that such velocity gradient is due to rotation (Chen & Ostriker 2018, also see e.g., Chen et al. 2020b;Hsieh et al. 2021 for related discussions in filaments).To test this, we follow the analysis conducted by Chen & Ostriker (2018) and measure the angle between the fitted core major axis  core and the averaged velocity gradient direction.If the velocity gradient across the core is due to rotation, we expect it to be roughly parallel to the core's major axis, because the classical theory of star formation envisioned the core to flatten along the rotational axis (see e.g., Shu et al. 1987).The results are illustrated in Fig. 12, which shows the cumulative distributions of the rotation-core alignment angle, ∡( core , ∇).Note that we limit this angle to be within [0 • , 90 • ].We see that there is no clear correlation between the core orientation (as measured by the  (Stephens et al. 2018)) and the derived mean velocity gradient ∇ within corresponding DiSCo cores.Our results suggest there is little or no correlation between the velocity structure at core scale and the outflow orientation, with a slight preference toward a perpendicular alignment between these two, which is the case described in the classical theory (e.g., Shu et al. 1987).
major axis) and the velocity gradient direction, as the step functions are all closer to the diagonal line (complete random distribution).This is a strong evidence suggesting that the velocity gradients measured within dense cores are probably not pure rotation but a mix between turbulence and real rotation, in good agreement with the simulation results in Chen & Ostriker (2018).
In addition, we compared the velocity gradient directions in our protostellar core targets to their outflow directions reported in MASSES (Stephens et al. 2018).The results are shown as cumulative distributions in Fig. 11.There seems to be no correlation between the mean velocity gradient direction and the outflow orientation, and thus the angle between these two directions has a relatively random distribution.This is again in contrast to the idealized assumption of solid-body rotation, because the protostellar outflows tend to be launched along the rotational axis of the core, and thus the observed outflows are expected to be aligned perpendicular to the core-scale velocity gradient direction.
If solid-body rotation is not a reasonable assumption, the so-called specific angular momentum  ≡ / is in fact simply  2 ∇ = Δ, where Δ = ∇ •  is the velocity difference across the core.The  −  correlation discussed in Sec.3.2,  ∝  1.5 , thus simply reduces to Δ() ∝  0.5 .As noted in Chen & Ostriker (2018), this resembles the Larson's law of cloud-scale turbulence,   (ℓ) ∝ ℓ 0.5 .Our results therefore support the hypothesis (see e.g., Chen & Ostriker 2018;Burkert & Bodenheimer 2000) that the observed velocity gradient across dense core has a tight connection to the larger-scale gas dynamics instead of local rotation.For example, Figures 2 and 3 show that the gas dynamics are spatially coherent (i.e., continuous and smooth between cores).This suggests that when we measure the velocity gradient within cores, we may be simply taking samples of the gas dynamics at the larger scale.If the largescale gas dynamics are directly related to cloud-scale turbulence, this could be a hint of the turbulence origin of core-scale velocity gradient.More detailed investigations are needed to provide better, quantitative evidences.

Connecting to Larger Scales
While the main focus of DiSCo is on dense core dynamics at 0.01 pc scale, the full-range coverage of physical scale and the large field of view (∼ 2 ′ ×2 ′ ) of GBT-Argus also allows DiSCo observations to recover gas dynamics at cloud (0.1 − 1 pc) scales.The comparison between DiSCo and previous studies with similar coverages but lower resolutions therefore could provide great insight into the evolution of gas dynamics across different scales during star formation.
The two surveys on which our target selection was based, CLASSy and GAS, provided the most direct comparison to our DiSCo data.While a detailed study on connecting DiSCo results to GAS results is still ongoing, we have compared N 2 H + data from our pilot study (Chen et al. 2019a) to the corresponding CLASSy regions, and concluded that the velocity structure probed by GBT-Argus is roughly consistent with CLASSy results.We note that one caveat of such comparison is that the N 2 H + data from CLASSy have ∼ 5× worse spectral resolution than GBT-Argus, which may affect the spectral fitting result and the subsequent velocity analysis.
In addition to the aforementioned CLASSy and GAS, Dhabal et al. (2019) shows high-resolution NH 3 maps of NGC 1333 using GBT and VLA combined.By looking at the DiSCo N 2 H + velocity map of NGC 1333 in Figure 3 (right panel) and the GBT+VLA NH 3 velocity map (Figure 8 in Dhabal et al. 2019), one can tell that these two datasets are highly consistent.While detailed, quantitative comparison is beyond the scope of this work, we stress that our DiSCo data also shows a huge velocity difference (Figure 3 (right panel)) across the SVS 13 region toward IRAS 4 region in NGC 1333 (see Figure 8 in Dhabal et al. 2019 for definitions), which is consistent with the turbulent cell colliding region as described in Dhabal et al. (2019).
Last but not the least, Hacar et al. (2017) conducted N 2 H + observations using the IRAM 30m telescope (30 ′′ resolution) and identified velocity-coherent structures (termed "fibers") in the main region of NGC 1333.By looking at dense core locations from Figure 3 (left panel) and comparing with Figure 13 in Hacar et al. (2017), there seems to be a trend for cores to form around the intersection of fibers.This is consistent with the results reported in Chen et al. (2019b), who identified pressure-bound cores using GAS data and pointed out that dense cores tend to form at the turning point along the velocity axis in position-position-velocity space (i.e., local extremes).Further studies are needed to provide more insight into this core-fiber connection.

SUMMARY
We present the first results from the DiSCo project, which is a GBT-Argus Large Program to survey the velocity structures within all starless and Class 0 cores in Perseus using the N 2 H + J=1-0 line.Data from the B1 and NGC 1333 regions in Perseus are reported here.We characterize core boundaries using either the Herschel column density map (most of the protostellar cores) or the integrated N 2 H + emission (starless cores) from this project.Line-of-sight velocity is measured by fitting all hyperfine structures of N 2 H + J=1-0 altogether.
We further calculate the gradient vector of the line-of-sight velocity within cores.The velocity gradient amplitude is then used to derive the specific angular momentum assuming it originates from core-scale rotation.Our results are highly consistent with previous observations and simulations; the specific angular momentum  roughly follows a power-law correlation with the core radius ,  ∝   where  ≈ 1.55 ± 0.15 according to a least-squares fit.
On the other hand, we use the velocity gradient direction to examine whether or not the core-scale velocity structure can be best explained by solid-body rotation.Since most of the DiSCo cores show large dispersion in the velocity gradient direction (≳ 30 • ), this suggests that the mean velocity gradient (that we used to calculate the specific angular momentum) is not truly representative, Hence, the community needs to be careful when adopting the linear approximation on core-scale velocity structure, especially on data with insufficient spatial or spectral resolution.Also, the velocity gradient direction shows no obvious correlation with either the core orientation or protostellar outflow direction, which means the velocity structure within cores may have external origins.
The key scientific goal of the DiSCo project is to produce finelyresolved velocity information at dense core scale and use it to provide justifications on theoretical models of the early stages of protostellar evolution.When complete, DiSCo will have high-quality velocity information on more than 100 cores.This sample will provide strong statistical evidence about whether rotation is a common feature of dense cores, and whether the core-scale dynamics are generally consistent with cloud-scale motions.While this is not the complete data set and more detailed investigations are needed to better utilize the spectral resolution and dynamic range of the DiSCo data, the simple analysis conducted here shows evidence supporting the turbulent origin of the core-scale velocity structure.

Figure 1 .
Figure 1.Herschel H 2 column density map (in cm −2 ) of the actively star-forming Perseus molecular cloud over-plotted with the 108 Class 0 protostellar (stars) and starless (circles) targets of the DiSCo survey.

Figure 2 .
Figure 2. DiSCo observations of B1.Left: selected Class 0 protostellar (stars) and starless (open circles) targets of the DiSCo survey overplotted on the NH 3 moment 0 map (in K km s −1 ) from GAS (Friesen et al. 2017), with white contours showing the coverage of DiSCo observations.Middle: DiSCo N 2 H + moment 0 map of B1 in K km s −1 .Right: the fitted line-of-sight N 2 H + central velocity in km s −1 from DiSCo.The 30 ′′ GAS beam and 9 ′′ GBT-Argus beam are shown in the lower left corner of the corresponding panels.Also plotted in the left panel are the Class I protostars (blue diamonds) that are not included in DiSCo.

Figure 4 .
Figure 4. DiSCo N 2 H + maps for cores in the B1 region.Identified core boundaries (dashed ellipses) with mean velocity gradient direction (arrows) overplotted on normalized line-of-sight velocity maps (Δ los ≡  los −  los | core in km s −1 ) and integrated N 2 H + intensity (grey contours; in the level of 1 K•km s −1 ).Green ellipses/arrows represent fitting results using Herschel column density contours, while purple ellipses/arrows showing cores identified using N 2 H + emission.Yellow stars mark the location of Class 0 protostars, while green circles represent dense cores (as N 2 H + peaks) identified in previous literature (see Sec. 2.1 for more details).Local maxima of Herschel column density (gray diamonds) and DiSCo integrated N 2 H + emission (cross) are also marked.The beam is shown in the upper left corner, and a physical scalebar of 0.02 pc is shown in the B1-2/B1-4 panel.The coordinates of protostars and N 2 H + peaks associated with starless cores are given inTable 1.

Figure 5 .
Figure 5. Similar to Fig. 4 but for NGC 1333 (except core N1333-3 and those in the SVS 13 region).A physical scalebar of 0.02 pc is shown in the N1333-23 panel.The coordinates of protostars and N 2 H + peaks associated with starless cores are given in Table2.MNRAS 000, 1-17 (2023)

Figure 7 .
Figure 7. Two previously-defined cores (left: Class 0 protostar Per-24, marked by a green star; right: a starless core, marked by a green circle) with no core-like structure in integrated N 2 H + emission (in K km s −1 , colormap and black contours in the level of 1 K•km s −1 ) from DiSCo.Column density ( H in cm −2 , in the level of Δ log(  H /cm 2 ) = 0.1) contours from Herschel are plotted as light gray contours.The beam is shown in the bottom left corner in each panel, and other DiSCo targets in the fields of view are also marked.

Figure 8 .
Figure 8. Normalized intensity variation plotted against the size of the corresponding cores, from both B1 (blue symbols) and NGC 1333 (yellow symbols) regions and for both starless (circles) and protostellar (stars) cores, single or multiple (open squares).The values of Δ/ peak offsets between protostellar cores measured in Herschel log  H 2 column densities (lower values of Δ/ peak ) and cores (mostly starless) measured in N 2 H + emission.

Figure 9 .
Figure 9. Left: The  −  correlation combining DiSCo data and previous observations (solid circles) and simulations (open circles).DiSCo objects for B1 and NGC 1333 regions are marked by different colors, and various symbols represent the three categories: protostellar cores (stars), starless cores (triangles), and multiples (squares).Dotted line shows the least-square fit using all data on the plot ( ∝  1.59±0.04with  2 = 0.14), which is very close to the fit using DiSCo data only ( ∝  1.55±0.15with  2 = 0.09).Right: The mean velocity gradient (top panel) and the dispersion of the velocity gradient direction (bottom panel) within the core, as functions of the core size.

Figure 10 .
Figure 10.Cumulative distributions of the dispersion of gradient velocity direction for different types of cores (different colors; see legend) and for separate regions (left and middle) and all regions combined (right).Most cores tend to have a dispersion of ∼ 30 − 60 • in the direction of the velocity gradient, indicating there is no prominent overall velocity gradient direction.

Figure 11 .
Figure 11.Cumulative distribution of the alignment between the protostellar outflows measured in MASSES(Stephens et al. 2018)) and the derived mean velocity gradient ∇ within corresponding DiSCo cores.Our results suggest there is little or no correlation between the velocity structure at core scale and the outflow orientation, with a slight preference toward a perpendicular alignment between these two, which is the case described in the classical theory (e.g.,Shu et al. 1987).

Figure 12 .
Figure12.Cumulative distributions of the rotation alignment, ∡ ( core , ∇), for protostellar (orange), starless (blue), and multiple (green) cores as well as all cores combined (black dashed lines) for both B1 (left panel) and NGC 1333 (middle panel) regions and the two regions combined (right panel).Most of the step functions are very close to that of a random distribution (gray dotted diagonal line), which suggests that there is no preferred alignment between the core and the velocity gradient direction within it.

Table 1 .
Summary of the targets in the B1 region.

Table 2 .
Johnstone et al. (2010)ellar cores (upper half in the table) are the protostar coordinates from the VANDAM survey(Tobin et al. 2016a).Starless core (lower half in the table) coordinates are the + source number inKirk et al. (2007)andJohnstone et al. (2010).Similar to Table1, but for single-peak targets in the NGC 1333 region.
♦The corresponding N 2 H * ‡ Here, the uncertainty of  ∇ is estimated as the circular mean of  ∇ −  ∇ .See Sec.4.1 for more discussions on  ∡ (∇,∇) .§B1-NE is too close to the B1-b region and thus is not resolved in Herschel.MNRAS 000, 1-17 (2023)

Table 3 .
Similar to Table