CARMA LARGE AREA STAR FORMATION SURVEY: PROJECT OVERVIEW WITH ANALYSIS OF DENSE GAS STRUCTURE AND KINEMATICS IN BARNARD 1

We mosaicked an approximately 8' × 20' area of Barnard 1 using CARMA. The total observing time (150 hr) was split evenly between DZ and EZ configurations, which have projected baselines from about 1–40 kλ and 1–30 kλ, respectively. DZ configuration observations occurred in the spring and fall of 2012, and EZ configuration observations occurred in the summer of 2012. Table 1 provides a summary of the observations and calibrators. The "Z" in the DZ and EZ configurations refers to use of the 23-element observing mode (CARMA23), which utilizes cross-correlations from all 23 antennas: six 10.4 m antennas, nine 6.1 m antennas, and eight 3.5 m antennas. (Standard CARMA15 observing only includes the 10.4 m and 6.1 m antennas.)

CARMA23 provides more baselines compared with the CARMA15 mode–this enhances imaging capabilities by filling in more of the uv-plane. It also offers shorter baselines, which are important for two reasons: (1) more extended emission, if present, is recovered, and (2) the increased uv-coverage improves the link between the zero spacing information from single-dish and the interferometer visibilities when doing joint deconvolution.

Our CARMA23 mosaic of Barnard 1 is made up of three adjoining rectangular regions, each at a position angle of 40° east of north; it contains 743 individual pointings in a hexagonal grid with 31'' spacing. See Figure 3 for a depiction of the pointing centers overlaid on a Herschel 250 μm image (André et al. 2010). The reference position of the map, which is the center of the northern rectangle, is at α = 03^h33^m20^s, δ = 31°08'45'' (J2000); it encompasses the Barnard 1 main core with the B1-b and B1-c continuum cores. The other two rectangles extending toward the southwest were chosen to follow the Herschel dust emission.

During each pass through the map, we integrated for 15 s on each mosaic position before moving to the next mosaic position. We used a newly commissioned "continuous integration" technique that improves on-source efficiency (on-source integration time compared to wall clock time) of large CARMA mosaics from about 35% to 48%. The improvement comes from continuously taking data, even while antennas slew to a new mosaic position; the data for a given antenna are automatically blanked while it is slewing. This technique removes the built-in software delay associated with a slew to a new target or mosaic position. The gain in on-source efficiency is critically important for large mosaic observations, with many moves and short integration times between moves. This method should not be confused with on-the-fly mosaicking, where the antennas are continuously moving while taking data.

In CARMA23 mode, the correlator has four spectral bands in the upper sideband. The lower sideband is not present on the 3.5 m telescopes; it is present but not used on the 10.4 and 6.1 m telescopes. We tuned to the J = 1 → 0 transitions of N₂H⁺, HCO⁺, and HCN, and used a 500 MHz wide band for calibration and continuum detections. Table 2 has an overview of the correlator setup. The molecular lines were each placed in 8 MHz bands, which have 159 channels that are 0.049 MHz wide in three-bit mode; this corresponds to ∼0.16 km s⁻¹ velocity resolution near 90 GHz. All hyperfine components of N₂H⁺ and HCN fit within the 8 MHz bands. However, there is high-velocity HCO⁺ and HCN emission from outflows that lies outside our bandwidth.

Table 2. Correlator Setup Summary

Line	Rest Freq.	No. Chan.	Chan. Width	Vel. Coverage	Vel. Resolution	Chan. rms	Synth. Beama
	(GHz)		(MHz)	(km s−1)	(km s−1)	(Jy beam−1)
N2H+	93.173704	159	0.049	24.82	0.157	0.14	76 × 65
Continuum	92.7947	47	10.4	1547	33.6	0.0013	80 × 62
HCO+	89.188518	159	0.049	25.92	0.164	0.12	78 × 68
HCN	88.631847	159	0.049	26.10	0.165	0.12	79 × 68

Note. ^aIn principle, the synthesized beam is slightly different for each pointing. Miriad calculates a synthesized beam for the full mosaic based on all of the pointings.

Download table as:  ASCIITypeset image

Two team members independently inspected, flagged, and calibrated each observing track during the 150 hr campaign using Multichannel Image Reconstruction, Image Analysis and Display (MIRIAD; Sault et al. 1995). The flags for each track were then copied to our MIRIAD Interferometric and Single-Dish (MIS; Teuben et al. 2013) pipeline that has the capability to calibrate all of the tracks with minimal user input and create combined visibility data sets and maps of all the observed sources. Standard interferometric calibration steps (e.g., passband, gain, flux calibration) were performed in the MIS pipeline calibration code. We observed a nearby quasar every 16 minutes for gain calibration; 3C84 was the main gain calibrator, and 3C111 was used when 3C84 rose above 80° elevation. 3C84 doubled as a passband calibrator. We observed Uranus for absolute flux calibration during each track and used the bootflux task to determine the absolute flux of the gain calibrators. The flux of 3C84 varied between 16.0 and 19.4 Jy over the course of the campaign, while 3C111 varied between 2.8 and 4.5 Jy. The uncertainty in absolute flux calibration in the 24 combined data sets is about 10%; hereafter, we only report statistical uncertainties in quoting errors in measured values.

Simultaneously with the interferometric observations, we obtained CARMA total power observations to recover the line emission resolved out by the interferometer (i.e., CARMA in standard interferometric observing mode). CARMA's single-dish mode utilizes the autocorrelation capability of the correlator and intersperses observations of an emission-free region between on-source integrations. A total power spectrum can then be constructed from knowledge of the system temperature (T_sys), the emission region (the ON position), and the emission-free region (the OFF position).

The OFF position for Barnard 1 was in a gap of ¹²CO and ¹³CO emission 28' west and 15 south of the Barnard 1 mosaic reference position. A hole in lower density CO gas ensured that there was no significant dense gas in that region. We integrated on the OFF position for 30 s every 3.5 minutes if the atmospheric opacity was stable on the 3.5 minute timescale; otherwise, we observed purely in interferometric mode. Fourteen of 24 tracks with passing weather grades were observed in single-dish mode. We flagged the autocorrelation data separately from the cross-correlation data to account for cases where the opacity deteriorated in a track that was started in single-dish mode.

