The Formation of a Stellar Association in the NGC 7000/IC 5070 Complex: Results from Kinematic Analysis of Stars and Gas

We searched for lists of candidate YSOs from published studies of the NAP region using the "YSO Corral" (Hillenbrand & Baliber 2015), a curated database for YSOs.

The previous studies have used a variety of techniques to identify candidate members, as listed in Table 1. These included selection based on Hα emission, infrared excess, X-ray emission, placement on color–magnitude diagrams (CMDs), and spatial clustering. Hα emission and infrared excess are mostly sensitive to disks or accretion, while X-ray emission can be used to detect pre–main-sequence stars both with and without disks. Each of these methods yields different types of contaminants and rates of contamination, so we aim to reassess memberships of their candidate members using astrometry from Gaia, as described below. Contaminants to YSO searches in the Galactic plane include both field stars and extragalactic sources; in particular, dusty (post-)asymptotic giant branch (AGB) stars and star-forming galaxies may be major sources of contaminants for selection based on infrared excess (Robitaille et al. 2008), while active galactic nuclei and foreground active M dwarfs may dominate X-ray contaminants (Broos et al. 2013).

We obtained stellar astrometry from Gaia's second data release (DR2; Gaia Collaboration et al. 2018), which provides five-parameter astrometric solutions for 1.3 billion objects to as faint as G ≈ 21 mag (Lindegren et al. 2018). These solutions provide source positions in R.A. (α) and decl. (δ), parallax (ϖ), and proper motion ⁵ (${\mu }_{{\alpha }^{\star }}$, μ_δ ). For DR2, most stars in the direction of NAP were observed during ∼15 visibility periods; typical formal uncertainties for a G = 15 mag star are ∼0.03 mas in parallax and ∼0.05 mas yr⁻¹ in proper motion. In addition to the formal uncertainty, ∼0.04 mas and ∼0.07 mas yr⁻¹ correlated systematic uncertainties affect these measurements (Lindegren et al. 2018). There have been a variety of efforts to estimate the systematic zero-point offsets (e.g., Leung & Bovy 2019, and references therein); in most cases we try to work with the observed quantities (i.e., ϖ), but when distances are required we apply a 0.0523 mas parallax correction estimated in the aforementioned paper.

We use the Gaia catalog both for cross-matching to the previously identified YSO candidates and to search for new member candidates. For cross-matching, we use a match radius of 12 and select the nearest source as a match, yielding 1939 matches out of ∼3500 candidates. This match radius is large enough that changes in position from proper motion are unlikely to affect whether a counterpart is found. The relatively low match rate is because, due to extinction, the NAP population extends to faint optical magnitudes, far below the Gaia detection limits. To estimate the rate of false matches, we artificially shifted the DR2 coordinates 2' north and ran a cross-match using these coordinates. We found 52 matches to the shifted coordinates, suggesting an incorrect match rate of up to 3%. This rate should be regarded as an upper limit on match errors, because a true counterpart, if it exists, would likely take precedence over a field star if both lie within the match radius.

We apply cuts to the Gaia data to ensure that the measurements are sufficiently precise to be useful for distinguishing between members and nonmembers. We follow the recommendations ⁶ of the Gaia Data Processing and Analysis Consortium and use only sources with renormalized unit weight error (RUWE) less than 1.4. The Gaia catalog also tabulates excess noise in the astrometric fit. Higher values of excess noise (e.g., >1 mas; Kuhn et al. 2019) could indicate some acceleration of a source (e.g., due to a binary) that might affect the astrometric solution, so we only include sources with excess noise ≤1.0 mas. And, finally, we require astrometric_sigma5d_max ≤ 0.5 mas to limit the maximum statistical measurement uncertainties on parallaxes and proper motions. These quality cuts only slightly affect the peak of the parallax distribution (shifting it by the equivalent of 4 pc), but the disadvantage of including all sources would be that the modes of the distribution become broader and source classifications become less certain. We also removed objects with ϖ > 4 mas because they are clearly in the foreground and interfere with the mixture model analysis in Section 3.

We obtained ¹³CO J = 1–0 (110.201354 GHz) observations of the molecular clouds associated with the NAP region using the Five College Radio Astronomy Observatory (FCRAO) 14 m telescope in 1998 June. The CO map covers most of the region of interest in this study.

The observations were obtained with the Second Quabbin Optical Imaging Array (SEQUOIA; Erickson et al. 1999). The array contains 16 pixel elements arranged in a 4 × 4 grid with separation of 88'' on the sky. At the time of the observations, 12 of the pixels were functional. The spectrometer contained 512 channels with a bandwidth of 40 MHz, which provided a channel spacing of 0.21 km s⁻¹. The FCRAO telescope has an FWHM angular resolution of 47'' at the ¹³CO J = 1–0 frequency. The forward scattering and spillover efficiency of the telescope at the observed frequency was η_FSS = 0.7, while the main-beam efficiency was η_mb = 0.45 (Heyer & Terebey 1998). The data were obtained in position-switching mode with spectra that were sampled every 44'' on the sky, or approximately the FWHM resolution. The typical rms noise in the spectra is ${\rm{\Delta }}{{T}_{A}}^{* }=0.22$ K per channel. The maps used for our analysis are in units of ${{T}_{R}}^{* }$, calculated from the atmosphere-corrected antenna temperature ${{T}_{A}}^{* }$ by dividing by a correction factor (Kutner & Ulich 1981). Here, we use the correction factor η_FSS since the ¹³CO emission is usually more extended than an arcminute in the NAP region.

Our ¹³CO data are similar in sensitivity and spatial resolution to a ¹³CO map published by Zhang et al. (2014) using the Purple Mountain Observatory Delingha 13.7 m telescope. They also provide ¹²CO and C¹⁸O maps.

Finally, high-dispersion optical spectra of the Bajamar Star were taken on 2019 December 1 and 2020 January 3 (UT) for the purposes of confirming the existing spectral typing that was based on low-resolution spectra and assessing the source radial velocity. These data were taken with the Keck I telescope and HIRES (Vogt et al. 1994) and cover ∼4800–9200 Å at R ≈ 37,000.