We calibrated the autocorrelation data from each antenna in MIRIAD with two steps: (1) sinbad calculated T_sys*(ON-OFF)/(OFF), and (2) sinpoly removed first-order polynomial baselines from the sinbad spectra. We averaged OFF scans on both sides of an ON scan in the sinbad calculation to reduce noise introduced from OFF scans.

We converted the calibrated single-dish spectra from each of the six 10.4 m antennas to six single-dish data cubes using the varmaps routine. We only used the 10.4 m antennas because they have the highest angular resolution and hence best overlap with the 3.5 m baselines in the uv-plane. The 10.4 m antennas have a halfpower beamwidth of 773 (N₂H⁺), 807 (HCO⁺), and 812 (HCN). The routine gridded each data cube onto 10'' × 10'' cells, and used a 50'' smoothing beam to calculate the emission in each cell according to the emission at each mosaic pointing. The size of the smoothing beam was determined empirically; a smaller beam increased the final noise in the maps, while a larger beam smoothed out structure seen in other single-dish maps of this region. The final halfpower beamwidth of the calibrated single-dish cubes is the quadrature sum of the original halfpower beamwidth and the smoothing beam: 921 (N₂H⁺), 949 (HCO⁺), and 954 (HCN).

We scaled the data cubes from all six antennas to a single reference antenna to account for systematic differences of ∼10% arising from physical differences in each antenna and from differences in bandpass shape of each antenna response. We calculated the mean of the six data cubes using imstack to improve the signal-to-noise ratio and limit antenna-based artifacts. The antenna temperature rms values in the final cubes are 0.025, 0.027, and 0.026 K for N₂H⁺, HCO⁺, and HCN, respectively.

The final data product for each observed molecular line transition is a spectral line cube produced from a joint deconvolution of the interferometric and single-dish data. The joint deconvolution was done in MIRIAD as summarized below.

We created the interferometric dirty cube and dirty beam from the calibrated visibility data set using invert with system temperature and antenna gain weighting, and Briggs's robustness parameter (Briggs et al. 1999) of −0.5. We de-selected baselines connecting 10.4 m and 3.5 m dishes due to the illumination of the first negative sidelobe of the 10.4 m beam by the 3.5 m beam. The dirty cube was then cleaned with mossdi, a Steer CLEAN algorithm, and the clean components were carried over to the joint deconvolution to aid in the convergence to a solution. The single-dish cube was regridded to the interferometric cube axes, and converted to Jy beam⁻¹ units with a 65 Jy beam⁻¹ K⁻¹ scaling factor.¹⁷

The joint deconvolution was done with a maximum entropy algorithm in MIRIAD, mosmem, that used the interferometric dirty cube, single-dish cube, interferometric dirty beam, single-dish beam, and interferometric clean components to solve for the maximum entropy model components. The final cubes are the dirty cube, minus the maximum entropy model components convolved by the dirty beam, plus the maximum entropy model components convolved by the synthesized beam. The noise levels and synthesize beams for the final data cubes are given in Table 2.

A 3 mm continuum map was created from the interferometric data in the 500 MHz window. We created and cleaned the dirty map in MIRIAD using invert and mossdi, and restored with a synthesized beam of 80 × 62. The rms in the calibrated continuum map is ∼1.3 mJy beam⁻¹.

The autocorrelation data from the 500 MHz band were not used in continuum mapping because single-dish continuum observations require fast chopping between on-position and off-position to cancel out variable sky emission; this observing technique is not available at CARMA.

The MCZ, shown in Figure 7, is the area of strongest dust emission seen by Herschel; it contains the largest cluster of young protostellar objects, and is the only region where we have detected compact continuum sources (see Section 3) and outflows (see Figure 8). The N₂H⁺ emission in this zone follows the overall structure of the dust, although our improved angular resolution reveals more detail than ever before. The HCN and HCO⁺ emission not associated with outflows is much weaker compared to N₂H⁺, and it does not follow the morphological structure of the N₂H⁺ or the dust. Figure 7 shows the locations of Spitzer YSO candidates (Jørgensen et al. 2006), Bolocam 1 mm clumps (Enoch et al. 2006), and SCUBA 850 μm clumps (Hatchell et al. 2005).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The spectra from this zone are the most complex in the Barnard 1 field—Figure 6 shows a sample N₂H⁺ spectrum from near the B1-b core that exhibits broad lines, while the HCN and HCO⁺ spectra from the same location show evidence of red-wing outflow emission overlapping with B1-b. The N₂H⁺ emission toward B1-b shows two resolved peaks. There are emission peaks in HCN and HCO⁺ toward B1-b South, but not B1-b North (see Figure 7). This agrees with findings from Huang & Hirano (2013) and supports the idea that carbon-bearing species are depleted around B1-b North, suggesting it is the younger source of the two.

The N₂H⁺ gas emission is very strong at the location of B1-c, while HCO⁺ is only detected in outflows; this agrees with observations from Matthews et al. (2006). HCN gas does exist at the location of B1-c, but is not a strong peak relative to the rest of the main core, like we see with N₂H⁺. HCN is also detected in outflow emission.

There are seven Spitzer Class 0/I YSOs in this zone. Two are known to be associated with the B1-c and Per-emb 30 compact continuum cores, which are sites of strong N₂H⁺ emission. There are four nearby the larger scale B1-a, B1-b, and B1-d dust clumps, which are also areas of strong N₂H⁺emission. The seventh YSO, which lies west of the strongest dense gas in this zone, only has weak dust and gas emission associated with it. Several of these YSOs are driving outflows, which are discussed later in this section, and identified in Figure 8 and Table 4.

There are five Bolocam 1.1 mm cores and six SCUBA 850 μm cores in this region. All of these lower resolution dust cores are near N₂H⁺ emission peaks, but not necessarily near HCN or HCO⁺ peaks. There are several peaks of HCN and HCO⁺ that are not associated with strong dust emission.

Figure 8 shows HCO⁺ and HCN dense gas outflows, which are only detected in this zone. In the HCO⁺ map, SSTc2d J033317.9+31092 (B1-c) is clearly driving a bipolar outflow. We detect another bipolar HCO⁺ outflow, likely associated with the IRAS 03301+3057 YSO based on a similar detection in CO by Hirano et al. (1997). There is a red outflow just to the west of SSTc2d J033314.4+310711, another to the north SSTc2d J033320.3+310721, and a third to the southeast of SSTc2d J033327.3+310710.

In the HCN map, we again see the bipolar outflow from SSTc2d J033317.9+31092 (B1-c), although the HCN additionally traces an extension of the red part of the outflow to the west, and an associated knot of emission to the east. The western extension flares out into two separate, high-velocity outflows that extend beyond the edge of our spectral window. The eastern knot is likely a bright bow shock from the edge of the outflow interacting with cloud material. We confirmed the association of these extended features with the main outflow using Spitzer images that trace the shocked gas. The outflows from IRAS 03301+3057 and SSTc2d J033320.3+310721 are also detected in HCN, while those from J033314.4+310711 and J033327.3+310710 are not.

The ridge zone, shown in Figure 9, lies southwest of the main core. The most striking feature is the backbone of the B1 Ridge that appears as a structure ("Ridge" in the figure) that measures approximately 260'' long by 80'' wide in the N₂H⁺ map. This corresponds to about 0.30 pc long by 0.09 pc wide at a distance of 235 pc. This structure has several Bolocam and SCUBA cores running along its spine, which are strongly clustered in the northeast half where gas emission from all molecules is the strongest. A single Class 0/I YSO sits at the northeast edge of the structure—none of the dust cores along the structure have associated protostars, and they are considered pre-stellar. There is a single Bolocam core at the southwestern end of the structure, where the Herschel map peaks in an area with little molecular emission; multi-wavelength Herschel images show that this region is brighter at shorter wavelengths compared to the northeastern section of the structure. This means it could be a region of slightly higher temperature, which would explain a lower abundance of dense molecular gas. Figure 6 shows sample spectra from along this structure.

Zoom In Zoom Out Reset image size

A second molecular feature in this zone is a newly discovered filament that runs parallel to the main ridge ("Narrow Fil" in Figure 9). This filament is extremely narrow and is offset from the rest of the gas in velocity space by about 1.5 km s⁻¹. Figure 10 shows molecular line contours overlaid on a Herschel 250 μm map. The filament is about 20'' wide and 25 long, with a small kink half-way along its length. At a distance of 235 pc, this corresponds to 0.022 pc wide by 0.17 pc long. The peak integrated intensity of the filament is 1.2 Jy beam⁻¹ km s⁻¹, 1.4 Jy beam⁻¹ km s⁻¹, and 0.9 Jy beam⁻¹ km s⁻¹, for N₂H⁺, HCN, and HCO⁺, respectively. It is not possible to identify the filament in Herschel maps because it is extremely narrow and lies along the same line of sight as the western edge of the B1 Ridge. We briefly discuss the kinematics and possible formation mechanisms of the filament in Section 5.

Zoom In Zoom Out Reset image size

The third molecular feature is gas in the northern part of the zone that provides a tenuous link between the main ridge structure and the MCZ. A Class 0/I YSO is located off the southwestern edge of this feature, but there are no dust condensations directly associated with it.

The SW clumps zone lies southwest of the B1 Ridge, and is shown in Figure 11. There are five distinct clumps of N₂H⁺ emission, with a sixth clump only detected in HCN and HCO⁺. There are three Class II YSOs in this zone that do not correlate with any gas or dust peaks, which is expected for more evolved protostars. There is one Class 0/I YSO in this zone that lies at the northeast tip of the western-most clump. The strongest integrated HCO⁺ emission across the entire Barnard 1 field is located in this clump—Figure 6 shows sample spectra from it. The gas morphology of this zone is very different from the MCZ and ridge zone. Those zones have dense gas peaks embedded within lower-intensity, larger-scale dense gas structures, while the dense gas peaks in this zone are not joined at weaker emission levels.

Zoom In Zoom Out Reset image size

We used our new non-binary dendrogram code to identify gas structures traced by the isolated hyperfine component of N₂H⁺. We chose to use the isolated hyperfine component of N₂H⁺ since it is sufficiently separated from other components in velocity space to prevent contaminating our object identification along the velocity axis. Before running the data through the dendrogram code, we binned the cube by two velocity channels to improve the signal-to-noise ratio. No velocity information was lost since there are no lines of sight that contain two or more independent structures within 1.0 km s⁻¹ of each other. We found that binning by two channels provides the most improvement to signal-to-noise ratio without biasing the maps toward wide-line regions. The 1σ sensitivity of the binned data cube is 0.094 Jy beam⁻¹.

Our non-binary dendrogram code takes the same input parameters as the standard IDL implementation discussed in Rosolowsky et al. (2008b). We ran it with the following critical inputs and parameters: (1) a masked cube containing all pixels greater than or equal to 4σ intensity, along with adjacent pixels of at least 2.5σ intensity, (2) a set of local maxima greater than or equal to all their neighbors in 10'' × 10'' by three channel (0.94 km s⁻¹) spatial-velocity pixels, (3) a requirement that a local maximum must peak 2σ above the intensity where it first merges with another local maximum for it to be considered a leaf, (4) a requirement of at least three synthesized beams of spatial-velocity pixels for a leaf to be considered real.