Out of the sample of YSO candidates, 395 objects are confirmed as members, while 814 objects are reclassified as nonmembers, and 2264 sources either do not have a Gaia counterpart or do not meet our quality criteria and therefore have uncertain classifications. The numbers of candidates in each category are tabulated in Table 1 for each of the published lists of YSO candidates. Table 2 provides a list of the astrometric members.

It is not surprising that nearly all historical studies have some contamination. For instance, out of the 68 objects from the pioneering study of the NAP region by Herbig (1958), Gaia provides sufficiently high quality measurements for 27 of them. Of these stars, our analysis has validated the membership of 21 objects, while rejecting membership for the other 6. The rejected members include LkHα 132, LkHα 133, LkHα 147, LkHα 183, LkHα 192, and LkHα 193, which all have parallax values that are too small to be a members of the region, suggesting that they are background Be stars.

For some categories of YSO candidate, Gaia's requirement that a source must be sufficiently bright in the optical may induce a selection bias toward higher contamination rates. For example, out of the 1286 Spitzer/MIPS sources classified as YSOs by Rebull et al. (2011) on the basis of infrared excess, only 30% met our Gaia criteria for inclusion in the analysis. Of those 30%, 131 were rejected and 251 were validated. However, this ratio is not representative of the whole sample of 1286 MIPS candidates because the YSOs with the most prominent infrared excesses tend to be very young and deeply embedded, hence not optically visible.

The spatial distribution of candidates classified as members, nonmembers, or uncertain membership is shown in Figure 3. These points are overplotted on an outline of the optical nebulosity making up the North America and Pelican components of the nebula. Nonmembers are widely dispersed throughout the whole region, while uncertain members are more tightly clustered, and validated members are the most tightly clustered. This suggests that the spatial distribution of stars is dominated by clustered groups, rather than distributed stars. Although a few members are located a degree or more away from the main groups, there is no significant nonclustered population scattered throughout the entire ∼3° diameter region.

Zoom In Zoom Out Reset image size

Our sample of 395 objects is clearly incomplete, missing both objects that were not selected in previous studies (e.g., pre–main-sequence stars without disks) and objects that could not be evaluated by Gaia. There is also evidence that our classifier has rejected a few legitimate YSOs. From a sample of 41 high-confidence YSOs studied by Findeisen et al. (2013), 4 were rejected by the classifier, likely spuriously. Altogether, it is likely that the NAP region contains at least several thousand stellar members. This assertion is based on Spitzer observations that have collectively identified ∼2000 YSO candidates (Guieu et al. 2009; Rebull et al. 2011), most of which are not included in our study because they are inaccessible to Gaia. Furthermore, in Section 8 we identify >1000 new candidate members based on Gaia astrometric data. Strictly speaking, our membership classification is for the Gaia source, so the chance alignment of a mid-infrared YSO with a Gaia field star would likely result in rejection of the source.

Both the Bajamar Star (the primary ionizing source) and HD 199579 pass our Gaia membership criteria in Section 3. However, the parallaxes and proper motions of both of these O stars place them near the edges of parameter space for objects identified as members. For HD 199579 (ϖ = 1.06 ± 0.06 mas, ${\mu }_{{{\ell }}^{\star }}$ = −1.4 ± 0.1 mas yr⁻¹, and μ_b  = −1.6 ± 0.1 mas yr⁻¹), its parallax places it behind most other cluster members. It lies on the northeast extreme of the NAP region and is moving in this direction relative to the rest of the NAP members.

The Bajamar Star (ϖ = 1.47 ± 0.08 mas, ${\mu }_{{{\ell }}^{\star }}=-3.4\,\pm 0.1$ mas yr⁻¹, and μ_b  = −2.9 ± 0.2 mas yr⁻¹) has a parallax measurement that is 2σ greater than the median parallax of the system, and the ${\mu }_{{{\ell }}^{\star }}$ proper motion is more negative than that of most other members. The Gaia data do not indicate any obvious problems with the astrometry; the RUWE is 1.04, which is within the recommended <1.4 range. The source has a small, but highly statistically significant, excess noise of 0.346 mas, which could be explained if it is a binary system. In Section 5 we argue that the star is unlikely to be in front of the complex, so the discrepant parallax measurement is likely to be the result of statistical measurement uncertainty.

From our HIRES spectrum of the Bajamar Star, we find an O4–O6 spectral type. This estimate is based on the presence of He ii λ5412 with W_λ  = 1.1 Å and the barely discernible He i λ4922 (W_λ  < 0.18 Å), He i λ5015 (W_λ  < 0.13 Å), and He i λ5047 (W_λ  < 0.08 Å).

We can also confirm the statement in Maíz Apellániz et al. (2016) that the different absorption-line species have different radial velocities. However, it is not clear based on our spectra that the source is a spectroscopic binary. While there is evidence for asymmetries in the line profiles, especially H i and He i, a high-velocity wind, as is expected from such a massive star, would produce similar features. Furthermore, a cross-correlation analysis of our HIRES spectrum does not reveal indications of two peaks, though the lines are extremely broad owing to the rapid rotation and possibly other effects, and it certainly would be possible to mask two sets of broad lines within the profiles.

Using τ Sco as a radial velocity standard, we derive v_helio = 11.6 ± 16.5 km s⁻¹ for our first observation and v_helio = −17.1 ± 16.8 km s⁻¹ for our second observation, both based on the single He ii λ5412 line. Although the errors are large, this is our best estimate of the stellar radial velocity. The velocities derived from the only other measurable and noninterstellar lines in the first spectrum are −100 km s⁻¹ from a single He i line and −90 and −35 from H i (Hα and Hβ respectively), very different from the He ii velocity. The second spectrum does show a radial velocity shift. Measured directly, rather than going through a narrow-line star, the shift in the strong He ii λ5412 line is by 28 ± 2 km s⁻¹. This certainly suggests binarity. However, an order containing He i λ5876 and C iv λλ5801, 5811 shows a velocity difference of around −40 km s⁻¹, while He ii λλ6406, 6527 suggest an 8–10 km s⁻¹ difference and both He i λ6678 and He ii λ6683 are closer to +40 km s⁻¹. In Hα and Hβ the shift between the spectra is about 15 ± 3 km s⁻¹. It remains unclear whether we should interpret the measurements as a single star having the same velocity as the molecular gas and strong wind that manifests in confusing absorption velocities, or as a binary with the more massive star producing much of the He ii absorption and both stars contributing to the relatively blueshifted metal, He i, and H i profiles.