The non-binary dendrogram shown in Figure 15 contains 41 leaves and 13 branches. The vertical axis of the dendrogram represents the intensity range of the pixels belonging to a leaf or branch. The horizontal axis does not normally carry physical meaning, but in this case we arranged the leaves and branches according to zone. Figure 16 shows two-dimensional representations of all leaves overlaid on the N₂H⁺ integrated intensity map; these representations were created by integrating over the velocity axis of the leaves and contouring the maximal R.A.–decl. extent of the integrated emission. The leaves that peak at least 6σ in intensity above their first branch are colored green (leaves 10, 25, and 39). Leaf 25 is the strongest, with a peak intensity of 1.91 Jy beam⁻¹ (in the binned data cube), and represents the gas around B1-b. Leaf 39 is the next strongest at 1.77 Jy beam⁻¹, and represents gas around B1-c. Leaf 10 peaks at 1.04 Jy beam⁻¹, and represents gas around the newly detected compact continuum source within Per-emb 30. As a reminder that the identification of leaves and branches was done in three dimensions, Figure 17 shows five N₂H⁺ velocity channels near the B1-b cores with the isocontours of leaves 24, 25, and 40, and branch 41, identified with distinct colors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We are using the 6σ contrast criteria to highlight the strongest leaves before we are able to confidently discuss the virial boundedness of actual "cores" and "clumps"—defining bound cores requires accurate measurements of mass to go along with our high-resolution kinematic information, which is beyond the scope of this paper. We chose a 6σ contrast to highlight the gas in Barnard 1 that is located near existing compact continuum objects—this gas is likely bound. We will apply the same contrast cut to other CLASSy regions. While we cannot discuss boundedness of individual gas cores at this stage, we can discuss the properties of the strongest, highest-column density, dendrogram features and compare them to weaker features in the field. Leaves with 6σ or greater contrast will be referred to as high-contrast leaves later in the paper, and the rest of the leaves will be called low-contrast leaves.

The tree structure presented here captures the qualitative hierarchical nature of the dense gas in Barnard 1, but it does not give quantitative information about the hierarchy or physical properties of the leaves and branches. Doing science with a dendrogram requires extracting tree statistics and physical properties of the leaves and branches, such as size, axis ratio, linewidth, mass, luminosity, etc. We present three simple tree statistics and the spatial and kinematic properties of leaves and branches in this paper. In future papers, we will perform more detailed comparisons between dendrogram properties across multiple CLASSy regions.

Computation of tree statistics from a non-binary dendrogram is one method of quantifying a tree structure. Houlahan & Scalo (1992) were the first to develop tree statistics for the astronomy community in order to quantify the hierarchical nature of molecular cloud structure (see their discussion of "merged trees," which are analogous to our non-binary trees). They describe several statistics that can be computed from non-binary trees; we describe three of the simplest statistics here.

The maximum branching level is the integer number of branching steps needed to reach the leaves at the top of the hierarchy. A region that has undergone a lot of hierarchical fragmentation will have more levels than a region that has not fragmented, or that has fragmented in a single step into several pieces with similar properties. Branching levels are defined upward from the base of the tree; the base (0 Jy beam⁻¹ intensity) is considered level 0. An isolated leaf that sprouts directly from the base is at level 0, while a leaf that grows from a branch one level above the base is at level 1. The maximum branching level in our Barnard 1 N₂H⁺ dendrogram is level 4, where leaves 24, 25, and 40 merge together into branch 41 (see Figure 15). These leaves are in the MCZ, which is the only zone with detections of 3 mm compact continuum cores and Class 0/I YSOs. It is reasonable that the largest amount of fragmentation has occurred in the region with the most young star formation activity. That being said, the ridge zone and SW clumps zone have maximum branching levels of three and one, respectively, so we are dealing with small variations of this statistic within Barnard 1.

A related tree statistic is the mean path length of all segments in a dendrogram. The path length of a segment is defined as the number of branching steps required to go from a leaf to the bottom of the tree (there are as many path lengths in a dendrogram as there are leaves). For example, the path length of leaf 25 is four, while the path length of leaf 37 is three—the path length of leaf 30 is zero because it directly sprouts from the bottom of the tree. The mean path length in the whole Barnard 1 N₂H⁺ dendrogram is 1.2 levels, while it is 2.0, 1.2, and 0.5 for the MCZ, ridge zone, and SW clumps zone, respectively. Again, it is logical that the MCZ, which is actively fragmenting into young stars, has the largest mean path length.

It is important to point out that these statistics would be different if we used the standard dendrogram algorithm instead of our non-binary one. The standard algorithm forces binary branching by introducing extra branching steps to ensure that only two objects merge at a time. For example, leaves 24, 25, and 40 would not all be allowed to merge into branch 41 even though that is what the data suggest is happening—leaves 24 and 25 would be merged into one branch, which would then be merged with leaf 40 into a second branch. This inflates the path length of leaves 24 and 25, and increases the maximum branching level of the tree. See Appendix B for more details about differences between the standard dendrogram algorithm and the new non-binary algorithm.

The third simple tree statistic is the mean branching ratio of all branches in a tree. The branching ratio for a single branch is defined as the number of substructures into which it fragments directly above it in the hierarchy. For example, the branching ratio of branch 41 is three, because it has three leaves directly above it. A very hierarchical region, where every object fragments into two nested objects, will have a mean branching ratio of 2, while a region that has fragmented in a single step into several pieces with similar properties will have a much larger mean branching ratio. The mean branching ratio in our Barnard 1 N₂H⁺ dendrogram is 3.9. We include the branching ratio for the base of the tree in this calculation, even though it is not explicitly defined with a branch number—this ensures that leaves that sprout directly from the base are factored in. (Note that the branching ratio above the tree base is always fixed at 2 for the standard dendrogram algorithm that forces binary branching.)

These simple statistics provide measures of the amount and type of fragmentation a given region has undergone; a region with a lot of hierarchical fragmentation will have more levels, a larger mean path length, and a smaller branching ratio, compared to a region that is just beginning to form strong overdensities. A caveat is that the absolute values of the branch statistics depend on algorithm parameters and choice of allowed branching steps (reminder that we restrict branching to intensity steps equal to integer values of the 1σ sensitivity of the data instead of allowing branching at infinitely small intensity steps). However, the statistics can be compared across trees as a differential measurement if the same algorithm parameters and allowed branching steps are used; as pointed out by Houlahan & Scalo (1992), the power of tree statistics is strongest when comparing regions with different star formation properties. We will compare the Barnard 1 tree statistics with the tree statistics of NGC 1333 and L1451 in an upcoming paper to explore whether the hierarchical complexity of dense gas represented by tree statistics correlates with the diversity of star formation activity in the western half of Perseus.

Another caveat arises from linking dendrogram hierarchical complexity to cloud fragmentation—tree statistics can be affected by projection effects. All observations suffer from projection effects when transforming the true position–position–position information of a molecular cloud to observable PPV information (for an in-depth discussion of projection effects in spectral line data, see Beaumont et al. 2013). Projection effects distort optically thick, space-filling emission more than optically thin emission concentrated in one part of the cloud. These effects are not overly concerning in our data because the isolated hyperfine component of N₂H⁺ is optically thin, and the N₂H⁺ emission is arising only from the densest parts of the cloud. If our data do suffer from projection effects, we assume that our large mapping areas create a big enough sample of dendrogram structures so that equal numbers of projection effects occur in each region we will compare.

The size of an object is one of its most fundamental properties, but defining the sizes of irregularly shaped objects in a uniform way is not trivial. This is one challenge of high-angular-resolution surveys of nearby molecular clouds—CLASSy resolves a rich morphology of structures with a range of irregular shapes.

When fitting for the size of a leaf or branch, we consider all pixels in its plane-of-sky footprint (as opposed to only considering the "onion-layer" emission of each individual branch, which is used to derive kinematic properties in the next section—see Figures 18 and 19 for examples). We use the regionprops program in MATLAB to fit an ellipse to the integrated intensity footprint of each dendrogram object. All pixels are given equal weight so that the ellipse is preferentially fit to the largest scale of the object—this prevents strong emission at the center of an object from driving the fit toward a smaller scale. The fit is defined with an R.A. centroid, decl. centroid, major axis, minor axis, and position angle. The location, major axis, minor axis, and position angle are directly determined by the fit, and are listed in Columns 2–6 in Table 5. We do not report formal uncertainties on these values since the spatial properties of irregularly shaped objects are dependent on the chosen method, and we do not want to mislead the reader as to the "certainty" of these values.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To quantify the shape of an object, we use the axis ratio and filling factor of the fit (Columns 7 and 8 of Table 5). The axis ratio is defined as the ratio of the minor-to-major axis; a circular fit will have an axis ratio of one, and a highly elongated fit will have an axis ratio approaching zero. The filling factor is defined as the area of emission within the fitted ellipse divided by the area of the fitted ellipse; a circular object and an elliptical object will both have a filling factor near one, but an irregular object will likely have a filling factor much less than one, regardless of the axis ratio of the fitted ellipse. We lastly define an object size (Column 9 of Table 5) as the geometric mean of the major and minor axes. The values reported in Table 5 have been converted to parsecs, assuming a distance of 235 pc. Figure 18 shows example spatial property fitting results for leaves 26 and 39 and branch 42.

Histograms of the size, filling factor, and axis ratio for all leaves and branches are plotted in the top row of Figure 20. The trend in size is obvious—branches represent larger structures than the leaves that peak within them. Leaf sizes range from about 0.01 pc to 0.07 pc, while branch sizes range from about 0.08 pc to 0.34 pc. The high-contrast leaves are larger than the majority of low-contrast leaves. Nearly 25% of leaves have filling factors greater than 0.8, while all branches have filling factors less than 0.8. This shows that leaves are better fit as regular shapes compared to branches, but there is still a large population of irregularly shaped leaves that are far from round or filamentary. The objects with the maximum axis ratios (leaf 14 and branch 46) are paired with filling factors below 0.8, indicating that they are not regularly shaped spheroids as the axis ratio alone might be interpreted. There are no obvious distinctions between high-contrast and low-contrast leaves in filling factor or axis ratio.

Zoom In Zoom Out Reset image size

One strength of CLASSy data is the kinematic information across a wide range of spatial scales. We can use the integrated intensity footprints of the dendrogram objects (leaf footprints shown in Figure 16) in combination with the centroid velocity and velocity dispersion maps in Figure 12 to determine the kinematic properties of leaves and branches. The four properties in Table 5 that we focus on are: the mean and rms centroid velocity (〈V_lsr〉 and ΔV_lsr, respectively), and the mean and rms velocity dispersion (〈σ〉 and Δσ, respectively).

We illustrate how we derive these properties for leaves and branches in Figure 19. For leaves, we mask the full N₂H⁺ V_lsr and σ maps with the integrated intensity footprint of each leaf, and calculate the error weighted mean and standard deviation of the masked pixels within the masked footprints. For branches, we mask the full N₂H⁺ V_lsr and σ maps with the integrated intensity footprint of each branch, excluding pixels associated with any leaves or branches inside the branch footprint. We calculate the error weighted mean and standard deviation of the remaining pixels belonging exclusively to the branch. By excluding leaf pixels, we are calculating the kinematic properties of branches only at the larger spatial scales that are not captured by the leaves within them—this method allows us to parse kinematics across different spatial scales, which enables a size–linewidth analysis in the next section. We only report kinematic properties of leaves and branches in Table 5 that have three or more synthesized beams worth of kinematic pixels; a low signal-to-noise feature may be strong enough to be selected as a dendrogram object, but still too weak to have enough kinematic data points based on our line fitting criteria discussed in Section 5.

Histograms of 〈V_lsr〉, ΔV_lsr, and 〈σ〉 are plotted in the bottom row of Figure 20. The distribution of 〈V_lsr〉 shows that most of the gas in Barnard 1 is near 6.5 km s⁻¹, while the gas from leaves 0 and 1 and branch 50 (representing the narrow filament in the ridge zone) is strongly redshifted.

There is a clear offset in the distribution of ΔV_lsr between leaves and branches: the majority of branches, which represent the largest objects, have larger ΔV_lsr than the smaller leaves. There are two effects that are likely causing this trend. The first effect is from organized velocity field motions (e.g., indicative of rotation or feedback from outflows) that exist toward a leaf and its surrounding branch. For example, the V_lsr field for the B1-c leaf and parent branch can be seen in Figure 19. The branch that directly surrounds B1-c has a larger ΔV_lsr compared to the B1-c leaf due to larger differences between the redshifted and blueshifted velocities in the rotating branch. However, the majority of leaf-branch pairs do not exhibit matching, well-ordered velocity fields. Therefore, turbulence is the second likely explanation for this trend. Turbulence has the property of causing scale-dependent spatial correlations within a velocity field—points separated by larger distances will have larger rms velocity differences between them than will points that are closer together (McKee & Ostriker 2007), meaning that branches will have larger ΔV_lsr values than leaves.