When analyzing the stellar population of the NAP region, it is convenient to divide the stars into subclusters. As for all cluster analysis problems, multiple possible strategies can achieve this, with no one approach being "best" (Everitt et al. 2011). We find that a Gaussian mixture model with hierarchical combination of clusters, outlined below, gives reasonable results that are astronomically useful.

We use a Gaussian mixture model analysis implemented by mclust (Scrucca et al. 2016) in R to identify groups of stars in four-dimensional position–proper-motion space. The analysis is performed on the variables ℓ, b, ${\mu }_{{{\ell }}^{\star }}$, and μ_b , each of which is normalized by subtracting the mean value and dividing by the standard deviation. The program then selects the number of clusters, the free parameters of the Gaussian components, and the values of the parameters via minimization of the BIC. For modeling the clustering of stars, the mixture model approach has several benefits, which include the ability to identify clusters that differ in both compactness and number of members, the ability to assign stars from the dense cluster core as well as the less dense cluster wings to the same cluster, and the ability to deal with overlapping clusters. On the other hand, stars clusters are typically not well fit with Gaussian distributions, which can lead to multiple Gaussian components used to fit a single cluster (Kuhn & Feigelson 2019). For the NAP stars, mclust selects a model with eight ellipsoidal components with equal shape and orientation, where several of the Gaussian components are significantly overlapping or nearly concentric.

The mclust package includes the function clustCombi to deal with such cases by hierarchically merging multiple Gaussian components into single clusters using an entropy criterion (Baudry et al. 2010). For our case, we find a large change in normalized entropy going from six to seven clusters, so we use this threshold to cut the hierarchical tree at six groups. The stars assigned to each of these six groups, labeled A through F, are overplotted on the optical DSS image of the NAP region in Figure 4.

The results of the cluster analysis algorithm can be qualitatively checked for reasonableness by examination of the group assignments in Figures 4–6. In any one of the projections, the edges of some groups overlap. However, the stars that overlap in one projection tend to be different from the stars that overlap in a different projection, and the groups are mostly well separated.

In projection on the sky (Figure 4), Group A is fairly isolated in the southwest corner of the region. Groups B, C, and D form a conglomeration in the northwestern portion of the nebula, more or less corresponding to the regions known as the "Pelican's Hat," the "Pelican's Neck," and the "Atlantic," respectively. Groups E and F lie in the southeastern/eastern region—Group E is composed of several clumps of stars on the dark "Gulf of Mexico" cloud, while Group F is a single clump to the northeast of the "Gulf of Mexico" superimposed on an area with bright nebulosity.

The plot of α versus ${\mu }_{{\alpha }^{\star }}$ (Figure 5, top left) shows the clearest separation of the stars into distinct groups. In the center of this diagram, Group D—the largest group—forms a diagonal swath from lower left to upper right. This group is clearly separate from Group E to its left, which follows a parallel diagonal track. On the right, Group C is located adjacent to Group D but is much more tightly clumped in α and ${\mu }_{{\alpha }^{\star }}$. On the right side of the diagram, Group B appears partially entangled with Group C, and Group A appears partially entangled with Group B. However, both Groups A and B are looser than Group C. The spatial separation of A from the rest of the groups makes it clear that A is distinct. However, whether B could be considered part of the same group as C is less clear. Group B contains only a few stars from our sample; however, these stars coincide with the deeply embedded "Pelican's Hat" cluster of stars discovered by Spitzer (Guieu et al. 2009; Rebull et al. 2011), most of which are unseen by Gaia.

In the plot of δ versus μ_δ (Figure 5, top right), the diagonal swath formed by Group D can still be seen, but the other groups appear distributed in a more vertical direction rather than diagonal. This plot enhances the separation between Groups A, B, and C. Groups C and D are less cleanly separated, but the stars in Group C do not continue the same diagonal trend as Group D. The bottom panels in Figure 5 show the other combinations of position and proper motion. Finally, on the plot of ${\mu }_{{\alpha }^{\star }}$ versus μ_δ (Figure 6), Groups D and E are in the center, Group A is at the top extreme, Groups B and C are at the bottom extreme, and Group E is at the rightmost extreme.

Intriguingly, the Bajamar Star is assigned to Group D even though it is spatially closer to stars from Group E in the "Gulf of Mexico." This is a result of the proper motion of the Bajamar Star, which is rapidly moving toward Group E away from the center of Group D. This proper motion does not match the proper motions of nearby stars in Group E; however, it does fit with the proper-motion gradient seen for stars in Group D (Figure 5).

Our division of stars into kinematic groups is meant to be useful for understanding the structure of the NAP region, but it is not the full description of the cluster structure. For example, nearby groups may have formed from related star formation events at different locations on the same cloud, and thus, from an astrophysical perspective, it is ambiguous whether these should be considered to be the same group. Furthermore, if the region were nearer (e.g., at the distance of Taurus; Luhman 2018; Fleming et al. 2019), we might be able to detect finer velocity differences with Gaia that could be used to subdivide the groups further, but if the region were more distant (e.g., at the distance of the Carina Nebula; Kuhn et al. 2019), the measurement uncertainties in the proper motion would be poorer and, thus, it might not have been possible to distinguish the groups.

There is spatial correspondence between the stellar groups and the clouds (Figure 10). Groups A and F are associated with the minor clouds 1 and 13, respectively. Group E is embedded in the "Gulf of Mexico," which consists of a major cloud component (Cloud 9) and several minor cloud components (Clouds 6, 10, 11, and 13) in the ¹³CO map. Group C lies on a bright rim cloud at the southeastern edge of Cloud 2. The stars in Group B are distributed in the northern part of Cloud 2, as well as Cloud 3.