The high-contrast leaves are preferentially shifted to larger ΔV_lsr relative to the low-contrast leaves. The high-contrast leaves all have known protostellar objects within them, so their gas motions are most likely dominated by infall, rotation, or feedback from outflows, all of which increase the variation in the velocity field. The low-contrast leaves with the largest ΔV_lsr tend to have more organized centroid velocity patterns (manifested as gradients) than the low-contrast leaves with the smallest ΔV_lsr, indicating that they may represent gas around pre-stellar cores that are more likely to form stars in the future.

The distribution of 〈σ〉 differs from the distribution of ΔV_lsr. The high-contrast leaves representing the gas around B1-b and B1-c have the largest 〈σ〉 in the entire field, and the majority of branches have 〈σ〉 that are similar to those of the majority of leaves. We would naively expect that if the branches represent the largest objects projected on the plane of the sky, then they should also represent the most extended objects along our lines of sight, and should therefore have the largest velocity dispersions (assuming velocity dispersion scales proportionally with length-scale, which is the case for a turbulent cloud; McKee & Ostriker 2007). We explore these results in more detail in Section 7, where we compare the sizes of structures to their ΔV_lsr and 〈σ〉 in order to probe the physical and turbulent nature of Barnard 1.

Our dendrogram analysis identified N₂H⁺ structures in the Barnard 1 hierarchy as either leaves or branches, and we evaluated the size and kinematic properties of those objects in Section 6. Since the spatial scales of dendrogram structures in Barnard 1 range from ∼0.01 to 0.3 pc, we can utilize these structures to investigate the turbulent properties of the cloud. Before discussing size–linewidth relations that use the full angular resolution of our observations, it is useful to consider a relation that uses the "total" linewidth of each dendrogram structure—this relation can be compared with classical size–linewidth relations from Larson (1981) and Solomon et al. (1987), in which the detailed kinematic variation within individual clouds was not resolved. We calculated a total linewidth for each structure by summing all the spectra assigned to it and fitting the velocity dispersion of the resulting spectrum (the σ, not FWHM, of the line); for a visualization, see Figure 19, which shows the kinematic maps of an example leaf and branch. In that figure, all the spectra within the red contour of leaf 39, and all the spectra within the white contours of branch 42 (including the leaf emission), are summed to calculate the total linewidths of those respective structures.

Figure 21 shows plots of projected size versus total linewidth for all dendrogram structures that have kinematic values listed in Table 5, excluding the objects associated with the peculiar narrow filament in the ridge zone (leaves 0 and 1, branch 50). We separate structures in the MCZ from structures in the ridge and SW clumps zones (RSWZs) because the MCZ has a cluster of protostars driving outflows. Keeping the zones separate allows us to see whether or not they have different turbulent or physical properties. We also separately identify the few structures surrounding compact continuum cores, since their non-thermal gas motions are likely influenced by their central source.

Zoom In Zoom Out Reset image size

The data in Figure 21 are scattered and could be fit by several functions. However, there does appear to be variation of total linewidth with size in both zones, with the linewidth varying more strongly with size in the MCZ. If we use a single power-law fit, the MCZ correlation is best fit with a slope of 0.37 ± 0.08, and the RSWZ correlation slope is 0.16 ± 0.06. The steeper slope in the MCZ could be due to outflows adding energy to the turbulence at these scales. The best-fit slope from the MCZ is consistent with the result (Larson 1981) that clouds ranging in size from 0.1 to 100 pc have a power-law slope of 0.38. However, we are probing the scaling relation across much smaller spatial scales than has typically been done, and unlike previous studies, we are using high-density molecular tracers that are sensitive to non-thermal motions near the sonic scale of the cloud. Therefore, directly comparing our results to studies that used larger scales and lower-density gas tracers, such as ¹²CO, should be considered with caution.

We have two kinematic measurements of leaves and branches in Table 5 that use the full angular resolution of our data set—they complement the total linewidth and enable a different type of size–linewidth analysis. The first measurement is the V_lsr variation, ΔV_lsr. For a leaf or branch of size L, ΔV_lsr(L) is computed by measuring the rms variation of the (beam-scale) centroid velocity over the whole structure. The second measurement is the mean non-thermal velocity dispersion, 〈σ〉_nt. To calculate 〈σ〉_nt from 〈σ〉 Table 5, we remove the thermal velocity dispersion, σ_th, of 11 K N₂H⁺ gas, which is the typical gas kinetic temperature across the entire Barnard 1 region (Rosolowsky et al. 2008a): $\langle \sigma \rangle _{\mathrm{nt}} = \sqrt{\vphantom{A^A}\smash{{{\langle \sigma \rangle ^{2} - \sigma _{\mathrm{th}}^{2}}}}}$. Recall that 〈σ〉 for a structure is defined by averaging all the individual velocity dispersions over the structure; each individual pointing captures the velocity dispersion along the line of sight for a given beam-width.

We compare the dependences of ΔV_lsr and 〈σ〉_nt on projected structure size in Figure 22. Size versus ΔV_lsr for the dendrogram structures in both zones is plotted on the left; the data appear to follow a power-law relation. There is a steep correlation in both regions, with power-law slopes of 0.80 ± 0.08 in the MCZ and 0.69 ± 0.17 in the RSWZ. The value of ΔV_lsr is very low for the smallest objects in all zones, and rises to the sound speed for the largest objects in the MCZ. Size versus 〈σ〉_nt is plotted on the right side Figure 22. There is a much flatter correlation between size and 〈σ〉_nt in both zones, with a power-law slope of 0.12 ± 0.11 in the MCZ and 0.04 ± 0.07 in the RSWZ. The mean non-thermal velocity dispersions are all sub-sonic to trans-sonic in the MCZ (though still superthermal). The two structures with largest 〈σ〉_nt are surrounding B1-b and B1-c; it is likely that the non-thermal motions of the gas around them is boosted by rotation and/or outflows and is not purely turbulent. The 〈σ〉_nt are all sub-sonic in the RSWZ.

Zoom In Zoom Out Reset image size

Why are the dependences of 〈σ〉_nt and ΔV_lsr on projected structure size so different? In the following section, we set up a theoretical framework that can explain this disparity in terms of differences between a structure's projected size and its depth into the plane of the sky.

In this section, we assume that the observed non-thermal gas motions in our dendrogram structures are generated by isotropic, three-dimensional turbulence. In the previous section, we noted a few structures surrounding compact continuum cores that likely have non-thermal gas motions due to the influence of their central source; we exclude them from this analysis.

The total turbulent linewidth of a structure is most strongly influenced by the largest scale an observation encompasses, under the assumptions that the gas we observe is isotropically turbulent, and optically thin with uniform excitation conditions. If the structure is wider across the plane of the sky than it is deep into the sky, then its turbulent linewidth will be most strongly influenced by its projected size. But if the structure is deeper into the sky than in either of the projected directions, then its turbulent linewidth will be most strongly influenced by its depth.

For a structure that is many resolution elements across (the case for all of our dendrogram structures), ΔV_lsr and 〈σ〉_nt probe turbulence in different ways from the total linewidth. The value of V_lsr on any given line of sight is a measure of the mean motion of gas along that line of sight, which varies from one line of sight to another line of sight as a consequence of turbulence. The V_lsr variation, ΔV_lsr, is expected to increase with the projected scale of a structure if turbulence increases with scale, as is the case for the observed and simulated turbulent power spectra in molecular clouds. For any given line of sight, the turbulent contribution to 〈σ〉 will be set by the larger of the beam size and the structure depth. Since our beam size is very small (∼0.008 pc), we expect that the linewidth along individual lines of sight will be dominated by structure depth, and thus 〈σ〉_nt for a structure will depend primarily on the mean (unseen) depth of that structure into the sky. Statistically, along many lines of sight for a resolved structure, the comparison of linewidths set by the projected size of the structure (ΔV_lsr) to linewidths set by the depth of the structure (〈σ〉_nt) should allow an estimation of the structure depth.

We can use a mathematical framework to explain how we can compare ΔV_lsr to 〈σ〉_nt to learn about the typical depth of gas structures. In the full three-dimensional turbulence, line-of-sight turbulent velocities (v_{z, turb}) can be written as

$$\begin{equation} v_{z,\mathrm{turb}}(x,y,z)=\sum _{k_x,k_y,k_z} v_{p}(k_x,k_y,k_z) \exp [i (k_x x + k_y y +k_z z)], \end{equation} \tag{ 4 }$$

where v_p(k_x, k_y, k_z) is the turbulent velocity power spectrum, k_x = (2π/L_x)n_x, k_y = (2π/L_y)n_y, and k_z = (2π/L_z)n_z. L_x is the size of an observed structure in the x direction, and n_x is the number of full waves of a given turbulent velocity component that fit in the x-direction (the same is true for the y and z components in these equations). If the turbulent motions are Gaussian and isotropic, then v_p(k_x, k_y, k_z) is drawn from a normal distribution whose standard deviation depends only on the magnitude of the k-vector, $\vert {\bf k} \vert =\sqrt{\vphantom{A^A}\smash{{k_{x}^2+k_{y}^2+k_{z}^2}}}$. We assign the z direction as the line of sight, so that in Equation (4), v_{z, turb}(x, y, z) is the z-component of the turbulence, and there are other independent components in the x and y directions.

In this framework, the total linewidth of a structure is the sum of components with non-zero k_x, k_y, or k_z. The resolved linewidths (ΔV_lsr and 〈σ〉_nt) are explained below with the aid of Figure 23.

Zoom In Zoom Out Reset image size

The observed centroid velocity at a resolved projected position in the structure, V_lsr(x, y), is the sum of components in Equation (4) that have k_z = 0. The V_lsr dependence on spatial scale via v_p(k_x, k_y; k_z = 0) is assumed to depend statistically only on |k| if the turbulence is isotropic. Therefore, V_lsr(x,y) has the same statistical dependence on $\vert {\bf k} \vert = \sqrt{\vphantom{A^A}\smash{{k_{x}^{2} + k_{y}^{2}}}}$ (which we can measure) that the underlying turbulent velocity, v_{z, turb}(x, y, z), has on $\vert {\bf k} \vert =\sqrt{\vphantom{A^A}\smash{{k_{x}^{2} + k_{y}^{2} + k_{z}^{2}}}}$ (which we cannot measure). The V_lsr(x, y) variation, ΔV_lsr, which is the standard deviation of V_lsr(x, y) over all resolved positions in a structure, will scale with the (x, y)-size of the structure according to the turbulent scaling relation of the cloud (see Dickman & Kleiner (1985) for a discussion of spatial correlation properties of centroid velocity fields). To visualize this, we consider the structure on the left side of Figure 23. We call the full four-cell structure, structure A, and any contiguous two-cell substructure within it, structure B. Structure A is wider than structure B in the (x, y) directions on the sky, so the minimum, non-zero k_x and k_y must be smaller in structure A than structure B. Assuming v_p has a (statistical) inverse power-law dependence on |k|, this implies higher power in structure A, with a greater contribution to the V_lsr variation across the (x, y) face of structure A compared to structure B.

The non-thermal velocity dispersion at each resolved position in a structure, σ_nt(x, y), is obtained from the sum of turbulent components with all possible k_x and k_y, and non-zero k_z, along each line of sight. For a given (x, y), the total line of sight velocity is a sum: v_z(z) = v_{z, turb}(z) + v_th(z), where v_th is the thermal component of the velocity. The mean value of v_th(z) is zero, and the mean value of v_{z, turb}(z) is V_lsr. Thus, assuming that v_{z, turb}(z) and v_th(z) are uncorrelated, for a given line of sight, the square of the total velocity dispersion, σ²(x, y), is