Group D, the most spatially extended stellar group, has the most complicated relationship to the ¹³CO clouds. This group is projected over a region that spans several clouds. However, Cloud 8 is the dominant cloud in this part of the NAP region and lies near the center of Group D. The stars selected in our Gaia analysis do not lie on top of Cloud 8, but instead the densest parts of Group D are next to the cloud. This indicates either that Group D is partially embedded in Cloud 8, which would imply that Cloud 8 is the remainder of the natal cloud that formed this group, or that Group D lies behind the cloud. Stars from this group are also superimposed on other nearby minor clouds (Clouds 5, 6, and 7), which may either be smaller remnants of the natal cloud or chance superpositions. The edges of Group D also intersect the edges of Cloud 2, Cloud 3, and even (in the case of the Bajamar Star) Cloud 9.

If the stellar groups and clouds are associated, then it is likely that they have similar mean motions. This enables us to associate mean proper motions with clouds and radial velocities with stars, creating a three-dimensional picture of the velocities in this region.

Line profiles for each of the clouds (Figure 11) were constructed by spatially integrating the ¹³CO gas detected within the contours on Figure 10. Most of these profiles display several distinct peaks. In Table 4, we report basic statistics of the velocity distributions, including the first (mean) and square root of the second (standard deviation) moments, the mode, and the FWHM of the tallest mode. The mode and FWHM (Δv_FWHM) may be more characteristic of the kinematics of the clouds because these statistics are less susceptible to merging components with very different velocities that may be from distinct clouds. Thus, we base our discussion of cloud dynamics (Section 6.3) on the latter quantities. We also estimate the characteristic size of each cloud, defined as the projected radius that contains half the integrated emission.

Zoom In Zoom Out Reset image size

Optical depth can affect line profiles, flattening the peaks and effectively broadening the lines (Hacar et al. 2016). To more accurately resolve the velocity structure, we integrated over the optical depths given by Equation (2) rather than ${{T}_{R}}^{* }$ (see Goldsmith & Langer 1999). This results in line profiles that are more sharply peaked with FWHM that are ∼20% narrower.

The interpretation of the multimodal structure becomes clearer with the position–velocity diagrams in Figure 12. Four of the major clouds (Clouds 2, 3, 8 and 9) are shown as examples, while the diagrams for the other clouds are accessible via the online figure set. Even some of the minor clouds (e.g., Clouds 5, 6, and 7) have complicated velocity structures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (13 panels)

Tarball of figure set images

In Cloud 2, v_lsr increases by several kilometers per second toward the southern end of the cloud, where the optically bright rim and stellar group are located. It is plausible that this could be an effect of the expansion of the H ii region pushing this end of the cloud away.

For Clouds 3 and 8, multiple components are visible, corresponding to the multiple peaks in the line profiles. In Cloud 3, there are two velocity components separated by ∼5 km s⁻¹, while in Cloud 8 there are three components, also separated by ∼5 km s⁻¹ from each other. The distinct components in Cloud 3 are also mostly separated along a north–south axis, possibly indicating that this cloud is composed of multiple subclouds. However, for Cloud 8 the main component is the dense northern end of the filament at v_lsr ≈ −5 km s⁻¹ (Figure 9), while the other components are related to overlapping filamentary structures at higher v_lsr. In our analysis, we tentatively associate the stars in Group D with the −5 km s⁻¹ gas because this is the dominant component around which the stars are conglomerated, but our data are insufficient to make a definitive conclusion (Appendix C).

The position–velocity diagram for Cloud 9 in Figure 12 shows that it has a clumpy structure, with clumps spanning the velocity range from ∼0 to 7 km s⁻¹. A small clump at −5 km s⁻¹ is probably not part of this cloud, but is instead a continuation of the long north–south filament seen in Figure 9 (top panels). Overall, Cloud 9 has the largest positive v_lsr of the complex, which, in combination with our evidence that the "Gulf of Mexico" and Group E are in front of the rest of the NAP region, suggests that Cloud 9 is infalling. However, the part of the cloud brightest in ¹³CO emission also has a fairly modest velocity of ∼1 km s⁻¹. This ¹³CO bright patch is fairly near the Bajamar Star and shows up as having a higher temperature in the ¹²CO map from Zhang et al. (2014).

From the kinematic measurements of the cloud in Table 4, several dynamical quantities can be calculated that are useful for understanding star formation in this region. For these calculations, we use the line widths for the most prominent peaks in Figure 12 to avoid combining unrelated velocity components. The dynamical quantities are given in Table 5.

Table 5. Cloud Dynamical Properties

Cloud	$\mathrm{log}\bar{\rho }$	τcross	τff	${ \mathcal M }$	Mdyn	αBM92
	(g cm−3)	(Myr)	(Myr)		(103 M⊙)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
1	−20.40	1.3	1.1	3.6	0.8	2.5
2	−20.17	1.6	0.8	4.7	2.2	1.1
3	−20.10	1.6	0.8	3.6	1.0	0.9
4	−20.01	1.8	0.7	2.5	0.3	0.6
5	−20.80	1.1	1.7	4.7	1.6	9.3
6	−19.98	0.7	0.7	4.4	0.8	3.5
7	−20.24	0.9	0.9	3.6	0.5	4.3
8	−20.65	1.6	1.4	6.1	4.7	3.2
9	−20.52	1.5	1.2	10.1	20.6	2.6
10	−19.67	0.6	0.5	2.8	0.2	2.2
11	−19.01	0.2	0.2	4.7	0.3	5.4
12	−18.38	0.1	0.1	5.3	0.2	6.4
13	−20.26	1.9	0.9	2.1	0.2	0.9

Note. Column (1): cloud designation. Column (2): mean density of the cloud, defined as $\bar{\rho }={M}_{{}^{13}\mathrm{CO}}/(4/3)\pi {r}_{\mathrm{cloud}}^{3}$. Column (3): crossing timescale, defined as τ_cross = r_cloud/Δv. Column (4): freefall timescale defined using $\bar{\rho }$ from Column (3). Column (5): Mach number. Column (6): dynamical mass of the cloud assuming α = 1 (Equation (9)). Column (7): Bertoldi & McKee (1992) estimate of the virial parameter, using ${M}_{{}^{13}\mathrm{CO}}$ for the mass, r_cloud for the radius, and σ_nt for the velocity dispersion.

Download table as:  ASCIITypeset image

If we assume that the kinetic temperature T_kin of the gas equals the excitation temperature measured for ¹²CO (∼14 K) and that the mean molecular mass is $\bar{m}=2.33$ amu, the one-dimensional sound speed in the molecular clouds would be c_s  ≈ 0.22 km s⁻¹. For each of the molecular clouds delimited in our analysis, the line width of the tallest mode is larger than can be accounted for by purely thermal broadening. The nonthermal component is

$\begin{eqnarray}&&{\sigma }_{\mathrm{nt}}=\sqrt{{\left(\displaystyle \frac{{\rm{\Delta }}{v}_{\mathrm{FWHM}}}{2.35}\right)}^{2}-\displaystyle \frac{{{kT}}_{\mathrm{kin}}}{{m}_{\mathrm{obs}}}},\end{eqnarray} \tag{ 4 }$

where m_obs = 29 amu is the mass of the observed molecule, ¹³CO. The Mach number is defined to be ${ \mathcal M }={\sigma }_{\mathrm{nt}}/{c}_{s}$. For most of the individual clouds, Mach numbers span $2\lesssim { \mathcal M }\lesssim 5$. Exceptions are Clouds 8 and 9—the most massive clouds—which have larger velocity dispersions of ${ \mathcal M }\approx 6$ and 10, respectively. The lowest velocity dispersion is found for Cloud 13 with ${ \mathcal M }\approx 2;$ this also happens to be the cloud with the lowest surface density. This range of Mach numbers is fairly typical for star-forming clouds, such as the Orion B cloud (${ \mathcal M }\sim 6;$ Orkisz et al. 2017) or Taurus (${ \mathcal M }$ < 2–3; Hacar et al. 2016). The position–velocity plots (Figure 12) show that some of this nonthermal broadening in the NAP clouds is due to cloud-scale velocity structure. In Cloud 8, the −5 km s⁻¹ component appears curved in position–velocity space, while, in Cloud 9, the large velocity dispersion is composed of multiple clumps with different velocities.

Several dynamically important quantities that can be estimated for each cloud include the mean density $\bar{\rho }$, freefall timescale τ_ff, and crossing timescale τ_cross defined as

$\begin{eqnarray}&&\bar{\rho }={M}_{{}^{13}\mathrm{CO}}/(4/3\pi ){r}_{\mathrm{cloud}}^{2}\end{eqnarray} \tag{ 5 }$

$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}=\sqrt{3\pi /32G\bar{\rho }}\end{eqnarray} \tag{ 6 }$

$\begin{eqnarray}&&{\tau }_{\mathrm{cross}}={r}_{\mathrm{cloud}}/{\sigma }_{\mathrm{nt}}.\end{eqnarray} \tag{ 7 }$

Both the crossing timescales and the freefall timescales range from 0.1 to 2 Myr. For the main clouds with ongoing star formation, the typical freefall timescale is τ_ff ∼ 1 Myr. We note that these timescales are similar to the ∼1 Myr ages inferred for the pre–main-sequence stars.

The velocity dispersions can be used to estimate cloud dynamical masses. However, these calculations depend on the virial ratio of the clouds, $\alpha =2{ \mathcal T }/| { \mathcal W }|$, where ${ \mathcal T }$ is kinetic energy and ${ \mathcal W }$ is potential energy. Previous studies typically use the definition from Bertoldi & McKee (1992) that assumes a spherical, uniform density cloud, giving

$\begin{eqnarray}&&{\alpha }_{\mathrm{BM}92}=\displaystyle \frac{5{\sigma }_{\mathrm{nt}}^{2}{r}_{\mathrm{cloud}}}{{{GM}}_{\mathrm{cloud}}}.\end{eqnarray} \tag{ 8 }$

However, the true value of α depends on the three-dimensional distribution of mass within the cloud (Singh et al. 2019; Guszejnov et al. 2020), with corrections for centrally concentrated mass distributions and filamentary structure typically being a factor of several (Bertoldi & McKee 1992). If we let α_BM92 = 1, we obtain a dynamical mass estimate

$\begin{eqnarray}&&{M}_{\mathrm{dyn}}=5\,{r}_{\mathrm{cloud}}{\sigma }_{\mathrm{nt}}^{2}/G.\end{eqnarray} \tag{ 9 }$