$$\begin{eqnarray} &&\sigma ^{2}(x,y) = \left\langle v_{z}^{2} \right\rangle - V_{\mathrm{lsr}}^2 = \left\langle v_{z,\mathrm{turb}}^2 \right\rangle - V_{\mathrm{lsr}}^2 + \left\langle v_{\mathrm{th}}^2 \right\rangle + 2 \langle v_{z,\mathrm{turb}}\rangle \langle v_{\mathrm{th}}\rangle.\nonumber \\ \end{eqnarray} \tag{ 5 }$$

The last term in Equation (5) is zero because the mean value of v_th(z) is zero. The second-to-last term in Equation (5) is equivalent to the square of the thermal velocity dispersion, $\sigma _{\mathrm{th}}^2$; $\sigma ^{2} - \sigma _{\mathrm{th}}^{2}$ is equal to the non-thermal component of the velocity variance, $\sigma _{\mathrm{nt}}^{2}$, which means that the first two terms Equation (5) contribute to the non-thermal velocity dispersion. The combination of the first two terms is

$$\begin{equation} \left\langle \left\vert \sum _{k_{x},k_{y};k_{z}\ne 0} v_{p}(k_{x},k_{y},k_{z}) \exp (i {\bf k}\cdot {\bf x}) \right\vert ^{2} \right\rangle, \end{equation} \tag{ 6 }$$

where the summation includes combinations with all k_x, k_y, and with all k_z except k_z = 0; k is the wavenumber vector, and x is the spatial vector. This shows that the mean non-thermal velocity dispersion for our structures, 〈σ〉_nt, defined as the mean of σ_nt(x, y) at each resolved position, will scale with the z-depth of the structure according to the turbulent scaling relation of the cloud. Along any given line of sight, a larger L_z implies that the minimum nonzero k_z = 2π/L_z is smaller, such that the contribution to Equation (6) would be larger (under the assumption that v_p(k) statistically increases with decreasing |k|). To visualize this, we can consider the structure on the right-side of Figure 23. We call the full four-cell structure structure A, and any contiguous two-cell substructure within it structure B. Structure A extends deeper into the sky than structure B, so the minimum, non-zero k_z in structure A is smaller than the minimum, non-zero k_z in structure B. This implies higher power in structure A, with a larger contribution to σ_nt(x, y) in structure A compared to structure B.

Based on the arguments above, the scaling relation between V_lsr variation and projected size, l, (ΔV_lsr∝l^q) should have the same power-law dependence as the scaling relation between mean non-thermal velocity dispersion and depth, d, into the sky (〈σ〉_nt∝d^q) if the turbulence is isotropic and we are observing all the gas along each line of sight. But we saw in Figure 22 that the dependences of ΔV_lsr and σ_nt on projected structure size are very different. The simplest explanation of this difference is that the projected size of an object need not be the same as its line-of-sight depth into the sky. The size axis in Figure 22 (right) is an estimate based on the geometric mean of the projected size. If every object, no matter its projected size, has the same depth, then we should measure a similar non-thermal velocity dispersion for those objects. Therefore, we interpret the shallow dependence of 〈σ〉_nt with size as an indication that our collection of dendrogram structures have similar depths.

To illustrate these effects, we created a numerical realization of a three-dimensional, isotropic turbulent power spectrum, with power spectrum scaling v²∝k⁻⁴. Using an arbitrary length scale, L, the position–position–position box had x, y, z dimensions of L × L × 2L, where L corresponds to 512 pixels. We used this realization along with three sub-boxes with different L_z to demonstrate how 〈σ〉_nt and ΔV_lsr depend on structure depth for a single turbulent realization. The original box had dimensions of L × L × 2L (the "double depth" box), the first sub-box was L × L × L (the "full depth" box), the second was L × L × L/2 (the "half depth" box), and the third was L × L × L/4 (the "quarter depth" box). For each of these sub-boxes (with uniform density), we created PPV cubes and then processed them similarly to the observations. The end products were centroid velocity and velocity dispersion maps across the L × L x, y surface of each numerical box. We segmented each L × L surface into equally sized square regions with side L, L/2, L/4, L/8, and calculated ΔV_lsr and 〈σ〉_nt within each segment to show how they scale with size. The results for each box are shown in Figure 24. Evidently, 〈σ〉_nt has no dependence on size, while ΔV_lsr increases with size, similar to the behavior seen in the observations in Figure 22. The size where ΔV_lsr crosses 〈σ〉_nt corresponds well with the depth of the numerical cube.

Zoom In Zoom Out Reset image size

If turbulence in the observed region has similar behavior to that in the isotropic turbulent realization, then we can use the correlation of V_lsr variation with size to estimate the depths of the Barnard 1 structures. For this correlation, we know the projected size is directly related to ΔV_lsr, so there is no size ambiguity as there is in the correlation between projected size and 〈σ〉_nt. To estimate the typical depth of the cloud, we can find the size scale at which the best-fit line for ΔV_lsr versus projected size is equal to the best-fit line for 〈σ〉_nt versus projected size. This corresponds to ∼0.11 pc for the MCZ, and ∼0.18 pc for the RSWZ. Under the assumption that turbulence is isotropic, this implies that Barnard 1 is on the order of 0.1–0.2 pc deep, comparable to the largest projected width (∼0.2–0.3 pc) of individual structures.

Since the largest-scale dust structure in Barnard 1 extends over 1 pc on the sky, the overall region may be more flattened than spherical at the largest scales. Numerical simulations over the past 15 years have consistently shown that high-density sheets and filaments are a generic result of strongly supersonic turbulence with parameters comparable to those in observed clouds (e.g., reviews by Elmegreen & Scalo 2004; Mac Low & Klessen 2004; McKee & Ostriker 2007). But measuring the depth of an individual structure in a molecular cloud is a difficult task. Comparison to realizations of turbulence suggests that the high angular resolution, large area observations from CLASSy allow us to estimate the physical depths of individual structures in molecular clouds, thus providing new observational insight into the true nature of those clouds.¹⁹ Our observational result suggests that Barnard 1 is an overdense region within the larger Perseus Molecular Cloud that could have formed a sheet-like geometry from supersonic turbulence. Lee et. al (2014, submitted) present the same size–linewidth relations using CLASSy observations of Serpens Main, and find similar behavior. We will compare the results from all our CLASSy regions in an upcoming cross-comparison paper.