In Figure 13, we plot $\mathrm{log}{M}_{\mathrm{dyn}}$ versus $\mathrm{log}{M}_{{}^{13}\mathrm{CO}}$ for the 13 clouds we have identified in the NAP region. These quantities are strongly correlated (p < 0.001 from Kendall's τ test), with an approximately proportional relationship over a dynamic range of ∼200 in cloud mass. However, the dynamical masses are systematically ∼0.4 dex higher, with ∼0.4 dex scatter in the relation. The similarity in masses derived by the two methods to within a factor of ∼3 helps corroborate our ¹³CO-based mass estimates in Table 4. Vázquez-Semadeni et al. (2019) point out that, observationally, it is nearly impossible to distinguish between the freefall and virial velocities, which only differ by a factor of 2 in mass or a factor of $\sqrt{2}$ in velocity. Table 5 provides both M_dyn, calculated from the velocity dispersions using Equation (9), and α_BM92, calculated from both the velocity dispersions and the integrated ¹³CO cloud masses using Equation (8).

Zoom In Zoom Out Reset image size

In Figure 10, we indicate the motions of NAP members using arrows, which, as in the previous figures, are color-coded by stellar group. The figure shows that stars near each other tend to have similar velocities, but the different groups appear to have fairly random motions relative to one another.

Figure 14 (top left) shows the relative 3D motions of the group centers inferred from both the mean proper motions of stars and the radial velocities of the associated clouds. This diagram shows no clear pattern of either convergence or divergence. Stars in Group A tend to be moving mostly north and slightly west—in a direction tangential to the NAP system center. Stars in both Groups B and C are moving southeast, inward toward the center of the system. Group E is drifting slowly east, away from the center of the NAP complex in the xy plane, but it is likely moving rapidly toward the system center in the third dimension. Stars in Group F are moving west—toward the center of the region. Finally, the center of Group D appears to be nearly stationary in the xy plane in this reference frame. However, most of the individual stars in this group are directed away from its own center as can be seen in Figure 10.

Zoom In Zoom Out Reset image size

The Bajamar Star is not stationary relative to the center-of-mass reference frame, but is instead rapidly traveling southeast away from Group D and toward Group E. The other O star, HD 199579, is traveling northeast away from the NAP complex. In this global reference frame, ⁹ both stars have velocities >6 km s⁻¹. Such speeds would classify both objects as "walkaways" (Renzo et al. 2019; Schoettler et al. 2020), although neither object has traveled outside the star-forming region.

In Kuhn et al. (2019), we presented several methods for characterizing stellar kinematics and testing for expansion, contraction, or rotation of a stellar system. Here, we apply these methods to each of the stellar groups.

To characterize total velocity dispersion in a group, we fit the (v_x , v_y ) distribution with a bivariate Gaussian, taking into account measurement uncertainties, which have the effect of artificially broadening the velocity dispersion. ¹⁰ The stellar velocities and the Gaussian fit are shown in Figure 15 for each of the six groups. We summarize these distributions using the quantity σ_1D (Table 6), defined as the square root of half the trace of the Gaussian's covariance matrix (Kuhn et al. 2019, their Equation (17)). For the six groups, σ_1D ranges from 1 to 2 km s⁻¹, being smallest for Groups B and F and largest for D and E. Group F is a small group in a small cloud, so it is unsurprising that this group has a small velocity dispersion. In contrast, Groups D and E are large groups, located in a massive cloud (E) or in a partially dispersed cloud (D), so it is unsurprising that they have larger velocity dispersions. However, it is notable that Group C, spatially adjacent to Group D, has a much lower velocity dispersion. In many cases, the velocity dispersions show signs of being anisotropic, similar to the results for other star-forming regions (Kuhn et al. 2019).

Zoom In Zoom Out Reset image size

We use several plots that are analogous to those in Kuhn et al. (2019), to investigate evidence for expansion. Figure 16 shows Group D, the largest group in our sample where expansion is most evident. In this group, the fastest-moving stars, with velocities >5 km s⁻¹, are located near the edges of the region, and they are generally directed away from the center of the group. The top right panel uses both arrow direction and hue to indicate the direction of motion, with color saturation indicating statistical uncertainty on direction. This use of color allows patterns of motions to be identifiable as a gradient in the color of the marks, in this case showing that stars at different positions in the group are oriented in different directions, all preferentially away from the center.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (14 panels)

Tarball of figure set images

The two bottom panels show statistical tests for expansion. The bottom left panel shows the distribution of v_out values; if a system is expanding, then these will be predominantly positive. For Group D, the weighted median v_out value is 1.9 ± 0.4 km s⁻¹, more than 3 standard deviations above 0, clearly demonstrating that the group is expanding. This expansion velocity is larger than any of those measured for the 28 systems investigated by Kuhn et al. (2019). The bottom right panel shows the mean expansion velocities for bins at different radii from the group center. Expansion velocity appears to increase with distance from the center, a trend that is seen in many of the expanding young star clusters and associations from Kuhn et al. (2019).

The other groups contain fewer stars than Group D, and the clear evidence for expansion is not found using these diagrams. For the other groups, arrow diagrams like those in the top two panels of Figure 16 are included in the figure set, but there are too few stars to make the other plots. For these groups, the weighted median v_out values are not significantly different from 0 (Table 6) owing to the large statistical scatter. Nevertheless, median v_out may not be the most sensitive test for expansion, given that this statistic combines multiple dimensions and does not take advantage of the tendency, seen in many expanding regions, for velocity to increase with distance from the center.

We apply a further test for expansion, using the nonparametric Kendall's τ test to test for correlation between x and v_x and between y and v_y (Table 6). Statistically significant results would indicate a velocity gradient, which would imply expansion if the gradient is positive and contraction if negative. We find strong statistical evidence for expansion in both dimensions for Group D, moderate statistical evidence for expansion in both dimensions for Group B, and strong statistical evidence for expansion in only one dimension (x) for Groups C and E.

Our results, along with the results from previous Gaia studies, suggest that stars in expanding young stellar groups often have velocities that are proportional to the distance from the groups' centers (Kuhn et al. 2019; Román-Zúñiga et al. 2019; Wright et al. 2019; Zamora-Avilés et al. 2019; Melnik & Dambis 2020). However, Wright et al. (2019) point out that the expansion may be anisotropic. In Appendix D, we describe a linear model for velocity as a function of position, using a method that can account for anisotropy and is independent of choice of coordinate system. Table 6 reports several of the interesting properties from the models, including velocity gradients in the directions of maximum and minimum expansion, the position angle of anisotropy, and the intrinsic scatter that remains after the velocity gradient is accounted for.

We find statistically significant expansion for Groups C, D, and E, but no group had statistically significant rotation. Group D has the strongest relation between position and proper motion, with expansion in all directions, but with a larger gradient in the east–west direction (0.5 km s⁻¹ pc⁻¹) than in the north–south direction (0.3 km s⁻¹ pc⁻¹). In contrast, Groups C and E had nonzero expansion gradients only in the east–west direction, not the north–south direction. We note that the anisotropy position angles of all three groups are similar. The regression analysis and the calculation of the model parameters are described in greater detail in the appendix.

Figure 17 shows scatter plots of position versus velocity with the model regression lines overplotted. The graphs of v_x versus x (top left) and v_y versus y (top right) exhibit positive correlations if there is expansion, while v_y versus x (bottom left) and v_x versus y (bottom right) exhibit correlations if there is rotation. The regression line and 95% credible interval from the linear model are shown in magenta. The intrinsic velocity scatter, applied equally to all points, is indicated by the red error bar (1σ), while the individual 1σ measurement uncertainties are indicated by the gray lines. For Group D, both the v_x versus x and v_y versus y plots show strong correlations, while for Groups C and E correlations are only strong on the v_x versus x plots owing to their θ ≈ 90° anisotropy position angles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (3 panels)

Tarball of figure set images

The O stars are both among the fastest-moving members of Group D. The Bajamar Star has a velocity of 6.6 ± 0.5 km s⁻¹ relative to the center of the group, while HD 199579 has a velocity of 6.4 ± 0.4 km s⁻¹. However, neither of these velocities is discrepant from the linear model. This suggests that whatever phenomenon accelerated these stars to their current velocities is also responsible for the expansion of the group as a whole. If we make the simplifying assumption that the velocities have been approximately constant, the Bajamar Star would have been nearest, in projection, to the center of Group D ∼1.5 Myr ago, while HD 199579 would have been nearest ∼1.8 Myr ago.

To ensure that our measurements of velocity gradients are robust, we examine how these would be affected by the systematic correlated errors in astrometry (Lindegren et al. 2018, their Section 14). They found that the covariance of quasar proper motion as a function of angular separation decreases from ∼4000 μas² yr⁻² at 007 separation to a local minimum of ∼40 μas² yr⁻² at 043 separation. This could induce an artificial proper-motion gradient of ∼0.17 mas yr⁻¹ deg⁻¹, which would look like a velocity gradient of 0.045 km s⁻¹ pc⁻¹ in a stellar association at a distance of 795 pc. Given that the velocity gradients that we detect are significantly larger than this value, they are most likely astrophysical.

Several mechanisms, including disruption by tidal forces, cloud dispersal, and sequential star formation, may play roles in the expansion of stellar groups in the NAP complex.

A clumpy mass distribution, as observed here, will generate internal tidal forces that affect stellar dynamics and influence the potential assembly of either bound clusters or unbound associations. Tidal fields would naturally produce an anisotropic velocity gradient similar to those observed in the NAP region. Furthermore, tidal forces would preserve clumpy substructure within a group (e.g., D and E) as it is stretched, while processes driven by dynamical relaxation would be more likely to erase these substructures. Tidal forces from giant molecular clouds have long been known to be a major cause of star cluster disruption (Spitzer 1958; Gieles et al. 2006), and these effects can occur in the star-forming environment itself, inhibiting the initial formation of a bound system (Kruijssen et al. 2011, 2012). Zamora-Avilés et al. (2019) have run numerical simulations in which gas clouds, driven away from the newly formed young stellar cluster by feedback, exert tidal forces on the cluster that pull it apart.

The tidal acceleration ${a}_{t,\mathrm{axial}}$ per unit of displacement Δr along the separation axis due to a mass M at distance r is

$\begin{eqnarray}&&{a}_{t,\mathrm{axial}}/{\rm{\Delta }}r=2\,G\displaystyle \frac{M}{{r}^{3}}.\end{eqnarray} \tag{ 12 }$

For Group D, the magnitude of the tidal acceleration from Cloud 2 (2.1 × 10³ M_⊙) located ∼5 pc to the west of the center of the group (assuming that the displacement in the third dimension is small) is 0.15 km s⁻¹ pc⁻¹ Myr⁻¹. This acceleration would be enough to induce the observed ∼0.5 km s⁻¹ pc⁻¹ east–west velocity gradient in ∼3 Myr. Other nearby clouds, notably Clouds 3, 4, and 8 (a potential remnant of Group D's natal cloud), could also contribute to the tidal forces on this group. These clouds would also exert tidal forces on the other nearby, expanding group, C, located at the edge of Cloud 2. However, precisely calculating the effect of tidal forces is challenging owing to uncertainties in the third dimension, velocity differences between the groups and clouds, and systematic uncertainties in cloud mass.

Expansion of Group E is more difficult to explain with tides, since its separation from other parts of the NAP region by ∼35 pc along the line of sight means that it would be much less strongly affected by tidal forces from other clouds. Thus, any tidal force would need to be generated by clumpy distributions of mass within Cloud 9 itself.

The ongoing dispersal of the molecular clouds due to O star winds and photoevaporation, combined with outflows from low-mass stars (Bally & Reipurth 2003; Bally et al. 2014), will weaken the binding energy of the stellar groups, causing them to expand and/or disperse (e.g., Tutukov 1978; Adams 2000; Kroupa et al. 2001, and many others). In the northwest, the shell-like structure may be the outer rim of a bubble blown by the Bajamar Star, which would have been closer to this part of the cloud ∼1.5 Myr ago. Cloud dispersal in this region may have affected the stellar groups that lie near the shell, particularly Groups C and D. In contrast, photoevaporation does not appear to be as significant in the Gulf of Mexico region that contains Group E. Nevertheless, Bally et al. (2014) found dozens of outflows here that are reprocessing a significant volume of the Gulf of Mexico cloud through supersonic shocks. Such outflows may dominate the injection of energy and momentum in the absence of massive stars (Li & Nakamura 2006; Nakamura & Li 2011) and may contribute to the cloud's mass loss (Arce et al. 2007).

We examine the scenario in which Group D formed as an embedded cluster in virial equilibrium in a molecular cloud that was subsequently dispersed. Given Group D's measured velocity dispersion of σ_1D = 2 km s⁻¹ and assuming an initial size of ∼1 pc, the corresponding virial mass of the system would have been ∼9 × 10³ M_⊙. Thus, the natal cloud mass would have been ∼6 times greater than the estimated mass of Cloud 8 and ∼1.5 times the sum of the masses of all the clouds in the northwestern part of the NAP complex. Wright & Parker (2019) have demonstrated that expansion in this scenario could be anisotropic if it were preceded by collapse and bounce in an asymmetric gravitational potential, setting up an anisotropic velocity dispersion.

A third explanation for expanding groups is that they are produced by a second generation of star formation in an expanding shell of material. In regions where star formation has been attributed to material swept up in shells around expanding H ii regions (e.g., Patel et al. 1998) the stars formed in this way would presumably inherit the velocities of the material in the expanding shells. In the NAP region, there is ongoing star formation in several of the clouds making up the shell structure, so it is plausible that such a scenario could be happening here too. In particular, Group C is associated with the bright rim cloud at the edge of Cloud 2. For several bright rim clouds in other star-forming regions, Getman et al. (2007, 2012) have shown that the stars in front of the cloud exhibit an age gradient, with the youngest stars nearest the cloud. In such a scenario, as a cloud is driven outward by an expanding shell, the sequential formation of stars could yield a velocity gradient.

Given the conditions in the NAP region, it seems likely that multiple scenarios operate simultaneously to induce the observed expansion patterns. For example, the lessening of the binding energy due to cloud dispersal would make groups more susceptible to disruption by interactions with neighboring clouds. Although theoretical work has shown that binary-star dynamics can preferentially eject O stars at high velocities (>30 km s⁻¹; Oh & Kroupa 2016), this mechanism does not appear to be required here because the O stars and low-mass stars follow the same velocity pattern.

The Gaia data show that individual groups are expanding, but the relative motions of the groups exhibit no clear pattern of either convergence or divergence (Figure 14, top left). This is similar to the situation found by Kuhn et al. (2019) in other large star-forming complexes like NGC 2264, NGC 6357, or the Carina Nebula and may help explain why expansion has been difficult to detect in young stellar associations analyzed as a whole (e.g., Wright et al. 2016; Ward & Kruijssen 2018; Wright & Mamajek 2018). The velocity dispersion of all stars in the NAP complex is σ_1D = 2.5 km s⁻¹ (Figure 14, bottom panels). This, combined with the separations of several to tens of parsecs between groups, means that the mass needed to gravitationally bind the system (>10⁵ M_⊙) is orders of magnitude greater than the mass in stars. Thus, it is impossible for the system to coalesce as a bound cluster. Nevertheless, it may be possible in some groups for bound remnants to remain behind even after most stars have escaped (Baumgardt & Kroupa 2007).

We can use our inferences about the dynamical state of the stars and gas, as well as the constraints on stellar ages, to test various theoretical scenarios for star formation as applied to the NAP region. These scenarios broadly break down into fast star formation, on the timescale of the freefall collapse of the molecular clouds (e.g., Elmegreen 2000; Hartmann et al. 2012), and slow star formation that persists over multiple freefall timescales (e.g., Tan et al. 2006; Krumholz & McKee 2020). In the former case, molecular clouds would have little pressure support (Vázquez-Semadeni et al. 2019), while in the latter case turbulent pressure support stabilizes the clouds against collapse and allows for an extended duration of star formation (Padoan & Nordlund 2011).

In the NAP region, we demonstrated that the velocity dispersion within individual clouds is consistent with the velocity dispersion that would be expected from either freefall collapse or virial equilibrium given the systematic uncertainties on the cloud masses (Section 6.3). When we consider the whole region (Figure 14), the combined velocity dispersion of all gas is larger than that of the individual clouds, with σ_all = 4.2 km s⁻¹. In the complex, half the gas is contained within a projected radius of 9.4 pc, so the dynamical mass for the whole system calculated using Equation (9) would be 1.9 × 10⁵ M_⊙, or about two and a half times the 7.2 × 10⁴ M_⊙ mass estimated from ¹³CO column density. Estimation depends on the virial parameter α, which we take to be 1 (virial equilibrium); however, a cloud with α = 2 (freefall) would yield a dynamical mass estimate half as large, making the differences smaller. As earlier, the systematic uncertainties make distinguishing between these scenarios difficult. Nevertheless, in either case, the velocity differences between the clouds, as well as the velocity dispersions within the clouds, can be accounted for by gravitation.

The freefall timescales for the clouds are ∼1 Myr (Table 5). This is similar to the ∼1 Myr ages of the stars (Figure 2). This implies that, in each cloud, most star formation occurred within the last 1–2 × τ_ff. There are two considerations that should be taken into account when interpreting freefall timescales. First, freefall timescales are typically calculated using the mean density of a cloud as we have done; a nonuniform density could change the mean τ_ff, but the correction factor is expected to be of order unity (Krumholz & McKee 2020). Second, we are calculating the present-day freefall timescales; however, in the past, when the stars observed today were forming, the freefall timescale would likely have been longer (Vázquez-Semadeni et al. 2019). Star cluster formation in a freefall time is expected for either clouds without turbulent support (Klessen et al. 1998) or clouds with decaying turbulence (Klessen 2000). Thus, the observation that most of the stars are only one to several τ_ff old is consistent with the rapid star formation scenarios.

The spatial clustering of star formation also favors a scenario of rapid star formation in a turbulent molecular cloud. Overall, the NAP cloud complex has features expected for a giant molecular cloud sculpted by turbulence, including a clumpy, hierarchical density structure and velocity structure that yields a larger velocity dispersion for the whole system than within the individual clouds, broadly consistent with the Larson (1981) relation. Given that the major stellar groups follow the present-day cloud structure, this supports the view that the stars froze out on a timescale less than the crossing time for the complex (Elmegreen 2000).

The similarity in age (∼1 Myr) of widely separated stellar groups raises the question of how star formation in the different clouds of the NAP complex became synchronized. Although many of the surveys preferentially selected stars that are young enough to still possess disks, even the X-ray-selected sample in the NAP region appears to follow the same age distribution (Appendix B). The sound crossing time for the complex is t_s  = 9.4 pc/c_s  = 40 Myr, meaning that the clouds would not be in thermal contact with each other on the collapse timescale, and we would not expect them to have similarly aged populations. A similar problem has been posed by Preibisch & Zinnecker (1999) for the Upper Scorpius OB association and by Herczeg et al. (2019) for the Serpens molecular clouds. A possible solution to this is an accelerating star formation rate, which is a feature of several models, including the conveyor belt model (Longmore et al. 2014) and global hierarchical collapse (Vázquez-Semadeni et al. 2017, 2019). With sufficiently high acceleration, most of the stars in star-forming clouds would have formed very recently, even if the first stars to form in those clouds formed at different times.

Although the NAP clouds appear to be forming stars on freefall timescales, fast star formation scenarios face challenges if applied to the entire Galaxy because they would imply a Galactic star formation rate much higher than observed (Zuckerman & Evans 1974; Krumholz et al. 2006). To overcome this, theoretical fast star formation scenarios invoke the quick disruption of the molecular clouds by stellar feedback to keep star formation inefficient and regulate the Galactic star formation rate (Chevance et al. 2020). In the NAP region, where the disruptive influences of the Bajamar Star and outflows from low-mass stars are already apparent, it appears likely that feedback processes could bring star formation to an end before a substantial fraction of the cloud mass is converted into stars.