The Formation of a Stellar Association in the NGC 7000/IC 5070 Complex: Results from Kinematic Analysis of Stars and Gas
Michael A. Kuhn, Lynne A. Hillenbrand, John M. Carpenter, Angel Rodrigo Avelar Menendez
DDraft version July 8, 2020
Typeset using L A TEX twocolumn style in AASTeX63
The Formation of a Stellar Association in the NGC 7000/IC 5070 Complex:Results from Kinematic Analysis of Stars and Gas
Michael A. Kuhn, Lynne A. Hillenbrand, John M. Carpenter, and Angel Rodrigo Avelar Menendez Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA Joint ALMA Observatory, Alonso de C´ordova 3107 Vitacura, Santiago, Chile (Received May 5, 2020; Revised June 15, 2020; Accepted June 24, 2020)
Submitted to ApJABSTRACTWe examine the clustering and kinematics of young stellar objects (YSOs) in the North Amer-ica/Pelican Nebulae, as revealed by
Gaia astrometry, in relation to the structure and motions of themolecular gas, as indicated in molecular line maps. The
Gaia parallaxes and proper motions allowus to significantly refine previously published lists of YSOs, demonstrating that many of the objectspreviously thought to form a distributed population turn out to be non-members. The members aresubdivided into at least 6 spatio-kinematic groups, each of which is associated with its own molecu-lar cloud component or components. Three of the groups are expanding, with velocity gradients of0.3–0.5 km s − pc − , up to maximum velocities of ∼ − away from the groups’ centers. Thetwo known O-type stars associated with the region, 2MASS J20555125+4352246 and HD 199579, arerapidly escaping one of these groups, following the same position–velocity relation as the low-massstars. We calculate that a combination of gas expulsion and tidal forces from the clumpy distributionof molecular gas could impart the observed velocity gradients within the groups. However, on a globalscale, the relative motions of the groups do not appear either divergent or convergent. The velocitydispersion of the whole system is consistent with the kinetic energy gained due to gravitational collapseof the complex. Most of the stellar population has ages similar to the free-fall timescales for the natalclouds. Thus, we suggest the nearly free-fall collapse of a turbulent molecular cloud as the most likelyscenario for star formation in this complex. INTRODUCTIONThe way that stars and gas are distributed in a star-forming region can provide useful constraints on the con-ditions in which the stars were formed (e.g., Elmegreen2002; Parker et al. 2014; Gouliermis 2018; Feigelson2018). For the nearest star-forming regions, astromet-ric measurements by ESA’s
Gaia spacecraft (Gaia Col-laboration et al. 2016) can provide a multi-dimensionalpicture of how the young stars are clustered. Spatialand kinematic clustering has already been examined inseveral of the major star-forming regions within 1 kpc,including Orion (Kounkel et al. 2018; Großschedl et al.2018; Getman et al. 2019), Taurus (Luhman 2018; Flem-ing et al. 2019; Galli et al. 2019), ρ Oph (C´anovas et al.2019), and Serpens (Herczeg et al. 2019). We aim to doa similar analysis for North America and Pelican Nebu-lae (hereafter NAP). Given the correlation between starsand gas within many star-forming regions (e.g., Guter-
Corresponding author: Michael A. [email protected] muth et al. 2011; Tobin et al. 2009; Lada et al. 2013),better understanding of the clustering of the stars canhelp decipher the complex velocity structures seen inradio molecular-line maps of star-forming clouds (e.g.,Larson 1981; Heyer & Dame 2015).The NAP region contains a molecular cloud complex(Bally & Scoville 1980), an H ii region (W80; Westerhout1958), and a population of YSOs, many of which havebeen extensively studied (Reipurth & Schneider 2008).This complex is fairly extended, with a diameter of ∼ ◦ ( ≈
40 pc) encompassing several sites of active star for-mation. The North America (NGC 7000) and Pelican(IC 5070) nebulae make up the east/west componentsof the H ii region, which is bisected in projection on thesky by a dark lane known as L933/935 (Lynds 1962),giving rise to the characteristic shape of the nebulaeas seen in optical light. The YSOs are predominantlylocated in the dark lane (Bally et al. 2014), with thesoutheastern portion known as the “Gulf of Mexico” andthe northern part of the cloud containing the “Atlantic a r X i v : . [ a s t r o - ph . GA ] J u l Kuhn et al.
Ocean” and the “Pelican” regions. Many of the first-identified NAP YSOs were optically visible emission-linestars (e.g., Herbig 1958). More recently,
Spitzer has un-covered thousands of embedded stars and protostars inthe “Gulf of Mexico,” “Pelican,” and “Pelican’s Hat”(Guieu et al. 2009; Rebull et al. 2011). Cambr´esy et al.(2002) identified several dense clusters of stars in theNAP region based on near-infrared star counts. Dozensof outflows from young stellar objects (YSOs) attest toongoing star formation throughout the complex (Bally& Reipurth 2003; Bally et al. 2014).Another intriguing aspect of this star-forming re-gion is its primary ionizing source, 2MASS J20555125+4352246 (dubbed the Bajamar Star; Ma´ız Apell´anizet al. 2016). This source was discovered by Comer´on& Pasquali (2005) lying to the northwest of the “Gulfof Mexico” behind ∼ Another O star,HD 199579, is projected on the North America Nebula.This star, with spectral type O6.5 V ((f))z (Sota et al.2011), has generally been regarded as being too far fromthe center of the H ii region to be a main ionizing source(e.g., Herbig 1958). Whether HD 199579 is a member ofthe NAP association has been uncertain (Wendker et al.1983).The NAP complex is in the Orion-Cygnus Arm of theGalaxy, positioned ahead of our Sun in the directionof Galactic rotation at a distance of about 795 ±
25 pc(Section 5). In projection on the sky, NAP is locatedadjacent to the massive star-forming regions of CygnusX, but it is generally thought that NAP is nearer thanCygnus X. A study of molecular clouds in the Solarneighborhood by Zucker et al. (2020) and Alves et al.(2020) has confirmed that this complex is part of a stringof star-forming regions stretching from the Orion Molec-ular Clouds to Cygnus X.The NAP region provides a useful laboratory forstudying both star formation in individual clouds aswell as the evolution of the whole complex. Section 2presents the YSO catalogs,
Gaia data, radio data, andspectroscopy used for this study. Section 3 describes themethods used to identify stellar members and reject con-taminants. Section 4 describes the properties of multiplestellar groups. Section 5 computes distance estimates forthe components of the association. Section 6 describesthe structure of the molecular clouds and their relationto the stars. Section 7 provides an analysis of the stellarkinematics. Section 8 identifies new candidate members The boundaries of the named subregions are given by Rebull et al.(2011) and Zhang et al. (2014). As a point of comparison, the nearest O4 star, ζ Pup, is less thanhalf as far as the Bajamar Star (Howarth & van Leeuwen 2019). identified from the
Gaia catalogs. Section 9 discussesthe formation and dynamical evolution of the system.And, Section 10 summarizes the main conclusions. DATA2.1.
Initial Stellar Catalog
We searched for lists of candidate YSOs from pub-lished studies of the NAP region using the “YSO Cor-ral” (Hillenbrand & Baliber 2015), a curated databasefor YSOs.The previous studies have used a variety of techniquesto identify candidate members, as listed in Table 1.These included selection based on H α emission, infraredexcess, X-ray emission, placement on color-magnitudediagrams, and spatial clustering. H α emission and in-frared 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 ofthese methods yields different types of contaminants andrates of contamination, so we aim to re-assess mem-berships of their candidate members using astrometryfrom Gaia , as described below. Contaminants to YSOsearches in the Galactic plane include both field starsand extragalactic sources, in particular dusty (post-)asymptotic giant branch (AGB) stars and star-forminggalaxies may be major sources of contaminants for se-lection based on infrared excess (Robitaille et al. 2008),while active galactic nuclei (AGN) and foreground ac-tive M dwarfs may dominate X-ray contaminants (Brooset al. 2013). 2.2.
Astrometric Data
We obtained stellar astrometry from
Gaia ’s seconddata release (DR2; Gaia Collaboration et al. 2018),which provides 5-parameter astrometric solutions for1.3 billion objects to as faint as G ≈
21 mag (Lindegrenet al. 2018). These solutions provide source positionsin right ascension ( α ) and declination ( δ ), parallax ( (cid:36) ),and proper motion ( µ α (cid:63) , µ δ ). For DR2, most stars inthe direction of NAP were observed during ∼
15 visibilityperiods; typical formal uncertainties for a G = 15 magstar are ∼ ∼ − inproper motion. In addition to the formal uncertainty, ∼ ∼ − correlated systematicuncertainties affect these measurments (Lindegren et al.2018). There have been a variety of efforts to estimatethe systematic zero-point offsets (e.g., Leung & Bovy2019, and references therein); in most cases we try towork with the observed quantities (i.e. (cid:36) ), but whendistances are required we apply a 0.0523 mas parallaxcorrection estimated in the aforementioned paper. We follow the convention of Perryman et al. (1997), and define µ α (cid:63) ≡ µ α cos δ . This quantity is called PMRA in the
Gaia
DR2tables. tar and Gas Kinematics in NGC 7000/IC 5070
Gaia catalog both for cross matching to thepreviously identified YSO candidates as well as to searchfor new member candidates. For cross matching, we usea match radius of 1.2 (cid:48)(cid:48) and select the nearest source as amatch, yielding 1939 matches out of ∼ Gaia detection limits. To estimate the rate of false matches,we artificially shifted the DR2 coordinates 2 (cid:48) north andran a cross match using these coordinates. We found52 matches to the shifted coordinates, suggesting an in-correct match rate of up to 3%. This rate should beregarded as an upper limit on match errors, because atrue counterpart, if it exists, would likely take prece-dence over a field star if both lie within the match ra-dius.We apply cuts to the
Gaia data to ensure that themeasurements are sufficiently precise to be useful fordistinguishing between members and non-members. Wefollow the recommendations of the Gaia
Data Pro-cessing and Analysis Consortium and use only sourceswith renormalized unit weight error (RUWE) less than1 .
4. The
Gaia catalog also tabulates excess noise inthe astrometric fit. Higher values of excess noise (e.g., > ≤ astrometric sigma5d max ≤ . (cid:36) > Millimeter Observations
We obtained CO J = 1–0 (110.201354 GHz) obser-vations of the molecular clouds associated with the NAPregion using the Five College Radio Astronomy Obser-vatory (FCRAO) 14-m telescope in 1998 June. The COmap covers most of the region of interest in this study.The observations were obtained with the SecondQuabbin Optical Imaging Array (SEQUOIA; Ericksonet al. 1999). The array contains 16 pixel elements ar-ranged in a 4 × (cid:48)(cid:48) on thesky. At the time of the observations, 12 of the pixels were functional. The spectrometer contained 512 chan-nels with a bandwidth of 40 MHz, which provided achannel spacing of 0.21 km s − . The FCRAO telescopehas a full width at half maximum (FWHM) angular res-olution of 47 (cid:48)(cid:48) at the CO J = 1–0 frequency. Theforward scattering and spillover efficiency of the tele-scope at the observed frequency was η FSS = 0 .
7, whilethe main-beam efficiency was η mb = 0 .
45 (Heyer & Tere-bey 1998). The data were obtained in position-switchingmode with spectra that were sampled every 44 (cid:48)(cid:48) on thesky, or approximately the FWHM resolution. The typi-cal RMS noise in the spectra is ∆ T A ∗ = 0 .
22 K per chan-nel. The maps used for our analysis are in units of T R ∗ ,calculated from the atmosphere-corrected antenna tem-perature 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 thanan arcminute in the NAP region.Our CO data are similar in sensitivity and spatialresolution to a CO map published by Zhang et al.(2014) using the Purple Mountain Observatory Delingha13.7 m telescope. They also provide CO and C Omaps.2.4.
Spectroscopic Observation of the Ionizing Source
Finally, high dispersion optical spectra of the BajamarStar were taken on 1 December 2019 and 3 January 2020(UT) for the purposes of confirming the existing spectraltyping that was based on low resolution spectra, andassessing the source radial velocity. These data weretaken with the Keck I telescope and HIRES (Vogt et al.1994) and cover ∼ − R ≈ , VALIDATION OF MEMBERSHIP WITHASTROMETRYWe can evaluate membership of candidate YSOs byobserving whether they are at the same distance andmoving with the same kinematics as the rest of the mem-bers in the NAP region. In Figures 1–3, we use the fol-lowing color scheme: gray – stars at the wrong distance;green – stars at a similar distance but with the wrongproper motion; magenta – probable astrometric mem-bers; goldenrod – stars with uncertain classification.The distribution of parallaxes of YSO candidates (Fig-ure 1, left panel) has at least two modes as well as a tailof objects with high parallax. The expected distance toNAP based on previous measurements is 500–1000 pc(e.g., Reipurth & Schneider 2008), which coincides withthe second peak on the histogram (magenta tic mark).The left-most peak can be attributed to contaminantswith small or zero parallax, including giant stars andextragalactic sources. The tail of objects to the rightincludes many bright objects with high-precision paral-lax measurements, so these are clearly foreground stars.To quantitatively divide stars based on parallax, wemodel the parallax distribution using a Gaussian mix-ture model, for which we use the Bayesian Information
Kuhn et al. −1 0 1 2 3 4 parallax [mas] nu m be r o f ob j e c t s a) l llll l llll ll l lll ll ll ll l ll ll lll lll l lll ll ll ll lll ll lll ll ll ll l ll ll ll ll ll lll ll ll l ll ll lll ll l ll ll ll l lllll llll l lll l l ll llll ll l ll l lll ll l llll lll l ll ll l lll l ll ll llll lll ll l lll lll ll ll l llll l lll ll ll lll l lll ll ll ll l
10 5 0 −5 −10 −15 − − PML [ mas yr - ] P M B [ m a s y r - ] b) Figure 1.
The panels illustrate the two step process for using
Gaia astrometry to validate YSO candidate membership. a) Thedistribution of parallaxes (histogram) is fit with a Gaussian mixture model (black line). The centers of the three Gaussians areindicated by the vertical marks above the curve. Objects with a probability >
50% of being a member of the middle component(associated with NAP) are selected as indicated by the striped green/magenta region. b) Proper motions µ (cid:96) (cid:63) (PML) and µ b (PMB) for the objects selected in step (a). The distribution of these sources is fit with a two-dimensional Gaussian mixturemodel, and the objects with a probability >
50% of being a member of the middle density peak are colored magenta. Thesecomprise our sample of 395
Gaia -validated NAP members.
Criterion (BIC; Schwarz et al. 1978) to determine thenumber of components. For Gaussian mixture mod-els the BIC can be viewed as an approximation for theBayes factor which is useful for the statistical problemof model selection (Everitt et al. 2011). We find thatthe distribution can be well modeled with 3 Gaussians,as marked in Figure 1, with ∆
BIC = 10 compared tothe next best model, implying strong preference for thissolution. From low parallax to high parallax, the firstcomponent models the peak that we associated withbackground contaminants, the second component cor-responds to the group of stars associated with NAP andhas a mean of (cid:36) = 1 .
25 mas, and the third compo-nent has a slightly higher mean parallax but a muchbroader distribution that appears to approximate theshape of the high-parallax tail. Objects that lie be-tween 0 . < (cid:36) < .
61 mas (green and magenta stripedregion) have a ≥
50% probability of belonging to the sec-ond component; all other sources are classified as non-members.In the parallax selected sample, cross-contaminationfrom field stars that happen to have similar parallaxesto the NAP region is inevitable. To reduce this con-tamination, we perform a second classification step, thistime using proper motion. Figure 1 (right panel) revealsa clump of sources with similar motions, and a halo ofobjects with significantly different motions. We showthese motions in Galactic coordinates ( µ (cid:96) (cid:63) , µ b ) to em-phasize that the proper-motion dispersion in the halo islargely parallel to the Galactic plane, and thus repre-sents orbital motions of field stars in the Galaxy whichare dynamically hotter than newly formed stars and thus have larger dispersions in dynamical phase space.We subdivide the sources in proper-motion space usinganother Gaussian mixture model, finding two compo-nents (∆BIC=7), corresponding to the clump (smallerdispersion) and the halo (larger dispersion). We clas-sify objects with ≥
50% probability of belonging to thefirst component as members (magenta points), while theothers are non members (green points). The membershave mean proper motions of µ (cid:96) (cid:63) = − .
35 mas yr − and µ b = − .
15 mas yr − ; other parameters of the classifierare given in Appendix A. The mixture model suggests a3% residual contamination rate among the stars classi-fied as members. We note that selection based on propermotion comes at the cost of omitting stars with highproper motions that could have been ejected from thestar-forming region.Figure 2 (left) shows the classified sources on a dia-gram of Gaia absolute magnitude versus color. Whencompared to theoretical PARSEC isochrones (Bressanet al. 2012), the members nearly all lie above the 3 Myrisochrone, and many of them lie above the 1 Myrisochrone. In contrast, the objects identified as non-members are scattered throughout the diagram, withmany lying below the 10 Myr isochrone, while othersare located in the region of post–main-sequence giantstars, and others lie in similar locations to the pre–main-sequence stars. This diagram confirms that nearlyall of the astrometrically validated members are pre–main-sequence stars with ages < Gaia ’s G and G RP bands, so we show three ap- tar and Gas Kinematics in NGC 7000/IC 5070 ll l l lll llll ll l llll ll lll l l ll l lll ll l ll l ll lll ll l ll lll l llll l l lll llllll ll l llll l ll ll lll lll ll lll l llll ll lll lll ll lll l ll l l lll l lll l ll l lll llll llll lll ll lll l lll llll l lllll ll l ll l ll lll llll ll l lllll llllll lll lll ll lll l l ll ll l ll l lllll l ll lll ll l lll ll lll l lllllll l ll ll l llll llll ll lll l lll l l ll ll lll lll l l ll lll ll ll l ll llll lll ll lll l lll lll lll ll l lll l lll lll lll llll lll ll lll ll ll ll lllll ll ll lll l ll l lll l ll llll l ll ll lll ll lll ll llll l l ll l llll l ll l l ll ll l ll lll ll lll l lll ll ll ll lll lllll lll l ll ll lll lll ll lll l lll ll lll lllll l l l ll lll lll l ll ll lll ll llll ll lll l lll ll lll l ll l llll l lll lll l ll lll ll l lll ll ll l ll ll lll lll ll ll l l lll ll lllll l ll lll ll lll l ll ll lll ll lll ll l l ll lll ll l l lll ll lll l ll l lll ll l lll lll l ll lll l ll ll ll l l ll ll l lll llll ll llll llll lll ll l ll l lll l ll ll l llll ll lllll l lll lll ll l lll lll ll ll lll ll l lll l ll l ll l l lll lll lll ll l ll l lll ll lll llll ll l llll lll ll l ll l lll l l ll l ll ll llll llll l l l l ll lll ll lll lll ll ll ll l ll ll ll ll l lll lll ll lll ll lll llll lll llll l ll l lll l lll l llll l ll l lll llllll ll l ll ll ll ll ll l ll ll lllll ll ll ll l lll ll ll lll l ll llllll l llll ll ll ll l ll l ll ll lllll lllll lll ll l l ll llll ll ll lll l l llll lll ll lll ll l l lll lll llll ll lll l ll ll l lll ll ll ll lll l lll lll llll lll llll l ll ll l l ll ll lll ll G − RP [mag] G a i a G + l og ( p l x [ m a s ] ) − l l ll llll ll l llll ll ll l l ll ll ll l ll ll l ll lllll lllll l l llll l ll ll lll llll lll llll llllll lll ll l lll l lll l ll l lll lll lll ll ll ll llll llll llll ll l ll ll lll llll ll l llll llllll lll lll ll ll l l ll ll ll l llll lll ll l lll lll l lllllll l lll lll llll ll lll l ll l l lll lll ll l ll lll lll l ll llll ll ll lll l ll ll llll l lll l lll lll lll ll lll l lll ll lll l lll ll ll l ll l lll ll l l ll l ll ll lll ll ll lll lllll lll l ll lll lll ll llll ll lll ll ll l ll ll lll l ll lll ll l ll ll lll ll ll lll ll ll l ll l l ll l l ll l llll ll l lll ll l l llll l llll ll ll l ll ll l l lll lllll l l ll ll l lll l llll l lll ll l ll ll lll l l ll ll ll ll ll lll ll l l l lll lll llll l ll l ll ll l lll lll ll llll ll lll lllll ll ll l lllll ll l ll ll ll ll ll llll ll lll lll lll ll ll ll ll ll l lll lll l lll l lll l l ll l lllll ll l lllll ll l lll ll l ll l l ll lllll lllll l llll l ll ll lll llll lll llll llllll lll ll l lll l lll l ll l lll lll lll ll ll ll llll llll llll ll l ll ll lll llll l lll lllll lll lll l ll l l ll ll ll l lll lll ll ll lll l lllllll l lll lll llll ll ll ll l l ll lll ll l ll lll lll l ll llll ll ll lll ll ll lll l lll l lll lll ll l ll lll ll ll lll ll ll l ll l lll ll l l ll l ll llll ll l llll l l lll lll llll ll lll lll ll llll ll lll ll ll ll ll ll ll ll ll l ll l l ll l ll l llll l ll l lll ll ll llll ll ll ll l l l ll lll l ll lll ll ll l l l A V = mag1 Myr3 Myr10 Myr l llll ll l l lll ll l ll lllll ll lll lll ll lll ll ll l ll lll ll ll l ll ll ll ll ll ll l ll ll l lllll ll ll ll lll lll l l l ll ll lll l llllll ll l l lll llll l l l ll ll l lll lllll l ll llll l ll ll ll l llll ll ll lll ll ll ll ll lllll l lll ll lll l llll lll ll ll lll lllll ll lll lll lll l l ll llll l lllll lll l lllll l ll ll l lll ll l ll lll lll lll llll l ll l ll l l llll lll ll ll ll ll l l l ll l ll l llll ll ll lll lll ll lll ll ll llll l ll l ll lll l l ll ll llll ll ll ll ll l lll lll l lll l lllll ll llll ll ll ll ll l lll ll ll l llllll l ll ll ll ll lll l lll ll l l ll lllll l lll l ll lll lllll ll lll ll l l ll ll ll ll l lll ll lll l ll ll lll l l lll ll ll l ll lll ll lll ll l ll lll ll l l l l lll llll ll ll ll l ll lll lll ll ll ll ll ll ll lll lll lll lll l ll l l lll l lll lllll ll lll l l ll l ll lll lll ll ll lll l ll lll ll ll ll llll lll ll ll l ll ll l lll lll ll l ll l l l ll ll l llll ll l ll lll lll ll ll ll ll llll l ll ll llll l ll lll l l ll l lll l lll lll l ll l lll lll lll lll l ll lll lll llll l l ll ll l lllllll ll ll ll ll lll lll ll l lll l l l lll ll ll l ll lll ll ll lll l ll l lll lll lll l ll llll l ll ll ll lll ll l lll l l lll l lll lll lll lll lll ll l lllll llll lll l l l ll lll ll ll ll l l l lllll ll l l llll ll lll llll ll l llllll ll lll lll ll ll l l llll l lll lll lll lll ll lll ll lllll lll ll ll l llll ll ll ll l l ll l lll lll l l llll lll lll l l llll l llll l l ll lll lll l lll l ll ll ll lll llll ll l l l lll ll lll l lll llll ll l lll lll lll ll ll lll ll ll llll lll ll ll ll l llll ll lllll l l llll lll ll lll lll llll lll lll ll l lll l ll lll lll ll ll ll l lll lll lll llll lll lllll ll lll l ll ll lllll . . . . . . H - K s [mag] J − H [ m ag ] l ll ll ll l ll lll ll ll ll ll lll lll lll lll ll l ll ll ll ll ll ll l ll l lll l lllllll l lll llll l ll l ll llll l lllll ll ll ll l llll l ll ll ll ll ll llll llll l ll llll ll l ll lll lllll ll ll lll ll l ll llll lllll ll lllll ll ll l lll lll lll lll lll ll l ll l l lll lll llll ll l l l ll l l l lll ll lll llll llll l llll l ll l lll l l ll ll ll lll l ll l l lll ll ll llll l ll lll ll l l ll lll lll l ll l l ll ll l ll ll ll ll ll ll l lll ll ll ll lll l ll lll lll l ll l l ll ll ll l ll ll ll l llll l ll ll ll llll lllll lll l ll llll ll ll l lll l ll l lll l l lll l ll lll ll l l lllll lll l ll ll l ll ll l l ll l lllll ll llllll llll ll ll l lll lll ll llll ll l ll ll ll ll lll ll ll l l l lll lll llll ll l lllll l lll lll ll l ll l ll l lll llll ll ll l ll lllll l lll l ll ll llll ll lll lll llll l ll l lll lll lll lll l llll lll l l ll l l l ll l ll lll ll l lll ll lll lll ll ll l ll ll ll ll ll ll l ll l lll l lllllll l lll llll l ll l ll llll l lllll ll ll ll l llll l ll ll ll ll ll llll llll l ll llll ll l ll lll llll l ll ll ll ll llll llll ll lllll ll ll l ll lll lll llll ll l ll l l lll lll llll ll l l l ll l llll lll llll llll l llll l ll l lll l l l ll l lll l ll l lll ll ll llll l lll l l ll ll lll l ll l l ll ll l ll lll ll ll ll lll ll ll lll llll l ll ll ll ll l ll l ll lll l lll ll lll llll l l ll lll ll ll l lll l ll lll l l lll l lll l lll lll llll l lll ll ll l ll lll l l ll lll ll ll ll l A V = mag Figure 2.
Left:
Gaia color-magnitude diagram for YSO candidates classified in this work. Members are magenta points,sources excluded based on parallax are gray points, and sources excluded based on proper motion are green points. Only sourceswith (cid:36)/σ (cid:36) > . A V = 2 mag for pre–main-sequence stars with unreddened colors of G − RP = 0, 0.6, and 1.3 mag. The two O stars stars are marked with the cyan asterisks; the Bajamar Star is the one with thevery red color. Right: 2MASS color-color diagram using the same symbols as the left panel. The black curve indicates the locusof unreddened stellar colors from the PARSEC models, whose shape is similar for young pre–main-sequence stars but extends ∼ J − H around and beyond the peak (Herczeg & Hillenbrand 2015).
316 315 314 313 312 . . . . . . . D e c [ deg ] R.A. [deg] llllllllllllllllllllllllllllllllllll lllllll ll lll lllll llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll l l lllllllllllllllllllllllllllllllllllllllllllllll lllllllllllllll l
Rejects
316 315 314 313 312 llllllllllllllllllllllllllllllllllllllllllllllllllll lllllllllllllll lllllllllllllllllllllllllllllllllllllllllllll l l lll l lll ll l llllll llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllll lll ll ll ll ll l l lll lll lllll l l ll ll ll llllll llllll lllll llll llllll l lll llllll llll l llllll llllllll llllll l llllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllllllllllllllllllllll l l lll l llllll llll l lll llll ll ll ll lllllll llll ll llll ll lllllll l lllllll ll lllll lll lll llllll lllllllllllll ll l llll l l l llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll l ll llll lll lllll llllll llll ll llll l lll lll llllll l lll lll llllllllllllll lllll l lllllllllllllllllll lllll llll llllllll llllll lll llllll l llllll lll llllll l llllllllllll lllll llllllllll lll lllll ll lllll lll l l lll lll ll lllllllllllllllll lll lllllll lllll lll l l lllll lll ll l lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll
Uncertain
316 315 314 313 312 llllllll llllll llllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll lllllllll lllllllll ll ll l llllllll lllllllllllll lll lllllllllllllllllllllllllllllllllllllllllllll lllllll l lllllllll l
Members
Figure 3.
Spatial distributions of YSO candidates stratified by their classification as contaminants or members. The left panelshows the distributions of different classes of contaminant: objects rejected by parallax in gray (background to the NAP as lightcrosses and foreground as dark diamonds) and objects rejected by proper motion in green. The middle panel shows objectswhere
Gaia -based classification is uncertain because the source does not exist in the
Gaia catalog or it does not meet our qualitycriteria. The right panel shows objects classified as bona fide members. In black we show contours of the optical nebulae fromthe DSS image. Much of the spatial clustering in the middle panel is reminiscent of structure in the right panel, indicatingthat our techniques are selecting only a portion of the true young star population, but there is also a broad spatial distributionreminiscent of the left panel, indicating a mix of contaminants in the uncertain population too.
Kuhn et al. proximate vectors for different G − RP colors assuminga typical pre–main-sequence stars spectrum. From the J − H vs. H − K s color-color diagram (Figure 2, right)using photometry from the Two Micron All Sky Survey(2MASS; Skrutskie et al. 2006), we estimate that typi-cal reddening is E ( J − H ) ∼ . A V ∼ Members and Contaminants
Out of the sample of YSO candidates, 395 objects areconfirmed as members, while 814 objects are reclassifiedas non-members, and 2264 sources either do not have
Gaia counterpart or do not meet our quality criteria andtherefore have uncertain classifications. The numbers ofcandidates in each category are tabulated in Table 1 foreach of the published lists of YSO candidates. Table 2provides a list of the astrometric members.It is not surprising that nearly all historical studieshave some contamination. For instance, out of the 68objects from the pioneering study of the NAP region byHerbig (1958),
Gaia provides sufficiently high qualitymeasurements for 27 of them. Of these stars, our analy-sis has validated the membership of 21 objects, while re-jecting membership for 6 others. The rejected membersinclude LkH α Gaia ’s re-quirement that a source must be sufficiently bright inthe optical may induce a selection bias toward highercontamination rates. For example, out of the 1,286
Spitzer /MIPS sources classified as YSOs by Rebull et al.(2011) on the basis of infrared excess, only 30% metour
Gaia criteria for inclusion in the analysis. Of those30%, 131 were rejected and 251 were validated. How-ever, this ratio is not representative of the whole sampleof 1,286 MIPS candidates because the YSOs with themost prominent infrared excesses tend to be very youngand deeply embedded, hence not optically visible.The spatial distribution of candidates classified mem-bers, non-members, or uncertain membership is shownin Figure 3. These points are over-plotted on an outlineof the optical nebulosity making up the North Americaand Pelican components of the nebula. Non-membersare widely dispersed throughout the whole region, whileuncertain members are more tightly clustered, and val-idated members are the most tightly clustered. Thissuggests that the spatial distribution of stars is domi-nated by clustered groups, rather than distributed stars.Although a few members are located a degree or moreaway from the main groups, there is no significant non-clustered population scattered throughout the entire ∼ ◦ diameter region. Our sample of 395 objects is clearly incomplete, miss-ing both objects that were not selected in previous stud-ies (e.g., pre–main-sequence stars without disks) as wellas objects that could not be evaluated by Gaia . There isalso evidence that our classifier has rejected a few legiti-mate YSOs. From a sample of 41 high-confidence YSOsstudied by Findeisen et al. (2013), 4 were rejected bythe classifier, likely spuriously. Altogether, it is likelythat the NAP region contains at least several thousandstellar members. This assertion is based on
Spitzer ob-servations that have collectively identified ∼ Gaia . Furthermore, in Section 8 we iden-tify > Gaia astro-metric data. Strictly speaking, our membership classifi-cation is for the
Gaia source, so the chance alignment ofa mid-infrared YSO with a
Gaia field star would likelyresult in rejection of the source.3.2.
Properties of the Ionizing Sources
Both the Bajamar Star (the primary ionizing source)and HD 199579 pass our
Gaia membership criteria inSection 3. However, the parallaxes and proper mo-tions of both of these O stars place them near theedges of parameter space for objects identified as mem-bers. For HD 199579 ( (cid:36) = 1 . ± .
06 mas, µ (cid:96) (cid:63) = − . ± . − and µ b = − . ± . − ), itsparallax places it behind most other cluster members. Itlies on the northeast extreme of the NAP region and ismoving in this direction relative to the rest of the NAPmembers.The Bajamar Star ( (cid:36) = 1 . ± .
08 mas, µ (cid:96) (cid:63) = − . ± . − and µ b = − . ± . − ) hasa parallax measurement that is 2 σ greater than the me-dian parallax of the system, and the µ (cid:96) (cid:63) proper motionis more negative than that of most other members. The Gaia data do not indicate any obvious problems withthe astrometry; the RUWE is 1.04 which is within therecommended < II W λ = 1 . I W λ < .
18 ˚A), He I W λ < .
13 ˚A), and He I W λ < .
08 ˚A).We can also confirm the statement in Ma´ız Apell´anizet al. (2016) that the different absorption line specieshave different radial velocities. However, it is not clearbased on our spectra that the source is a spectroscopicbinary. While there is evidence for asymmetries in theline profiles, especially H I and He I , a high velocity tar and Gas Kinematics in NGC 7000/IC 5070 IC 5070NGC 7000 LDN 933LDN 935 1 degHD 199579Group F Group E Group D Group CGroup BGroup A NEV1057 CygHBC 722 Bajamar Star
Figure 4.
Members, as indicated in the right panel of Figure 3, are over-plotted on a DSS red-band image of the region andcolor coded based on the kinematic group they belong to. Note that several outlying members are cropped out of the image(see Figure 3) to show enough detail in the central region. Several interesting stars are indicated, including the two O stars andtwo well-studied outbursting stars. The nebulous regions are those that were outlined in Figure 3 – on the left the nebulosityresembles the shape of North America and on the right it resembles a pelican with its beak facing the Atlantic Coast of NorthAmerica. The naming of subregions (e.g., Gulf of Mexico, Atlantic, Pelican’s Neck) has historically followed these analogies. wind, as is expected from such a massive star, wouldproduce similar features. Furthermore, a cross corre-lation analysis of our HIRES spectrum does not revealindications of two peaks, through the lines are extremelybroad due to the rapid rotation and possibly other ef-fects, and it certainly would be possible to mask two setsof broad lines within the profiles.Using τ Sco as a radial velocity standard, we derive v helio = 11 . ± . − for our first observation and v helio = − . ± . − , both based on the singleHe II −
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 veloc-ity shift. Measured directly, rather than going througha narrow-line star, the shift in the strong He II ± − . This certainly suggestsbinarity. However, an order containing He I IV −
40 km s − while He II − difference and while both He I II − . In H α and H β the shift between the spectra is about 15 ± − . Itremains unclear whether we should interpret the mea-surements as a single star having the same velocity as themolecular gas and strong wind that manifests in confus-ing absorption velocities, or as a binary with the moremassive star producing much of the He II absorptionand both stars contributing to the relatively blue-shiftedmetal, He I and H I profiles. KINEMATIC GROUPSThe distribution of
Gaia -validated NAP members isneither smooth nor centrally concentrated, but has amore complicated, clumpy structure that is apparent
Kuhn et al. both in the stars’ spatial positions (Figure 4) and intheir distributions in proper-motion space (Figures 5and 6). Spatial clustering and sub-clustering of youngstars has been appreciated for decades (e.g., Carpenter2000; Allen et al. 2007), but only recently have the kine-matic data achieved comparable fidelity (e.g., Gonz´alez& Alfaro 2017; Da Rio et al. 2017). This analysis is nowpossible in a larger number of star-forming regions with
Gaia
DR2+ astrometry.4.1.
Cluster Analysis Algorithm
When analyzing the stellar population of the NAP re-gion, it is convenient to divide the stars into subclusters.As for all cluster analysis problems, multiple possiblestrategies can achieve this, with no one approach being“best” (Everitt et al. 2011). We find that a Gaussianmixture model with hierarchical combination of clus-ters, outlined below, gives reasonable results that areastronomically useful.We use a Gaussian mixture model analysis imple-mented by mclust (Scrucca et al. 2016) in R to identifygroups of stars in 4-dimensional position–proper motionspace. The analysis is performed on the variables (cid:96) , b , µ (cid:96) (cid:63) , and µ b , each of which is normalized by subtractingthe mean value and dividing by the standard deviation.The program then selects the number of clusters, thefree parameters of the Gaussian components, and thevalues of the parameters via minimization of the BIC.For modeling the clustering of stars, the mixture-modelapproach has several benefits, which include the abilityto identify clusters that differ in both compactness andnumber of members, the ability to assign stars from thedense cluster core as well as the less-dense cluster wingsto the same cluster, and the ability to deal with over-lapping clusters. On the other hand, stars clusters aretypically not well fit with Gaussian distributions, whichcan lead to multiple Gaussian components used to fit asingle cluster (Kuhn & Feigelson 2019). For the NAPstars, mclust selects a model with eight ellipsoidal com-ponents with equal shape and orientation, where severalof the Gaussian components are significantly overlap-ping or nearly concentric.The mclust package includes the function clustCombi to deal with such cases by hierarchically merging mul-tiple Gaussian components into single clusters using anentropy criterion (Baudry et al. 2010). For our case, wefind a large change in normalized entropy going from 6to 7 clusters, so we use this threshold to cut the hier-archical tree at 6 groups. The stars assigned to each ofthese six groups, labeled A through F, are over-plottedon the optical DSS image of the NAP region in Figure 4.4.2. Groups in Position–Proper-Motion Space
The results of the cluster analysis algorithm can bequalitatively checked for reasonableness by examinationof the Group assignments in Figures 4–6. In any oneof the projections, the edges of some groups overlap. However, the stars that overlap in one projection tendto be different from the stars that overlap in a differentprojection, and the groups are mostly well separated.In projection on the sky (Figure 4), Group A is fairlyisolated in the southwest corner of the region. GroupsB, C, and D form a conglomeration in the northwest-ern portion of the nebula, more-or-less corresponding tothe regions known as the “Pelican’s Hat,” the “Pelican’sNeck,” and the “Atlantic,” respectively. Groups E andF lie in the southeastern/eastern region – Group E iscomposed of several clumps of stars on the dark “Gulfof Mexico” cloud, while Group F is a single clump tothe northeast of the “Gulf of Mexico” superimposed onan area with bright nebulosity.The plot of α vs. µ α (cid:63) (Figure 5, upper left) shows theclearest separation of the stars into distinct groups. Inthe 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 muchmore tightly clumped in α and µ α (cid:63) . On the right side ofthe diagram, Group B appears partially entangled withGroup C, and Group A appears partially entangled withGroup B. However, both Groups A and B are looser thanGroup C. The spatial separation of A from the rest ofthe groups make it clear that A is distinct. However,whether B could be considered part of the same groupas C is less clear. Group B contains only a few starsfrom our sample; however, these stars coincide with thedeeply embedded “Pelican’s Hat” cluster of stars discov-ered by Spitzer (Guieu et al. 2009; Rebull et al. 2011),most of which are unseen by
Gaia .In the plot of δ vs. µ δ (Figure 5, upper right), thediagonal swath formed by Group D can still be seen, butthe other groups appear distributed in a more verticaldirection rather than diagonal. This plot enhances theseparation between Groups A, B, and C. Groups C andD are less cleanly separated, but the stars in Group C donot continue the same diagonal tend as Group D. Thelower panels in Figure 5 show the other combinations ofposition and proper motion. Finally, on the plot of µ α (cid:63) vs. µ δ (Figure 6), Groups D and E are in the center,Group A is at the top extreme, Groups B and C are atthe bottom extreme, and Group E is at the right-mostextreme.Intriguingly, the Bajamar Star is assigned to Group Deven though it is spatially closer to stars from Group Ein the “Gulf of Mexico.” This is a result of the propermotion of the Bajamar Star, which is rapidly movingtoward Group E away from the center of Group D. Thisproper motion does not match the proper motions ofnearby stars in Group E; however, it does fit with theproper-motion gradient seen for stars in Group D (Fig-ure 5).Our division of stars into kinematic groups is meantto be useful for understanding the structure of the NAP tar and Gas Kinematics in NGC 7000/IC 5070 l l ll l l l lllllllllllllllllllllllll llll l ll l llll l ll ll l l ll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll l ll llll lll llllllll ll lll ll l l lll l l l l l l llll ll lll l lllllllllllllllll lllllllllllll l l l ll ll llll ll llllllllllllllllllllllllllllllllllll lll l ll lllllll lllllll ll ll ll llllllllllll ll
315 314 313 312 − − − R.A. [deg] P M R A [ m a s y r - ] l l lll ll llll ll ll lll ll ll lll lll ll l ll ll l ll l l ll lll ll l llll l l ll ll l lll ll l lll lll l lll llllll lllll lll lll lll ll lll l ll lll l lllll l ll ll llll ll lll l lll l ll lllll llll l llll ll lll lll lll ll ll llll lll ll l ll lllll ll lllll ll lll l ll lll lll lll lllll lll l ll llll ll llllll l llll llll llll l ll llll lll l lllllll ll ll ll ll l lll ll ll lll ll l lllllllll lll ll llllll lll ll l ll ll l lll ll ll lll lllllll llllll lll l llll lll lll ll ll ll lll l l lll ll lll ll llll l l lll lllllll ll lll ll lllll − − − − − Dec. [deg] P M D E C [ m a s y r - ] l l l l ll l l l lllllllllllllllllllllllll llll l ll l llll l ll ll l l ll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll l ll llll lll llllllll ll lll ll l l lll l l l l l l llll ll lll l lllllllllllllllll lllllllllllll l l l ll ll llll ll llllllllllllllllllllllllllllllllllll lll l ll lllllll lllllll ll ll ll llllllllllll ll
315 314 313 312 − − − − − R.A. [deg] P M D E C [ m a s y r - ] lll ll llll ll ll lll ll ll lll lll ll l ll ll l ll l l ll lll ll l llll l l ll ll l lll ll l lll lll l lll llllll lllll lll lll lll ll lll l ll lll l lllll l ll ll llll ll lll l lll l ll lllll llll l llll ll lll lll lll ll ll llll lll ll l ll lllll ll lllll ll lll l ll lll lll lll lllll lll l ll llll ll llllll l llll llll llll l ll llll lll l lllllll ll ll ll ll l lll ll ll lll ll l lllllllll lll ll llllll lll ll l ll ll l lll ll ll lll lllllll llllll lll l llll lll lll ll ll ll lll l l lll ll lll ll llll l l lll lllllll ll lll ll lllll − − − Dec. [deg] P M R A [ m a s y r - ] Figure 5.
Plots of position coordinates α and δ vs. proper motion coordinates µ α (cid:63) (PMRA) and µ δ (PMDEC). NAP membersare color coded by group using the same hues as Figure 4. The black points in the top panels with error bars indicate themedian formal uncertainties in µ α (cid:63) and µ δ . region, but it is not the full description of the clusterstructure. For example, nearby groups may have formedfrom related star-formation events at different locationson the same cloud, and thus, from an astrophysical per-spective, it is ambiguous whether these should be consid-ered to be the same group. Furthermore, if the regionwere nearer (e.g., at the distance of Taurus; Luhman2018; Fleming et al. 2019) we might be able to detectfiner velocity differences with Gaia that could be usedto subdivide the groups further, but if the region weremore distant (e.g., at the distance of the Carina Nebula;Kuhn et al. 2019) the measurement uncertainties in theproper motion would be poorer and, thus, it might nothave been possible to distinguish the groups. THREE DIMENSIONAL STRUCTURE OF THESTELLAR POPULATIONThe region is sufficiently close that differences in dis-tance to the groups begin to become detectible by
Gaia .However, these differences are subtle enough that caremust be taken in estimating mean parallaxes of eachgroup. We calculate the mean parallaxes using a max-imum likelihood method that takes into account boththe heteroscedastic measurement errors as well as thetruncation of the sample at (cid:36) min = 0 .
86 mas and (cid:36) max = 1 .
61 mas (Section 3). The log-likelihood as0
Kuhn et al. l l llll ll lll lll ll l lll l llll l l lll ll ll lll ll ll l ll l l llll ll l llll lll ll ll l llll lllllll l ll lll l llll lll lll lllll l lll ll llll lll ll llll ll lllll ll lll l ll ll lllll ll lll ll l lll lll l l llllll lllllll ll ll l lll ll ll ll l ll l lll l ll lll lll ll llll lll ll l lll ll llll ll l l l ll lll ll ll llllllll ll ll ll l ll l l l ll ll lll l l ll ll lll l ll lll lll ll lll ll ll l ll ll lll l lllll l llll l ll ll l ll lll llll ll l lllll ll ll lll lllll lll lll lllll l lllllllll ll l l l ll ll l l lll l lll llll ll ll l l ll ll l − − − − − PMRA [ mas yr - ] P M D E C [ m a s y r - ] Figure 6.
Scatter plot showing µ α (cid:63) vs. µ δ . NAP membersare color coded by group as in Figures 4 and 5. a function of mean parallax, ¯ (cid:36) , is given by L ( ¯ (cid:36), σ ) = (cid:88) i ln φ ( (cid:36) i ; ¯ (cid:36), σ i + σ ) (cid:82) (cid:36) max (cid:36) min φ ( (cid:36) i ; ¯ (cid:36), σ i + σ )d (cid:36) (1)where φ denotes a Gaussian distribution and (cid:36) i are themeasured parallaxes with uncertainties σ i . This modelallows stellar groups to have an intrinsic depth, char-acterized by the parameter σ . We use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm to find thevalue of ¯ (cid:36) that maximizes this likelihood, and use theHessian matrix calculated by the R function optim toestimate confidence intervals.The mean parallaxes and their confidence intervals arereported in Table 3 and plotted in Figure 7. All groupsin the northern and western part of the NAP (GroupsA–D) have statistically indistinguishable parallax mea-surements. Of these, the value for Group D is the mostprecise, with (cid:36) ≈ .
21 mas (corresponding to ∼
795 pc).The parallax of Group F differs the most significantly,with a value of (cid:36) = 1 . ± .
02 mas, placing it behindthe other groups by ∼
130 pc. Group E has a mean par-allax of 1.27 ± ∼
35 pc, but this difference is onlysignificant at slightly over the 2 σ level, meaning that theparallax data are not absolutely conclusive. The modelalso suggests non-zero depths to the groups in paral-lax ranging from 0.04–0.12 mas. These seem unphysi-cally large when translated to physical sizes (25-80 pc)and may result from outliers with small error bars (Fig-ure 7). It is possible that some of these outliers could bemisclassified stars, either spurious members or memberswith a misassigned group.Uncertainty in our estimate of the absolute distancewould be dominated by the systematic ∼ ±
25 pc.Extinction and nebulosity seen in optical images pro-vide additional evidence to support the order in dis-tance suggested by parallax measurements. Group Fis superimposed on a bright region of the nebula, eventhough the stars in this group appear reddened on color-magnitude diagrams, suggesting that the H α nebulosityis in front of the group. The stars in Group E are associ-ated with cloud clumps in the “Gulf of Mexico” region,which are visible as dark dust lanes in the optical im-ages, suggesting that both the cloud and group are infront of the H ii region.The Bajamar Star is obscured by ∼
10 mag of cloudnear the periphery of the “Gulf of Mexico.” Thus, itis implausible that this star could be nearer than therest of the complex as its measured parallax indicates.More likely, the Bajamar Star is at a similar distance asthe complex, and its slightly larger measured parallaxis a ∼ ii region and the bright rim clouds. The “Gulfof Mexico” is superimposed between the Bajamar Starand the southeast bright rim cloud (“Pacific Coast ofMexico”), requiring the cloud to be in front of the lineof sight between this star and this rim.The order in distance of the other components is moredifficult to ascertain given the data. The cloud in the“Atlantic” region (possibly associated with Group D) isdark in the optical image, suggesting that it is mostlyin front of the H ii region. Group C sits atop the north-west bright rim (“Pelican’s Neck”). Thus, a clear line ofsight from the Bajamar Star to this bright rim, wouldput the “Atlantic Cloud” and Group D in front of the“Pelican’s Neck Cloud” and Group C. Both the opti-cally dark “Pelican’s Hat” cloud, associated with GroupB, as well as the cloud around Group A can be seen inabsorption, suggesting that they are also mostly in frontof the H ii region.Other group properties listed in Table 3 include aver-age proper motions, approximate centers of the groups,and the radii that contain half the groups’ stars. To cal-culate the proper motion, we use the same strategy asKuhn et al. (2019) and calculated the weighted medianwith the conventional 1 / error weights. To estimate thegroup centers, we calculate the sigma-clipped mean α and δ values, with a three standard deviation thresholdto avoid strong influences from outlying stars. STRUCTURE OF THE MOLECULAR CLOUDThe CO map provides kinematic information aboutthe system that is complementary to the stellar propermotions measured by
Gaia . Figure 8 shows the inte- tar and Gas Kinematics in NGC 7000/IC 5070 mem$phot_g_mean_mag[grp] m e m $pa r a ll a x [ g r p ] l ll lll l lllll ll l lllll l l ll lll l lllll ll l lllll l . . . . . Group A mem$phot_g_mean_mag[grp] m e m $pa r a ll a x [ g r p ] l l ll lll ll lll ll ll lll l l ll lll ll lll ll ll lll Group B mem$phot_g_mean_mag[grp] m e m $pa r a ll a x [ g r p ] l ll ll ll l ll lllllll ll ll l l lll ll ll ll lll ll l lll ll ll l l l ll ll ll l ll lll llll ll ll l l lll ll ll ll lll ll l lll ll ll l l Group C mem$phot_g_mean_mag[grp] m e m $pa r a ll a x [ g r p ] ll ll ll ll ll ll ll ll llll lll ll llll lll l ll l lll ll ll llllll l l llll lll llll lllll l ll ll ll ll llll ll lllll lll lll ll lll ll l l llllll l lll lll llll ll llll l lll lll ll l ll lll lll lll llll ll ll lll ll lll lll l ll l lll ll l ll lll l ll lll l llll lll lll ll ll ll l ll l l ll l ll l llll l ll lll l l llll ll ll ll ll ll ll ll ll ll ll l lll lll ll llll lll l ll l lll ll ll ll llll l l llll lll llll lllll l ll ll ll ll llll ll lllll lll lll ll lll ll l l llllll l lll lll llll ll llll l lll lll ll l ll lll lll ll l llll ll ll lll ll lll lll l ll l lll ll l ll lll l ll lll l l lll lll lll ll ll ll l ll l l ll l ll l llll l ll lll l l llll ll ll
10 12 14 16 18 20 . . . . . pa r a ll a x [ m a s ] Group D mem$phot_g_mean_mag[grp] m e m $pa r a ll a x [ g r p ] lll ll ll lll ll ll llll lll lll l ll ll llll l llllll ll ll llll l lllll l ll ll lll ll ll lll ll ll llll lll lll l ll ll llll l llll ll ll ll llll l lllll l ll ll
10 12 14 16 18 20
Gaia G [mag]Group E mem$phot_g_mean_mag[grp] m e m $pa r a ll a x [ g r p ] l llll llll ll l ll l llll llll ll l ll
10 12 14 16 18 20
Group F
Figure 7.
The panels show the mean parallax (colored horizontal line) and scatter plot of each star’s measured parallax vs.magnitude (points) for stars in each group. Uncertainty on the mean parallax is indicated by the colored band around the line,which shows the 2 σ confidence interval. The gray regions at the top and bottom of each graph indicate the range of parallaxesthat were excluded by the selection cut made when identifying members. grated CO emission over a velocity range of v lsr = − − , which encompasses nearly all the line emis-sion associated with the star-forming region. Dominantclouds are located in the “Gulf of Mexico,” “Pelican,”“Pelican hat,” and “Atlantic” regions. These areas tendto have gas at different velocities, with differences rang-ing from several to tens of km s − . Figure 9 shows inte-grated intensity maps and first-moment maps for threeadjoining velocity ranges. Clouds in the southern “Gulfof Mexico” region tend to have more positive velocities,while gas in the centrally located “Atlantic” region formsa ∼
15 pc long filamentary structure with more negativevelocities, and the eastern “Pelican” and “Pelican Hat”regions contains gas with both positive and negative ve-locities.An early radio study of the region by Bally & Scoville(1980) suggested that the complex is in the process ofexpansion. Zhang et al. (2014) also attribute the cloudmorphology in the “Pelican” region to an expandingshell around the H ii region. We find that some partsof the cloud have kinematics consistent with expansion,but that the expansion is not global. The shell-likestructure is most distinct at − − in the north-west. In contrast, in the south, the “Gulf of Mexico”cloud appears to be infalling, not outward-moving. The ∼
15 pc long filament at − − , which stretches from the north to the south, appears to be moving away fromthe H ii and towards us with a fairly coherent radial ve-locity across its whole length.To investigate the molecular cloud structure in moredetail, in Figure 10 we subdivide the system into indi-vidual clouds with contours at 6 K km s − from theintegrated map. This threshold picks out most of themain cloud components, which are labeled 1 through 12in both the figure and in Table 4. We also include aCloud 13 defined by contours at 3 K km s − that is as-sociated with stellar Group F, but does not meet the6 K km s − threshold.From the brightness of the CO J =1–0 emission line,it is possible estimate cloud masses if we make certainassumptions. To do so, we must accept that, as in mostCO studies, there are implied systematic effects on theresults. For example, the CO emission is less sen-sitive to diffuse gas traced by CO, but cannot tracedense gas where the CO line becomes too opticallythick. To calculate column densities, we assume localthermodynamic equilibrium (LTE), an excitation tem-perature, and the ratio of CO to H . Zhang et al.(2014) use a CO map from the Purple Mountain Ob-servatory Delingha 13.7 m telescope to estimate excita-tion temperature throughout the cloud complex. Theyfind temperatures ranging from 5–25 K, with the bulk of2
Kuhn et al. (cid:176) (cid:176) (cid:176) (cid:176) FWHM3 pc0 5 30K km/s (cid:176) (cid:176) (cid:176) (cid:176) FWHM3 pc−5 1 7km/s
Figure 8.
Left: integrated CO J = 1–0 emission over the velocity range v lsr = −
17 to 10 km s − . Right: The color scalefor the CO map uses intensity to indicate integrated emission and hue to indicate the first moment. The FWHM of theobservations and the physical scale of the map are indicated in the upper left. We also include the outline of the optical emissionfrom the DSS image (white contours). the cloud around ∼ ∼
14 K. How-ever, several areas have higher temperatures (20–25 K),particularly some edges of the clouds in the north of theNAP complex, including the cloud behind the bright rimin the “Pelican’s Neck” and a cloud filament in the “At-lantic” region. Consistent with the higher inferred tem-peratures, these regions appear geometrically orientedso they would be illuminated by the Bajamar Star. Wenote that the measured CO excitation temperatureswould be representative of the cloud surfaces, but thecloud interiors may not be the same.We follow the method outlined by Mangum & Shirley(2015) using the coefficients from the JPL MolecularSpectroscopy Database . Assuming a CO excitationtemperature of T ex = 14 K (as justified above), we es-timate the optical depth of the CO J =1–0 transitionwith the equation τ J =1 ( v ) = − ln (cid:40) − T , R ( v )5 .
29 K (cid:20) .
29 K /T ex ) − − . (cid:21) − (cid:41) , (2)using T R ∗ as an approximation for the radiation tem-perature T R in this equation. This equation is appliedto each pixel, channel-by-channel, in the emission-linedata cube. From the distribution of derived opticaldepths, ∼
95% of the CO gas by mass has opticaldepths τ < .
5, while ∼
98% has optical depth τ < https://spec.jpl.nasa.gov Although CO is not sensitive to the highest densityparts of the cloud, and thus cannot account for mass inthese regions, the shape of the distribution suggests thatmuch of the cloud is only moderately optically thick.The CO column density would then be N tot ( CO) ≈ (cid:0) . × cm (cid:1) (cid:18) T ex .
64 K + 13 (cid:19) × (cid:20) − exp (cid:18) − .
29 K T ex (cid:19)(cid:21) − (cid:82) τ J =1 ( v )d v km s − . (3)Finally, we assume that N (H ) ≈ . × N ( CO)(Bolatto et al. 2013) to calculate N (H ) column density.Integrating column density over the entire FCRAO CO map yields a total mass of 7 . × M (cid:12) . Massesof individual clouds are given in Table 4. Uncertaintyin these quantities is likely to be dominated by system-atic errors in the assumptions, particularly the CO-to-H ratio (Bolatto et al. 2013). In addition, clumpi-ness in the clouds that is not resolved by the FCRAObeam could cause us to systematically underestimate theamount of gas in regions with high optical depth.6.1. Relation of Clouds to Stellar Groups
There is spatial correspondence between the stellargroups and the clouds (Figure 10). Groups A and F areassociated with the minor clouds 1 and 13, respectively.Group E is embedded in the “Gulf of Mexico”, whichconsists of a major cloud component (Cloud 9) and sev-eral minor cloud components (Clouds 6, 10, 11, and 13)in the CO map. Group C lies on a bright rim cloud atthe southeastern edge of Cloud 2. The stars in Group Bare distributed in the northern part of Cloud 2 as wellas Cloud 3. tar and Gas Kinematics in NGC 7000/IC 5070 −10.9 to −3.8 km/s HPBW3 pc y −10.9 to −3.8 km/s HPBW3 pc0 5 16K km/s −10 −7 −4km/s −3.8 to 2.4 km/s
HPBW3 pc y −3.8 to 2.4 km/s HPBW3 pc0 6 25K km/s −3.0 −0.5 2.0km/s
HPBW3 pc y HPBW3 pc0 4 16K km/s 3 5 7km/s
Figure 9.
Integrated intensity maps (left panels) and first moment maps (right panels) over three velocity ranges, from top tobottom in order of increasing velocity. These three ranges highlight different aspects of the cloud’s structure. The top panel( − . − . − ) reveals a ∼
15 pc long filamentary structure, the middle panel ( − . − ) reveals structuresspanning the whole region, including a vaguely shell-like structure in the north, and the bottom panel (2.4 to 10.2 km s − )reveals a large cloud complex concentrated in the “Gulf or Mexico” region. Kuhn et al.
Figure 10.
Map of star positions and velocities, indicated by arrows, superimposed on the CO map. Arrows are drawnrelative to the median proper motions of the system, which we define to be ( µ α (cid:63) , µ δ ) = ( − . , − .
2) mas yr − . Arrows arecolor coded by the stellar groups A–F, and the outlines of the molecular clouds 1–13 are also shown. The two O stars areindicated by the cyan asterisks. Group D, the most spatially extended stellar group,has the most complicated relationship to the COclouds. This group is projected over a region that spansseveral clouds. However, Cloud 8 is the dominant cloudin this part of the NAP region and lies near the centerof Group D. The stars selected in our
Gaia analysis donot lie on top of Cloud 8, but instead the densest partsof Group D are next to the cloud. This indicates thatGroup D is either partially embedded in Cloud 8, whichwould imply that Cloud 8 is the remainder of the na-tal cloud that formed this group, or that Group D liesbehind the cloud. Stars from this group are also su-perimposed on other nearby minor clouds (Clouds 5, 6, and 7) which may either be smaller remnants of the na-tal cloud or chance superpositions. The edges of GroupD also intersect the edges of Clouds 2, 3, and even (inthe case of the Bajamar Star) Cloud 9.If the stellar groups and clouds are associated, then itis likely that they have similar mean motions. This en-ables us to associate mean proper motions with cloudsand radial velocities with stars, creating a three dimen-sional picture of the velocities in this region.6.2.
Velocity Structure of Clouds
Line profiles for each of the clouds (Figure 11) wereconstructed by spatially integrating the CO gas de- tar and Gas Kinematics in NGC 7000/IC 5070 v FWHM ) may be more characteristic of the kinematicsof the clouds because these statistics are less suscepti-ble to merging components with very different velocitiesthat may be from distinct clouds. Thus, we base ourdiscussion of cloud dynamics (Section 6.3) on these lat-ter quantities. We also estimate the characteristic size ofeach cloud, defined as the projected radius that containshalf the integrated emission.Optical depth can affect line profiles, flattening thepeaks and effectively broadening the lines (Hacar et al.2016). To more accurately resolve the velocity structure,we integrated over the optical depths given by Equa-tion 2 rather than T R ∗ (see Goldsmith & Langer 1999).This results in line profiles that are more sharply peakedwith FWHM that are ∼
20% narrower.The interpretation of the multimodal structure be-comes clearer with the position–velocity diagrams inFigure 12. Four of the major clouds (Clouds 2, 3, 8and 9) are shown as examples, while the diagrams forthe other clouds are accessible via the online figure set.Even some of the minor clouds (e.g., Clouds 5, 6, and7) have complicated velocity structures.In Cloud 2, v lsr increases by several km s − toward thesouthern end of the cloud where the optically bright rimand stellar group are located. It is plausible that thiscould be an effect of the expansion of the H ii regionpushing 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 separatedby ∼ − , while in Cloud 8 there are 3 components,also separated by ∼ − from each other. The dis-tinct components in Cloud 3 are also mostly separatedalong a north–south axis, possibly indicating that thiscloud is composed of multiple sub-clouds. However, forCloud 8 the main component is the dense northern endof the filament at v lsr ≈ − − (Figure 9), while theother components are related to overlapping filamentarystructures at higher v lsr . In our analysis, we tentativelyassociate the stars in Group D with the − − gasbecause this is the dominant component around whichthe stars are conglomerated, but our data are insufficientto make a definitive conclusion (Appendix C).The position-velocity diagram for Cloud 9 in Figure 12shows that it has a clumpy structure, with clumps span-ning the velocity range from ∼ − . A smallclump at − − is probably not part of this cloud,but is instead a continuation of the long north–southfilament seen in Figure 9 (top panels). Overall, Cloud 9has the largest positive v lsr of the complex, which, incombination with our evidence that the “Gulf of Mex- ico” and Group E are in front of the rest of the NAPregion, suggests that Cloud 9 is infalling. However, thepart of the cloud brightest in CO emission also has afairly modest velocity of ∼ − . This CO brightpatch is fairly near the Bajamar Star and shows up asa having a higher temperature in the CO map fromZhang et al. (2014).6.3.
Dynamical Properties of the Clouds
From the kinematic measurements of the cloud in Ta-ble 4, several dynamical quantities can be calculatedthat are useful for understanding star formation in thisregion. For these calculations, we use the line-widths forthe most prominent peaks in Figure 12 to avoid com-bining unrelated velocity components. The dynamicalquantities are given in Table 5.If we assume that the kinetic temperature T kin ofthe gas equals the excitation temperature measured for CO ( ∼
14 K) and that the mean molecular mass is¯ m = 2 .
33 amu, the one-dimensional sound speed in themolecular clouds would be c s ≈ .
22 km s − . For eachof the molecular clouds delimited in our analysis, theline width of the tallest mode is larger than can be ac-counted for by purely thermal broadening. The non-thermal component is σ nt = (cid:115)(cid:18) ∆ v FWHM . (cid:19) − kT kin m obs , (4)where m obs = 29 amu is the mass of the observedmolecule, CO. The Mach number is defined to be M = σ nt /c s . For most of the individual clouds, Machnumbers span 2 (cid:46) M (cid:46)
5. Exceptions are Clouds 8 and9 – the most massive clouds – which have larger velocitydispersions of
M ≈
M ≈
M ∼
6; Orkisz et al. 2017) or Taurus ( M < − − component appears curved inposition-velocity space, while, in Cloud 9, the large ve-locity dispersion is composed of multiple clumps withdifferent velocities.Several dynamically important quantities that can beestimated for each cloud include the mean density ¯ ρ ,free-fall timescale τ ff , and crossing timescale τ cross de-fined as ¯ ρ = M CO / (4 / π ) r (5) τ ff = (cid:112) π/ G ¯ ρ (6) τ cross = r cloud /σ nt . (7)6 Kuhn et al. −15 −10 −5 0 5 10 . . . . . . no r m a li z ed f l u x North and West
123 −15 −10 −5 0 5 10 v lsr [ km s − ] Center
South and East
Figure 11.
Spatially integrated line flux as a function of v lsr for each of the clouds in Figure 10. The three panels show velocityprofiles for clouds in three different areas of the NAP region. Molecular gas spans the velocity range −
15 to 10 km s − , andseveral of the clouds have multimodal structure. − − − v l s r [ k m s - ] − − − v l s r [ k m s - ] − − − v l s r [ k m s - ] − − − v l s r [ k m s - ] Figure 12.
Position–velocity diagrams of CO J = 1–0 emission are shown for Clouds 2, 3, 8, and 9 (upper left, upper right,lower left, and lower right, respectively).(The complete figure set (13 images) is available.) Both the crossing timescales and the free-fall timescalesrange from 0.1 to 2 Myr. For the main clouds withongoing star formation, the typical free-fall timescale is τ ff ∼ ∼ α = 2 T / |W| , where T is kinetic energy and W is potential energy. Previ-ous studies typically use the definition from Bertoldi &McKee (1992) that assumes a spherical, uniform densitycloud, giving α BM92 = 5 σ nt r cloud GM cloud . (8) tar and Gas Kinematics in NGC 7000/IC 5070 Figure 13.
Dynamical cloud masses versus mass estimatesfrom CO column density. The dashed line indicates equiv-alence between these quantities and the points are labeled bythe cloud they represent. Systematic uncertainties on CObased mass estimates may be a factor of several as discussedin the text.
However, the true value of α depends on the three-dimensional distribution of mass within the cloud (Singhet al. 2019; Guszejnov et al. 2020), with correctionsfor centrally concentrated mass distributions and fila-mentary structure typically being a factor of several(Bertoldi & McKee 1992). If we let α BM92 = 1, weobtain a dynamical mass estimate M dyn = 5 r cloud σ nt /G. (9)In Figure 13, we plot log M dyn versus log M CO forthe 13 clouds we have identified in the NAP region.These quantities are strongly correlated ( p < .
001 fromKendall’s τ test), with an approximately proportionalrelationship over a dynamic range of ∼
200 in cloudmass. However, the dynamical masses are systemati-cally ∼ ∼ ∼ CO-based mass estimates in Table 4. V´azquez-Semadeniet al. (2019) point out that, observationally, it is nearlyimpossible to distinguish between the free-fall and virialvelocities, which only differ by a factor of 2 in mass or afactor of √ M dyn , cal-culated from the velocity dispersions using Equation 9,and α BM92 , calculated from both the velocity dispersionsand the integrated CO cloud masses using Equation 8. STELLAR KINEMATICSFor the analysis of stellar kinematics, we convert theobserved astrometric quantities of parallax and propermotion into physical velocities with units of km s − . Forthis we establish a Cartesian coordinate system with po-sitions ( x, y, z ) and velocities ( v x , v y , v z ), where we se-lect the origin to be located at the center of the systemthat we are analyzing and moving with the median ve-locity of the stars in the system. Given the nature ofthe data, we are most interested in examining kinemat-ics in the two-dimensional xy plane. We define x to beparallel to α at the origin and y to be parallel to δ atthe origin. The transformations, including orthographicprojections, corrections for perspective expansion , andshifts to the rest frame of the system are described byKuhn et al. (2019, their Equations 1–4).We are also interested in the component of a star’sprojected velocity in the direction outward/inward rel-ative to the center of the group ( v out ) as well as theperpendicular azimuthal component ( v az ). These aredefined by v out = v · ˆr (10) v az = v · ˆ ϕ , (11)where ˆr and ˆ ϕ are the radial and azimuthal unit vectorsrelative to the group center in the xy plane. Uncertain-ties in these velocities are calculated by propagation ofthe Gaia astrometric errors (see procedure in Kuhn et al.2019). 7.1.
Global Kinematics
In Figure 10, we indicate the motions of NAP mem-bers using arrows, which, as in the previous figures, arecolor-coded by stellar group. The figure shows that starsnear each other tend to have similar velocities, but thedifferent groups appear to have fairly random motionsrelative to one another.Figure 14 (upper left) shows the relative 3D motions ofthe group centers inferred from both the mean propermotions of stars and the radial velocities of the asso-ciated clouds. This diagram shows no clear pattern ofeither convergence or divergence. Stars in Group A tendto be moving mostly north and slightly west – in a di-rection tangential to the NAP system center. Stars inboth Groups B and C are moving southeast, inward to-ward the center of the system. Group E is drifting slowlyeast, away from the center of the NAP complex in the xy plane, but it is likely moving rapidly toward the systemcenter in the third dimension. And, stars in Group F aremoving west – toward the center of the region. Finally, To compute the perspective expansion correction we assume thatthe stellar group is moving with the radial velocity of the mostsignificant of the cloud components associated with each group(Table 4) and transformed to heliocentric coordinates. Kuhn et al. −10 0 10 20 − ∆ x [ pc ] ∆ y [ p c ] group centers +1.0 km/s−1.9 km/s−1.1 km/s−5.3 km/s+4.0 km/s−35 pc+4.2 km/s+100 pc 10 km/s −15 −10 −5 0 5 10 . . . . . . v lsr [ km s − ] no r m a li z ed f l u x all gas σ all = km s − −10 −5 0 5 10 15 v x [ km s − ] c oun t s all stars σ = km s − −10 −5 0 5 10 15 v y [ km s − ] c oun t s all stars Figure 14.
Diagrams illustrating global motions in theNAP region. Upper left: Relative 3D motions of the groupcenters (reference frame µ α (cid:63) , = − . − , µ δ, = − . − ). Velocities in the plane of the sky from Gaia are shown by the arrows, while velocities from themolecular gas along the line of sight ( v lsr ) and offsets indistance are indicated by the labels. Upper right: Spatiallyintegrated CO J = 0–1 line flux for the whole NAP com-plex. The velocity dispersion is 4.2 km s − . Lower pan-els: Velocity distributions ( v x and v y ) for the entire NAPstellar population. The total velocity dispersion for stars is σ D = 2 . ± . − . The scale in km s − is the same forall three plots of velocity distributions. the center of Group D appears to be nearly stationaryin the xy plane in this reference frame. However, mostof the individual stars in this group are directed awayfrom its own center as can be seen in Figure 10.The Bajamar Star is not stationary relative to thecenter of mass reference frame, but is instead rapidlytraveling southeast away from Group D and towardsGroup E. The other O star, HD 199579, is travelingnortheast away from the NAP complex. In this globalreference frame , both stars have velocities > − .Such speeds would classify both objects as “walkaways”(Renzo et al. 2019; Schoettler et al. 2020), although nei-ther object has traveled outside the star-forming region.7.2. Kinematics within Groups
In Kuhn et al. (2019), we presented several methodsfor characterizing stellar kinematics and testing for ex-pansion, contraction, or rotation of a stellar system. In Section 7.3, we recalculate the O stars’ velocities with respectto Group D.
Here, we apply these methods to each of the stellargroups.To characterize total velocity dispersion in a group, wefit the ( v x , v y ) distribution with a bivariate Gaussian,taking into account measurement uncertainties, whichhave the effect of artificially broadening the velocity dis-persion. The stellar velocities and the Gaussian fit areshown in Figure 15 for each of the six groups. We sum-marize these distributions using the quantity σ D (Ta-ble 6), defined as the square root of half the trace of theGaussian’s covariance matrix (Kuhn et al. 2019, theirEquation 17). For the six groups, σ D ranges from 1 to2 km s − , being smallest for Groups B and F and largestfor D and E. Group F is a small group in a small cloud,so it is unsurprising that this group has a small veloc-ity dispersion. In contrast, Groups D and E are largegroups, located in a massive cloud (E) or in a partiallydispersed cloud (D), so it is unsurprising that they havelarger velocity dispersions. However, it is notable thatGroup C, spatially adjacent to Group D, has a muchlower velocity dispersion. In many cases, the velocitydispersions show signs of being anisotropic, similar tothe results for other star-forming regions (Kuhn et al.2019).We use several plots, analogous to those in Kuhn et al.(2019), to investigate evidence for expansion. Figure 16shows Group D, the largest group in our sample whereexpansion is most evident. In this group, the fastestmoving stars, with velocities > − are located nearthe edges of the region, and they are generally directedaway from the center of the group. The upper right plotuses both arrow direction and hue to indicate the direc-tion of motion, with color saturation indicating statis-tical uncertainty on direction. This use of color allowspatterns of motions to be identifiable as a gradient inthe color of the marks, in this case showing that stars atdifferent positions in the group are oriented in differentdirections, all preferentially away from the center.The two lower panels show statistical tests for expan-sion. The plot on the lower left shows the distributionof v out values; if a system is expanding then these willbe predominantly positive. For Group D, the weightedmedian v out value is 1.9 ± − , more than threestandard deviations above 0, clearly demonstrating thatthe group is expanding. This expansion velocity is largerthan any of those measured for the 28 systems investi-gated by Kuhn et al. (2019). The plot on the lower rightshows the mean expansion velocities for bins at differentradii from the group center. Expansion velocity appears The classification of sources using proper motion in Section 3 re-stricts the allowable range of velocities, which truncates the tailsof the velocity distributions and will slightly decrease the esti-mated σ D . We simulate this effect and find that it will producea <
8% effect on the results for velocity dispersions < − ,which is smaller than the estimated statistical uncertainty on thederived values. tar and Gas Kinematics in NGC 7000/IC 5070 Group A − Group B Group CGroup D v y [ k m s − ] − −5 0 5 Group E v x [ km s − ] −5 0 5 Group F−5 0 5 Figure 15.
Scatter plots of star velocities. The magenta ellipses indicate the best-fit Gaussians; the ellipses are drawn atMahalanobis distances of 1 and the shaded region shows the 95% confidence intervals on the ellipses. to increases with distance from the center, a trend thatis seen in many of the expanding young star clusters andassociations from Kuhn et al. (2019).The other groups contain fewer stars than Group D,and the clear evidence for expansion is not found usingthese diagrams. For the other groups, arrow diagramslike those in the top two panels of Figure 16 are in-cluded in the figure set, but there are too few stars tomake the other plots. For these groups, the weightedmedian v out values are not significantly different from 0(Table 6) owing to the large statistical scatter. Never-theless, median v out may not be the most sensitive testfor expansion, given that this statistic combines multipledimensions and does not take advantage of the tendency,seen in many expanding regions, for velocity to increasewith distance from the center.We apply a further test for expansion, using the non-parametric Kendall’s τ test to test for correlation be-tween x and v x and between y and v y (Table 6). Statis-tically significant results would indicate a velocity gra-dient, which would imply expansion if the gradient ispositive and contraction if negative. We find strong sta-tistical evidence for expansion in both dimensions forGroup D, moderate statistical evidence for expansionin both dimensions for Group B, and strong statisticalevidence for expansion in only one dimension ( x ) forGroups C and E.7.3. Velocity Gradients in Expanding Groups
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 thedistance from the groups’ centers (Kuhn et al. 2019;Rom´an-Z´u˜niga et al. 2019; Zamora-Avil´es et al. 2019;Wright et al. 2019; Melnik & Dambis 2020). However,Wright et al. (2019) point out that the expansion may beanisotropic. In Appendix D, we describe a linear modelfor velocity as a function of position, using a methodthat can account for anisotropy and is independent ofchoice of coordinate system. Table 6 reports severalof the interesting properties from the models, includ-ing velocity gradients in the directions of maximum andminimum expansion, the position angle of anisotropy,and the intrinsic scatter that remains after the velocitygradient is accounted for.We find statistically significant expansion forGroups C, D, and E, but no group had statisticallysignificant rotation. Group D has the strongest relationbetween position and proper motion, with expansion inall directions, but with a larger gradient in the east-westdirection (0.5 km s − pc − ) than in the north-south di-rection (0.3 km s − pc − ). In contrast, Groups C and Ehad non-zero expansion gradients only in the east-westdirection, but not the north-south direction. We notethat the anisotropy position angles of all three groupsare similar. The regression analysis and the calculationof the model parameters are described in greater detailin the appendix.Figure 17 shows scatter plots of position versus ve-locity with the model regression lines overplotted. Thegraphs of v x vs. x (upper left) and v y vs. y (upper right)exhibit positive correlations if there is expansion, while0 Kuhn et al. −10 −5 0 5 10 15 − − − x [pc] y [ p c ] +43 30+44 00+44 30+45 00 485220 56 ll ll lll ll ll l ll l lll l lll ll ll l l ll ll ll l ll llll ll lll llll lll l ll lll ll l lllll l l ll l ll l lll lll l lll ll lll l llll lll ll lllll ll ll l ll lll ll lll ll ll lll llll l lll lll l ll l ll llll ll llll ll ll ll lll ll lll ll l lll l ll lll l lll ll l lll ll l ll ll ll lll ll lll l lll l lll l ll ll lll l lll ll lll ll Group D −10 −5 0 5 10 − − − x [pc] y [ p c ] +43 30+44 00+44 30 5220 56 Group D NE S W−10 −5 0 5 10 . . . . . . v out [ km s - ] s m oo t hed d i s t r i bu t i on Group D 0 2 4 6 8 10 12 14 radius [pc] m ed i an v ou t [ k m s - ] l l l l l l l l Group D y = x(0.5±0.1) − 0.4±0.7
Figure 16.
Plots of stellar kinematics for stars in Group D. Top left: Vectors show motions of stars in the ( x, y ) plane relativeto the center of the group. Stars with the most precise velocity estimates ( < − ) are shown in black, while stars with lesscertain velocity estimates are shown in gray. Top right: Arrow heads are color coded by the direction of motion of a star asindicated by the color wheel. Bulk motions such as expansion, contraction, or rotation can be seen as gradients in color acrossthe group. Bottom left: Distribution of the outward velocity component v out . The median value is indicated by the magentaline and the 3 σ confidence interval is shown by the shaded area. Bottom right: Median value of v out as a function of distancefrom the center of the region. An error weighted least-squares regression fit is shown.(The complete figure set (14 images) is available.) tar and Gas Kinematics in NGC 7000/IC 5070 −10 −5 0 5 10 − x [pc] v x [ k m s − ] expansionplot −15 −10 −5 0 5 − y [pc] v y [ k m s − ] expansionplot −10 −5 0 5 10 − x [pc] v y [ k m s − ] rotationplot −15 −10 −5 0 5 − y [pc] v x [ k m s − ] rotationplot Figure 17.
Plots of velocity coordinates v x and v y versus position coordinates x and y for stars in Group D. The top rowshows v x vs. x and v y vs. y , so expansion/contraction would produce correlation on these diagrams. The bottom row shows thecombinations v y vs. x and v x vs. y , so correlation here would indicate rotation. We have over-plotted graphical illustrations of ourMCMC model fit for comparison to the data. On the top two plots we show lines indicating the relations v x = v x, + xA (left)and v y = v y, + yA (median and 95% credible region shown in magenta), while the bottom two plots show v y = v y, + xA (left) and v x = v x, + yA (right). The contribution from the intrinsic velocity dispersion, described by the covariance matrixΣ scatter , is indicated by the red bars.(The complete figure set (3 images) is available.) v y vs. x (lower left) and v x vs. y (lower right) exhibitcorrelations if there is rotation. The regression line and95% credible interval from the linear model are shown inmagenta. The intrinsic velocity scatter, applied equallyto all points, is indicated by the red error bar (1 σ ), whilethe individual 1 σ measurement uncertainties are indi-cated by the gray lines. For Group D, both the v x vs. x and v y vs. y plots show strong correlations, while forGroups C and E correlations are only strong on the v x vs. x plots due to their θ ≈ ◦ anisotropy position an-gles.The O stars are both among the fastest moving mem-bers of Group D. The Bajamar Star has a velocity of6.6 ± − relative to the center of the group, while HD 199579 has a velocity of 6.4 ± − . However,neither of these velocities are discrepant from the lin-ear model. This suggests that whatever phenomenonaccelerated these stars to their current velocities is alsoresponsible for the expansion of the group as a whole. Ifwe make the simplifying assumption that the velocitieshave been approximately constant, the Bajamar Starwould have been nearest, in projection, to the center ofGroup D ∼ ∼ Effects of Spatial Covariances in GaiaAstrometric Errors Kuhn et al.
To ensure our measurements of velocity gradients arerobust, we examine how these would be affected by thesystematic correlated errors in astrometry (Lindegrenet al. 2018, their Section 14). They found that the co-variance of quasar proper motion as a function of angu-lar separation decreases from ∼ µ as yr − at 0.07 ◦ separation to a local minimum of ∼ µ as yr − at 0.43 ◦ separation. This could induce an artificial proper mo-tion gradient of ∼ − deg − , which would looklike a velocity gradient of 0.045 km s − pc − in a stel-lar association at a distance of 795 pc. Given that thevelocity gradients that we detect are significantly largerthan this value, they are most likely astrophysical. NEW MEMBER CANDIDATES FROM
GAIA
In addition to validating previously proposed YSOsin the NAP region,
Gaia astrometry can be used tosuggest new member candidates. In most earlier stud-ies, candidates were identified based on criteria that re-quired a star to possess a circumstellar disk (e.g., strongH α emission, infrared excess, large amplitude variabil-ity). However, many star-forming regions are dominatedby pre–main-sequence stars for which disks are not de-tected (Broos et al. 2013). Thus, astrometric selection ofcandidates can reveal whether any populations of NAPmembers have remained undetected in previous analy-sis, which, in addition to their possible biases, are alsoincomplete in their selection of candidate members.We start with a sample of all Gaia
DR2 sources within3 ◦ of the Bajamar Star that pass the same quality cri-teria from Section 2.2, and lie within the same par-allax range and the ( µ (cid:96) (cid:63) , µ b ) region identified in Sec-tion 3. However, we apply further parallax cuts to se-lect only objects that are consistent within 2 standarddeviations of the median parallax of the groups, i.e. (cid:36) + 2 σ (cid:36) ≥ .
06 mas and (cid:36) − σ (cid:36) ≤ .
24 mas. Thesecriteria greatly reduce the number of
Gaia sources inthe region, from ∼ Gaia sources in the originalparallax range down to ∼ Gaia sources whose absolute G -band magnitude and G − RP color are two standard deviations above the 3 Myr Bres-san et al. (2012) isochrone with A V = 1 mag (Fig-ure 18, left). However, remaining contaminants can in-clude main-sequence stars with high extinction as wellas post-main-sequence stars. To reduce contaminationfrom reddened main-sequence stars, we use the 2MASS J − H vs. H − K s diagram to estimate extinction usingthe reddening laws from Rieke & Lebofsky (1985). Typ-ical reddening is ∆( J − H ) ∼ . ≈ . V band), but ranges from 0 to 0.5 mag, with a tail out to1.5 mag. Dereddening on the near-infrared color-colordiagram can lead to degeneracies between low-mass andintermediate/high-mass stars. In most cases we pick the low-mass solution. However, when a star could below-mass with ∆( J − H ) < . J − H ) ≥ . Gaia isochrone with thecorresponding reddening. To reduce contamination bygiants, we selectively remove objects from the regionof 2MASS color magnitude diagram where these objectsare likely to dominate the sample: J − K s > . M K s < . × ( J − K s ) − . k -nearest neighbor( k NN) algorithm implemented in the R package class (Venables & Ripley 2002). For the training examples, weused the 395 previously assigned members of groups A–F in addition to an equal number of randomly selected“distributed” stars from
Gaia in the same parallax andproper motion range that were not selected as new can-didates. The k NN classifier is run using the normalized (cid:96) , b , µ (cid:96) (cid:63) , and µ b variables, and classifications are basedon the 5 nearest neighbors. We assign 52% of the candi-dates to groups and 48% to the “distributed” category.A previously unidentified group of ∼
17 new membersis centered at 21:04:44 +43:55:30. We call these starsGroup G and list them in Table 7. This group is lo-cated next to the third-magnitude red-supergiant star ξ Cyg. However, this star has a parallax of ∼ tar and Gas Kinematics in NGC 7000/IC 5070 G − RP (Gaia) [mag] G + l og ( p l x [ m a s ] ) − A V = mag1 Myr3 Myr10 Myr 0.0 0.5 1.0 1.5 2.0 2.5 − J − Ks (2MASS) [mag] K s + l og ( p l x [ m a s ] ) − A V = magNew CandidatesValidated YSOsField Stars Figure 18.
Color magnitude diagrams in the
Gaia bands (left) and near-infrared (right) showing how the distributions of newcandidates NAP members (black points) compare to other categories of sources. The green shaded regions show the density of
Gaia sources with the same parallax range, with darker colors indicating higher number densities of sources. On both plots,field stars in the red clump and the giant branch can be seen near the top of the diagram. We use cuts on the near-infrareddiagram (indicated by dotted lines) to remove sources from our candidate sample that have a high probability of being giantstars. On the near-infrared diagram, we also show YSOs from our
Gaia validated literature sample (magenta points). Thedistribution of these points is shifted to slightly higher K s band luminosities and redder J − K s colors. This may be the resultof earlier studies selecting stars with disks, and thus more likely to have K s -band excess emission, while the Gaia selection isindependent of whether the object has a disk.
Group E, so we classify the new stars as a continuationof this group. DISCUSSIONObservations of the NAP complex reveal that it con-tains multiple groups of very young stars ( ∼ Gaia parallaxes suggestthat there are several conglomerations along the line ofsight spanning up to 150–200 pc. The motions of thegroups, as inferred from
Gaia astrometry, appear to beapproximately random relative to one another, and sev-eral of the groups are rapidly expanding.9.1.
Expanding Stellar Groups
Several mechanisms, including disruption by tidalforces, cloud dispersal, and sequential star formation,may play roles in the expansion of stellar groups in theNAP complex.A clumpy mass distribution, as observed here, willgenerate internal tidal forces that affect stellar dynam-ics and influence the potential assembly of either boundclusters or unbound associations. Tidal fields would naturally produce an anisotropic velocity gradient sim-ilar to those observed in the NAP region. Furthermore,tidal forces would preserve clumpy substructure withina group (e.g., D and E) as it is stretched, while processesdriven by dynamical relaxation would be more likely toerase these substructures. Tidal forces from giant molec-ular clouds have long been known to be a major cause ofstar-cluster disruption (Spitzer 1958; Gieles et al. 2006),and these effects can occur in the star-forming environ-ment itself, inhibiting the initial formation of a boundsystem (Kruijssen et al. 2011, 2012). Zamora-Avil´eset al. (2019) have run numerical simulations in whichgas clouds, driven away from the newly formed youngstellar cluster by feedback, exert tidal forces on the clus-ter that pull it apart.The tidal acceleration a t, axial per unit of displacement∆ r along the separation axis due to a mass M at dis-tance r is a t, axial / ∆ r = 2 G Mr . (12)For Group D, the magnitude of the tidal accelera-tion from Cloud 2 (2 . × M (cid:12) ) located ∼ Kuhn et al.
10 pc
Gaia candidate member
21 10 21 00 20 50 20 40424446
Figure 19.
Spatial distributions of new candidate members. Left: All new candidates within 3 ◦ of the Bajamar Star. Thedashed box outlines the blown-up central region shown in the other panel. Right: New candidates are shown as +’s if they areassigned to existing groups or (cid:5) ’s otherwise, while the validated YSOs from previous studies are shown as circles. Color-coding ofgroup members is the same as in previous figures. This plot shows that most new candidates in the central region are distributedsimilarly to the distribution of the previously identified stars. − pc − Myr − . This acceleration would beenough to induce the observed ∼ − pc − east–west velocity gradient in ∼ ∼
35 pc along the line of sight means that itwould much less strongly affected by tidal forces fromother clouds. Thus, any tidal force would need to be gen-erated by clumpy distributions of mass within Cloud 9itself.The ongoing dispersal of the molecular clouds due toO star winds and photoevaporation, combined with out-flows from low-mass stars (Bally & Reipurth 2003; Ballyet al. 2014), will weaken the binding energy of the stellargroups, causing them to expand and/or disperse (e.g.,Tutukov 1978; Adams 2000; Kroupa et al. 2001, andmany others). In the northwest, the shell-like structuremay be the outer rim of a bubble blown by the BajamarStar, which would have been closer to this part of the cloud ∼ σ D = 2 km s − and as-suming an initial size of ∼ ∼ × M (cid:12) . Thus,the natal cloud mass would have been ∼ ∼ tar and Gas Kinematics in NGC 7000/IC 5070 ii regions (e.g., Patel et al.1998) the stars formed in this way would presumablyinherit the velocities of the material in the expandingshells. In the NAP region, there is ongoing star for-mation in several of the clouds making up the shell-structure, so it is plausible that such a scenario could behappening here too. In particular, the Group C is asso-ciated with the bright-rim cloud at the edge of Cloud 2.For several bright rim clouds in other star-forming re-gions, Getman et al. (2007, 2012) have shown that thestars in front of the cloud exhibit an age gradient, withthe youngest stars nearest the cloud. In such a scenario,as a cloud is driven outward by an expanding shell, thesequential formation of stars could yield a velocity gra-dient.Given the conditions in the NAP region, it seems likelythat multiple scenarios operate simultaneously to inducethe observed expansion patterns. For example, the less-ening of the binding energy due to cloud dispersal wouldmake groups more susceptible to disruption by inter-actions with neighboring clouds. Although theoreticalwork has shown that binary-star dynamics can prefer-entially eject O stars at high velocities ( >
30 km s − ; Oh& Kroupa 2016), this mechanism does not appear to berequired here because the O stars and low-mass starsfollow the same velocity pattern.The Gaia data show that individual groups are ex-panding, but the relative motions of the groups exhibitno clear pattern of either convergence or divergence (Fig-ure 14, upper left). This is similar to the situation foundby Kuhn et al. (2019) in other large star-forming com-plexes like NGC 2264, NGC 6357, or the Carina Nebula,and may help explain why expansion has been difficult todetect in young stellar associations analyzed as a whole(e.g., Ward & Kruijssen 2018; Wright et al. 2016; Wright& Mamajek 2018). The velocity dispersion of all stars inthe NAP complex is σ D = 2 . − (Figure 14, lowerpanels). This, combined with the separations of severalto tens of parsecs between groups, means that the massneeded to gravitationally bind the system ( > M (cid:12) )is orders of magnitude greater than the mass in stars.Thus, it is impossible for the system to coalesce as abound cluster. Nevertheless, it may be possible in somegroups for bound remnants to remain behind even aftermost stars have escaped (Baumgardt & Kroupa 2007).9.2. Star-Formation Scenarios
We can use our inferences about the dynamical stateof the stars and gas, as well as the constraints on stel-lar ages, to test various theoretical scenarios for starformation as applied to the NAP region. These sce-narios broadly break down into fast star formation, onthe timescale of the free-fall collapse of the molecularclouds (e.g., Elmegreen 2000; Hartmann et al. 2012),and slow star-formation that persists over multiple free- fall timescales (e.g., Tan et al. 2006; Krumholz & McKee2020). In the former case, molecular clouds would havelittle pressure support (V´azquez-Semadeni et al. 2019),while in the latter case turbulent pressure support sta-bilizes the clouds against collapse and allows for an ex-tended duration of star formation (Padoan & Nordlund2011).In the NAP region, we demonstrated that the veloc-ity dispersion within individual clouds is consistent withthe velocity dispersion that would be expected from ei-ther free-fall collapse or virial equilibrium given the sys-tematic uncertainties on the cloud masses (Section 6.3).When we consider the whole region (Figure 14), thecombined velocity dispersion of all gas is larger thanthat of the individual clouds, with σ all = 4 . − . Inthe complex, half the gas is contained within a projectedradius of 9.4 pc, so the dynamical mass for the whole sys-tem calculated using Equation 9 would be 1 . × M (cid:12) ,or about two and a half times the 7 . × M (cid:12) massestimated from CO column density. Estimation de-pends on the virial parameter α , which we take to be1 (virial equilibrium); however, a cloud with α = 2(free fall) would yield a dynamical mass estimate half aslarge, making the differences smaller. As earlier, the sys-tematic uncertainties make distinguishing between thesescenarios difficult. Nevertheless, in either case, the ve-locity differences between the clouds, as well as the ve-locity dispersions within the clouds, can be accountedfor by gravitation.The free-fall timescales for the clouds are ∼ ∼ × τ ff . There aretwo considerations that should be taken into accountwhen interpreting free-fall timescales. First, free-falltimescales are typically calculated using the mean den-sity of a cloud as we have done; a non uniform den-sity could change the mean τ ff , but the correction factoris expected to be of order unity (Krumholz & McKee2020). Second, we are calculating the present day free-fall timescales; however, in the past, when the stars ob-served today were forming, the free-fall timescale wouldlikely have been longer (V´azquez-Semadeni et al. 2019).Star cluster formation in a free-fall time is expected foreither clouds without turbulent support (Klessen et al.1998) or with decaying turbulence (Klessen 2000). Thus,the observation that most of the stars are only one toseveral τ ff old is consistent with the rapid star-formationscenarios.The spatial clustering of star formation also favors ascenario of rapid star formation in a turbulent molecu-lar cloud. Overall, the NAP cloud complex has featuresexpected for a giant molecular cloud sculpted by turbu-lence, including a clumpy, hierarchical density structureand velocity structure that yields a larger velocity dis-persion for the whole system than within the individualclouds, broadly consistent with the Larson (1981) re-6 Kuhn et al. lation. Given that the major stellar groups follow thepresent-day cloud structure, this supports the view thatthe stars froze out on a timescale less than the crossingtime for the complex (Elmegreen 2000).The similarity in age ( ∼ t s = 9 . /c s = 40 Myr, meaning that the clouds wouldnot be in thermal contact with each other on the col-lapse timescale, and we would not expect them to havesimilarly aged populations. A similar problem has beenposed by Preibisch & Zinnecker (1999) for the UpperScorpius OB association and by Herczeg et al. (2019)for the Serpens molecular clouds. A possible solutionto this is an accelerating star-formation rate, which isa feature of several models, including the conveyor beltmodel (Longmore et al. 2014) and global hierarchicalcollapse (V´azquez-Semadeni et al. 2017, 2019). Withsufficiently high acceleration, most of the stars in star-forming clouds would have formed very recently, even ifthe first stars to form in those clouds formed at differenttimes.Although the NAP clouds appears to be formingstars on free-fall timescales, fast star-formation scenar-ios face challenges if applied to the entire Galaxy be-cause they would imply a Galactic star-formation ratemuch higher than observed (Zuckerman & Evans 1974;Krumholz et al. 2006). To overcome this, theoreticalfast star-formation scenarios invoke the quick disrup-tion of the molecular clouds by stellar feedback to keepstar-formation inefficient and regulate the Galactic star-formation rate (Chevance et al. 2020). In the NAPregion, where the disruptive influences of the BajamarStar and outflows from low-mass stars are already ap-parent, it appears likely that feedback processes couldbring star formation to an end before a substantial frac-tion of the cloud mass is converted into stars. SUMMARYIn this study we have examined the kinematics of starsand gas in the NAP region using previously publishedcandidate YSOs,
Gaia astrometry, and a newly pre-sented molecular gas map. Our main results are thefollowing. • On the basis of
Gaia parallax and proper mo-tions, we identify a sample of 395 stars as high-probability members of the complex (estimated3% residual contamination). We also reclassifyhundreds of previously cataloged candidates ascontaminants. In addition, ∼ • The locations of the confirmed members onthe
Gaia color-magnitude diagram suggests thatnearly all stars are < ∼ • Most of the widely dispersed YSO candidates fromprevious studies are identified as contaminants,while the confirmed members tend to be moretightly clustered. We identify 6 groups using un-supervised cluster analysis in position and propermotion. • The NAP region is ∼
795 pc from the Sun. How-ever, parallax distributions show slight differencesin distances to the individual kinematic groups.Notably, the stars in the “Gulf of Mexico” regionappear to be ∼
35 pc closer than the rest of thesystem, and a small group, containing the famousV1057 Cyg star, is ∼
130 pc farther away. • The locations of each of the stellar groups are spa-tially correlated with the main components of themolecular cloud. This implies that all parts of theNAP cloud complex, if sufficiently massive, are ac-tively forming stars. • Most of the clouds have complex, multimodal ve-locity structures. We use the CO map to esti-mate various cloud properties, including massesand velocity dispersions. We find high correla-tion between the masses estimated from integrated CO column densities to be in remarkable agree-ment with dynamical mass estimates, suggesting astrong connection between gravitation and veloc-ity. Mean free-fall times for individual clouds are (cid:46) • Relative velocities of the different groups appearrandomly oriented, showing no sign of either globalexpansion or contraction. The radial velocity ofthe “Gulf of Mexico” region implies that it isplunging inward. On the other hand, in the north-west of the NAP complex (“Atlantic” and “Peli-can” regions) where the morphology of the CO gassuggests an expanding H ii region, the relative mo-tions of the stellar groups seem unaffected. • In contrast to the lack of global expansion, severalstellar groups are individually rapidly expanding,with velocity gradients of 0.3–0.5 km s − pc − .The expansion gradients are anisotropic, and weargue that these could be, in part, attributedto tidal forces from within the clumpy molecularcloud complex. tar and Gas Kinematics in NGC 7000/IC 5070 • The primary ionizing source in the region, theearly-O Bajamar Star, lies between two groups.We suggest that it likely originated as part of agroup (which we call Group D) centered in the“Atlantic” part of the NAP region because its tra-jectory would trace back to this group. The star’svelocity ( ∼ − ) is consistent with the expan-sion velocity seen for low mass stars of the samegroup. Another O-type star, HD 199579, appearsto be ejected in a different direction from the samegroup. • We identify > Gaia catalog, including anew seventh stellar group (Group G) located eastof the complex. Slightly over half these new can-didates are associated with Groups A–G, whilethe others are more broadly distributed through-out the 6 ◦ -diameter selection area of investiga-tion. These objects would require follow-up ob-servations to validate them.Several lines of evidence suggest a scenario of rapidstar formation in a free-fall time (e.g., Elmegreen 2000;Hartmann et al. 2012; V´azquez-Semadeni et al. 2019) in the NAP complex. The structure of the stellar groupsis spatially correlated with the structure of molecularclouds, stellar ages are similar to the free-fall timescalesof the star-forming clouds, and the gas velocity disper-sions in the clouds are consistent with the velocities ex-pected from gravitational collapse. Furthermore, an ac-celerating star-forming rate, as predicted by some of themodels, would mean that regardless of how long starshave been forming, the majority of stars would havebeen formed during the last few free-fall timescales. Thiscan explain why several distinct clouds in the NAP re-gion that are sufficiently far apart to be out of directthermal contact could all have very young stellar popu-lations of nearly the same age.Nevertheless, some important properties of the starforming region remain relatively poorly determined. Forexample, the census of stars in the region is still highlyincomplete even after our study, making it difficult toconstrain the total stellar mass. While, our sample of395 Gaia -validated YSOs is useful for understanding thespatial and kinematic distributions of stars in this re-gion, it is only representative of the optically brightesttail of the full stellar population.APPENDIX A. PARAMETERS OF THE NAP MEMBERSHIPCLASSIFIERBefore applying the final step of the NAP membershipclassifier (i.e. classification using proper motions, de-scribed by the equations below) we first remove sourcesthat either do not meet our
Gaia quality criteria or haveparallaxes outside the range 0.86–1.61 mas (Section 3).For YSO candidates that pass these steps, membershipprobability, p mem ( µ (cid:96) (cid:63) , µ b ) is calculated using the Gaus-sian mixture model (Figure 1, right panel) with the fol-lowing parameters: d mem ( µ i ) = 0 . φ ( µ i ; µ , Σ ) (A1) d field ( µ i ) = 0 . φ ( µ i ; µ , Σ ) (A2) p mem ( µ i ) = d mem ( µ i ) / [ d mem ( µ i ) + d field ( µ i )] (A3) µ = ( − . , − .
15) (A4)Σ = (cid:34) .
48 0 . .
06 0 . (cid:35) (A5) µ = ( − . , − .
79) (A6)Σ = (cid:34) . − . − . . (cid:35) , (A7)where the d ’s represents the density of stars in proper-motion parameter space, φ denotes the bivariate normaldistribution, µ = ( µ (cid:96) (cid:63) , µ b ) are proper motions in Galac- tic coordinates in units of mas yr − , and Σ i are thecovariance matrices of the Gaussian components. B. SELECTION EFFECTS RELATED TOSTELLAR AGESSelection of YSOs using features that are connectedto disks and accretion means that pre–main-sequencestars that have lost their disks will be missed. Thisimposes the disk survival function, often modeled asan exponentially decreasing function with an e -foldingtimescale of 2–4 Myr (Mamajek 2009; Ribas et al. 2015;Richert et al. 2018), as a bias on the observed age dis-tribution. In contrast, X-rays can identify pre–main-sequence stars both with and without disks (Feigelson &Montmerle 1999; Feigelson 2018). Low-mass pre–main-sequence stars maintain high X-ray emission even after10 Myr (e.g., Argiroffi et al. 2016; Preibisch et al. 2017),so X-ray selection should detect such a population if itexists.The study by Damiani et al. (2017) provides an X-ray sample for the NAP region. In Section 3, we usedthe positions of the NAP members on the Gaia color-magnitude diagram as evidence that the NAP membersare < G vs. G − RP dia-gram, the X-ray sources lies within the distribution ofstars selected by other methods. This suggests that bothsamples have similar age distributions. On the J vs.8 Kuhn et al. J − H diagram, relatively few X-ray sources have high J − H values, indicating that the X-ray selected sampleis not as highly reddened as the disk-selected sample.Nevertheless, nearly all X-ray sources are above or tothe right of the 1 Myr isochrone. Thus, we concludethat there is no evidence for an older pre–main-sequencestellar population within the nebula. C. STAR FORMATION WITHIN CLOUD 8Cloud 8 is superimposed near the center, and densestpart, of the expanding stellar group D. Given this con-figuration, it appears that Cloud 8 could be the mainremnant of the cloud responsible for producing this stel-lar group. However, it is alternatively plausible thatthis superposition is coincidental, and that the cloudis closer to us than the stellar group is. It is particu-larly difficult to distinguish between these two possibil-ities because our sample does not include many objectswithin Cloud 8, presumably due to difficulty detectingobscured YSOs. However, Cambr´esy et al. (2002) iden-tified a dense, embedded cluster (number 6 in their cata-log; hereafter Cambr´esy 6) near the center of this cloud.Figure 21 shows the
Spitzer /IRAC 3.6 µ m image alongwith our Gaia members (green circles),
Spitzer
YSOsfrom Guieu et al. (2009) and Rebull et al. (2011) notincluded in our sample (yellow squares), the contoursof Cloud 8 (outer white curves), and the location ofCambr´esy 6 (white circle). A star cluster can clearly beseen in the
Spitzer image at the location of Cambr´esy 6,but almost none the individual cluster members wereidentified as YSOs by any survey. This suggests thateither the stars did not have infrared excess or the mid-infrared nebulosity in the region prevented reliable de-tection of infrared excess. Nevertheless, we suspect thatCambr´esy 6 is young because, in
XMM Newton images(not shown), a groups of spatially confused X-ray pointsources can be seen at this location. The high absorp-tion of this group and its spatial coincidence with themolecular gas suggests that it lies within Cloud 8 andthat the cloud is actively forming stars. Nevertheless,the relation between Cambr´esy 6 and Group D remainsa matter for future studies. D. MODELING LINEAR EXPANSIONHere we describe a linear model to relate the positionof a star to its velocity. This model can account forexpansion (homologous or anisotropic), contraction, andeven some rotational effects.We model the velocity of the i -th star, v i = ( v x,i , v y,i ),as being related to its position, x i = ( x i , y i ), by thelinear regression equation v i = v + M x i + v scatter ,i , (D8)where v is a constant velocity shift, M is a 2 × v scatter ,i is a random vector drawn froma bivariate normal distribution with covariance matrix Σ scatter . We use v scatter to represent the intrinsic ran-dom deviations in velocity. In addition to the intrinsicscatter in velocity, in the observed velocities, v obs ,i , wemust also account for Gaia ’s measurement uncertainty.Thus, v obs ,i is related to its actual velocity by the equa-tion v obs ,i = v i + ε i , (D9)where measurement errors ε i are drawn from Gaus-sian distributions with covariance matrices Σ ε,i obtainedfrom the Gaia
DR2 catalog. Putting this together giveus the likelihood equation, v obs ,i ∼ N ( v + M x i , Σ scatter + Σ ε,i ) , (D10)from which we estimate the model parameters using aBayesian approach.For the Bayesian model fitting, we use Markov chainMonte Carlo (MCMC) to sample from the posterior dis-tribution, which is implemented by the software pack-age “Just another Gibbs sampler” (JAGS) version 4.3.0(Plummer 2017), which is run using rjags (Plummer2019) in R version 3.6.0. We use “non-informative” pri-ors for the parameters, including uniform distributionsfor v ,x and v ,y from − − , and broad dis-tributions for M and Σ scatter that cover the reasonablevalues for these parameters. The statistical model wasimplemented in the BUGS language. For each dataset,three independent chains were run for 5,000 iterations,retaining every fifth sample, after an initial burn-in of1,000 iterations. Convergence was assessed from inspec-tion of the trace plots as well as the Gelman-Rubinstatistic < .
001 (Gelman & Rubin 1992). Through ex-perimentation with a variety of priors, we found thatchanges to the functional form for the priors has lit-tle effect on the results, provided that the priors aresufficiently broad. We also find that the results of theBayesian analysis are approximately consistent with theresults from numerical maximum-likelihood estimation.The 2 × M describes the dependence of veloc-ity on position, encoding phenomena such as expansion,contraction, or rotation, which may be either radiallysymmetric or anisotropic. To examine these effects, wedecompose M using singular value decomposition (SVD)to write M = U DV ∗ , (D11)where U and V ∗ are unitary matrices and D is a diagonalmatrix. By multiplying by the identify matrix, writtenas the product of two reflection matrices, without loss The prior distribution for the matrix Σ scatter is made up of chi-squared χ distributions, scaled by a factor of 9, for the diagonalentries, and uniform distributions from − A , we base our priorson the decomposition of this matrix described later on: D , and D , use Gaussian priors with a mean of 0 and standard deviationof 1.5, and θ and ϕ use uniform priors from 0 to 2 π . tar and Gas Kinematics in NGC 7000/IC 5070 l lllll l l l lll ll ll l lll ll l ll l l ll lllll lllll l llll l ll ll lll llll lll l lll llllll lll ll l lll l lll l ll l lll ll l ll l ll ll ll l lll llll llll ll l ll ll lll llll l lll ll lll lll lll l ll l l ll ll ll l lll lll ll ll lll l lllllll l l ll lll llll ll ll ll l l ll lll ll l ll lll lll l ll llll ll ll lll ll ll lll l lll l lll lll ll l ll l ll ll ll lll ll ll l ll l lll ll l l ll l ll llll ll l llll l l lll l ll llll ll lll lll ll ll ll ll lll ll ll ll ll ll ll ll ll l ll l l ll l ll l lll l l ll l lll l l ll llll ll ll ll l l l ll lll l ll lll ll l G − RP (Gaia) M G l l A V = mag l Gaia selected membersnon X−rayX−ray − J − H (2MASS) M J l A V = mag Figure 20.
Optical (left) and near-infrared (right) color-magnitude diagrams of the astrometrically validated members, com-paring stars selected using X-ray emission (black x’s) to stars selected by other criteria (magenta points). The X-ray candidatesfrom Damiani et al. (2017) are some of the only stars in our literature-based sample selected with a method that is independentof whether a stars has a disk. Nevertheless, the X-ray sample does not appear older than the other NAP members. We showthe same isochrones and reddening laws as in Figure 2. The error bars illustrate typical 2MASS photometric uncertainties.
Figure 21.
Spitzer /IRAC 3.6 µ m image of the region aroundCloud 8 (white contour). Identified YSO candidates areshown, including sources from Group D (green circles) as wellas Spitzer candidates from Rebull et al. (2011) not includedin our study (yellow squares). A 2.5 (cid:48) circle encompasses thecluster Cambr´esy 6. of generality, we can take U and V ∗ to be rotation ma-trices that rotate the coordinate systems by angles ϕ and θ , respectively, and D to be a diagonal matrix withentries D , and D , , where the diagonal elements canbe positive, negative, or zero, but | D , | ≥ | D , | . Thus, M first rotates the coordinate system by angle θ – theposition angle of the velocity anisotropy, then appliesthe velocity gradient D , , with units of km s − pc − ,along the major axis of the anisotropy and D , alongthe minor axis, and finally rotates the coordinate systemagain by ϕ . The coordinate system rotations mean thatthe component of the velocity in the outward directionwill be cos( ϕ + θ ) and the component of the velocity inthe azimuthal direction will be sin( ϕ + θ ).To interpret these results, we note that statisticallysignificant non-zero values of D , or D , imply depen-dence of velocity on position. If the velocity dependenceis primarily in the radial direction (i.e. cos[ θ + ϕ ] ≈ D , , D , > D , , D , < θ + ϕ ] ≈ ± D , and D , implyrotation.Figure 22 shows the marginal distributions for sev-eral of the interesting parameters in the linear models.For Group D, the allowed area of parameter space isvery small, so constraints on the velocity gradients aretight. Both velocity gradients are positive; with D , being noticeably larger than D , . For Groups C andE, the values of D , are constrained to be positive, butthe uncertainties in D , are large enough that this pa-rameter could equal zero. The marginal distributions ofthe position angles θ show that these are fairly tightlyconstrained for all three groups, C, D, and E. We alsoinclude marginal distributions of the normalized radial0 Kuhn et al. −0.5 0.0 0.5 1.0 − . . . . D [ km s − pc − ] D [ k m s − p c − ] Group C . . . . . position angle θ [ degrees ] p r obab ili t y den s i t y . . . . . . . normalized radial component p r obab ili t y den s i t y . . . . . . . σ scatter [ km s − ] p r obab ili t y den s i t y −0.5 0.0 0.5 1.0 − . . . . D [ km s − pc − ] D [ k m s − p c − ] Group D . . . . position angle θ [ degrees ] p r obab ili t y den s i t y normalized radial component p r obab ili t y den s i t y σ scatter [ km s − ] p r obab ili t y den s i t y −0.5 0.0 0.5 1.0 − . . . . D [ km s − pc − ] D [ k m s − p c − ] Group E . . . . . position angle θ [ degrees ] p r obab ili t y den s i t y normalized radial component p r obab ili t y den s i t y . . . . . . σ scatter [ km s − ] p r obab ili t y den s i t y Figure 22.
Marginal distributions for parameters of the linear velocity model in Groups C, D and E. These are the three stellargroups with enough sources to provide meaningful constraints on the model. For each fit, distributions are shown for the valuesof D (upper left), the position angle θ of the anisotropy (upper right), the radial component of the velocity gradient (lower left),and the random velocity scatter (lower right). The black contour lines (upper left) and vertical gray lines (other graphs) enclosethe regions of parameter space with 95% of the probability. component defined as cos( θ + ϕ ). These distributions allfavor values close to 1, showing that the velocity gradi-ents imply expansion, not rotation. E. NEARBY CLUSTERSExamination of the
Gaia µ α (cid:63) – µ δ distribution in thevicinity of the NAP region reveals two salient featuresthat are not directly related to the NAP region. One ofthese, NGC 6997, is a star cluster superimposed on the“North America” region but generally thought to be at a different distance. This group of ∼ Gaia sourcesstands out clearly in proper-motion space because it hasa smaller velocity dispersion, but it also has significantlydifferent mean proper motions of µ α (cid:63) = − . − and µ δ = − . − , meaning that confusion be-tween its members and those of the NAP complex isminimal. The group has a parallax of 1.12 mas, meaningthat it is ∼
100 pc farther than the main NAP complex. tar and Gas Kinematics in NGC 7000/IC 5070 ∼ Gaia sources cen-tered at 21:01:53 +45:12:00 with µ α (cid:63) = − . − , µ δ = − . − , and (cid:36) = 1 . Facility:
FCRAO,
Gaia , Keck:I (HIRES)
Software: astropy (Astropy Collaboration et al.2013; Price-Whelan et al. 2018), BUGS (Lunn et al.2009), class (Venables & Ripley 2002), JAGS (Plummer2017), mclust (Scrucca et al. 2016), numpy (van der Waltet al. 2011), pandas (Wes McKinney 2010; pandas devel-opment team 2020), postgres, R (R Core Team 2018),rjags (Plummer 2019), SAOImage DS9 (Joye & Man-del 2003), scipy (Virtanen et al. 2020), scikit-learn (Pe-dregosa et al. 2011), TOPCAT (Taylor 2005)REFERENCES
Adams, F. C. 2000, ApJ, 542, 964, doi: 10.1086/317052Allen, L., Megeath, S. T., Gutermuth, R., et al. 2007,Protostars and Planets V, 361Alves, J., Zucker, C., Goodman, A. A., et al. 2020, Nature,578, 237, doi: 10.1038/s41586-019-1874-zArce, H. G., Shepherd, D., Gueth, F., et al. 2007, inProtostars and Planets V, ed. B. Reipurth, D. Jewitt, &K. Keil, 245. https://arxiv.org/abs/astro-ph/0603071Argiroffi, C., Caramazza, M., Micela, G., et al. 2016, A&A,589, A113, doi: 10.1051/0004-6361/201526539Armond, T., Reipurth, B., Bally, J., & Aspin, C. 2011,A&A, 528, A125, doi: 10.1051/0004-6361/200912671Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,et al. 2013, A&A, 558, A33,doi: 10.1051/0004-6361/201322068Bally, J., Ginsburg, A., Probst, R., et al. 2014, AJ, 148,120, doi: 10.1088/0004-6256/148/6/120Bally, J., & Reipurth, B. 2003, AJ, 126, 893,doi: 10.1086/376599Bally, J., & Scoville, N. Z. 1980, ApJ, 239, 121,doi: 10.1086/158094Baudry, J.-P., Raftery, A. E., Celeux, G., Lo, K., &Gottardo, R. 2010, Journal of computational andgraphical statistics, 19, 332Baumgardt, H., & Kroupa, P. 2007, MNRAS, 380, 1589,doi: 10.1111/j.1365-2966.2007.12209.xBertoldi, F., & McKee, C. F. 1992, ApJ, 395, 140,doi: 10.1086/171638Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A,51, 207, doi: 10.1146/annurev-astro-082812-140944 Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS,427, 127, doi: 10.1111/j.1365-2966.2012.21948.xBroos, P. S., Getman, K. V., Povich, M. S., et al. 2013,ApJS, 209, 32, doi: 10.1088/0067-0049/209/2/32Cambr´esy, L., Beichman, C. A., Jarrett, T. H., & Cutri,R. M. 2002, AJ, 123, 2559, doi: 10.1086/339830C´anovas, H., Cantero, C., Cieza, L., et al. 2019, A&A, 626,A80, doi: 10.1051/0004-6361/201935321Carpenter, J. M. 2000, AJ, 120, 3139, doi: 10.1086/316845Chevance, M., Kruijssen, J. M. D., Vazquez-Semadeni, E.,et al. 2020, arXiv e-prints, arXiv:2004.06113.https://arxiv.org/abs/2004.06113Cohen, M., & Kuhi, L. V. 1979, ApJS, 41, 743,doi: 10.1086/190641Comer´on, F., & Pasquali, A. 2005, A&A, 430, 541,doi: 10.1051/0004-6361:20041788Corbally, C. J., Straiˇzys, V., & Laugalys, V. 2009, BalticAstronomy, 18, 111. https://arxiv.org/abs/0911.4287Da Rio, N., Tan, J. C., Covey, K. R., et al. 2017, ApJ, 845,105, doi: 10.3847/1538-4357/aa7a5bDamiani, F., Pillitteri, I., & Prisinzano, L. 2017, A&A, 602,A115, doi: 10.1051/0004-6361/201629402Elmegreen, B. G. 2000, ApJ, 530, 277, doi: 10.1086/308361—. 2002, ApJ, 577, 206, doi: 10.1086/342177Erickson, N. R., Grosslein, R. M., Erickson, R. B., &Weinreb, S. 1999, IEEE Transactions on MicrowaveTheory and Techniques, 47, 2212Everitt, B. S., Landau, S., Leese, M., & Stahl, D. 2011,John Wiley & SonsFeigelson, E. D. 2018, in Astrophysics and Space ScienceLibrary, Vol. 424, The Birth of Star Clusters, ed.S. Stahler, 119, doi: 10.1007/978-3-319-22801-3 5 Kuhn et al.
Feigelson, E. D., & Montmerle, T. 1999, ARA&A, 37, 363,doi: 10.1146/annurev.astro.37.1.363Findeisen, K., Hillenbrand, L., Ofek, E., et al. 2013, ApJ,768, 93, doi: 10.1088/0004-637X/768/1/93Fleming, G. D., Kirk, J. M., Ward-Thompson, D., &Pattle, K. 2019, arXiv e-prints, arXiv:1904.06980.https://arxiv.org/abs/1904.06980Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al.2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al.2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051Galli, P. A. B., Loinard, L., Bouy, H., et al. 2019, A&A,630, A137, doi: 10.1051/0004-6361/201935928Gelman, A., & Rubin, D. B. 1992, Statistical Science, 7,457, doi: 10.1214/ss/1177011136Getman, K. V., Feigelson, E. D., Garmire, G., Broos, P., &Wang, J. 2007, ApJ, 654, 316, doi: 10.1086/509112Getman, K. V., Feigelson, E. D., Kuhn, M. A., & Garmire,G. P. 2019, MNRAS, 487, 2977,doi: 10.1093/mnras/stz1457Getman, K. V., Feigelson, E. D., Sicilia-Aguilar, A., et al.2012, MNRAS, 426, 2917,doi: 10.1111/j.1365-2966.2012.21879.xGieles, M., Portegies Zwart, S. F., Baumgardt, H., et al.2006, MNRAS, 371, 793,doi: 10.1111/j.1365-2966.2006.10711.xGoldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209,doi: 10.1086/307195Gonz´alez, M., & Alfaro, E. J. 2017, MNRAS, 465, 1889,doi: 10.1093/mnras/stw2855Gouliermis, D. A. 2018, PASP, 130, 072001,doi: 10.1088/1538-3873/aac1fdGroßschedl, J. E., Alves, J., Meingast, S., et al. 2018, A&A,619, A106, doi: 10.1051/0004-6361/201833901Guieu, S., Rebull, L. M., Stauffer, J. R., et al. 2009, ApJ,697, 787, doi: 10.1088/0004-637X/697/1/787Guszejnov, D., Grudi´c, M. Y., Offner, S. S. R., et al. 2020,MNRAS, 492, 488, doi: 10.1093/mnras/stz3527Gutermuth, R. A., Pipher, J. L., Megeath, S. T., et al.2011, ApJ, 739, 84, doi: 10.1088/0004-637X/739/2/84Hacar, A., Alves, J., Burkert, A., & Goldsmith, P. 2016,A&A, 591, A104, doi: 10.1051/0004-6361/201527319Hartmann, L., Ballesteros-Paredes, J., & Heitsch, F. 2012,MNRAS, 420, 1457,doi: 10.1111/j.1365-2966.2011.20131.xHerbig, G. H. 1958, ApJ, 128, 259, doi: 10.1086/146540Herczeg, G. J., & Hillenbrand, L. A. 2015, ApJ, 808, 23,doi: 10.1088/0004-637X/808/1/23Herczeg, G. J., Kuhn, M. A., Zhou, X., et al. 2019, ApJ,878, 111, doi: 10.3847/1538-4357/ab1d67 Heyer, M., & Dame, T. M. 2015, ARA&A, 53, 583,doi: 10.1146/annurev-astro-082214-122324Heyer, M. H., & Terebey, S. 1998, ApJ, 502, 265,doi: 10.1086/305881Hillenbrand, L. A., & Baliber, N. 2015, in IAU GeneralAssembly, Vol. 29, 2253943Howarth, I. D., & van Leeuwen, F. 2019, MNRAS, 484,5350, doi: 10.1093/mnras/stz291Joye, W. A., & Mandel, E. 2003, in Astronomical Society ofthe Pacific Conference Series, Vol. 295, AstronomicalData Analysis Software and Systems XII, ed. H. E.Payne, R. I. Jedrzejewski, & R. N. Hook, 489Klessen, R. S. 2000, ApJ, 535, 869, doi: 10.1086/308854Klessen, R. S., Burkert, A., & Bate, M. R. 1998, ApJL,501, L205, doi: 10.1086/311471Kounkel, M., Covey, K., Su´arez, G., et al. 2018, AJ, 156,84, doi: 10.3847/1538-3881/aad1f1Kroupa, P., Aarseth, S., & Hurley, J. 2001, MNRAS, 321,699, doi: 10.1046/j.1365-8711.2001.04050.xKruijssen, J. M. D., Maschberger, T., Moeckel, N., et al.2012, MNRAS, 419, 841,doi: 10.1111/j.1365-2966.2011.19748.xKruijssen, J. M. D., Pelupessy, F. I., Lamers, H.J. G. L. M., Portegies Zwart, S. F., & Icke, V. 2011,MNRAS, 414, 1339,doi: 10.1111/j.1365-2966.2011.18467.xKrumholz, M. R., Matzner, C. D., & McKee, C. F. 2006,ApJ, 653, 361, doi: 10.1086/508679Krumholz, M. R., & McKee, C. F. 2020, MNRAS,doi: 10.1093/mnras/staa659Kuhn, M. A., & Feigelson, E. D. 2019, Applications inAstronomy, ed. S. Fruhwirth-Schnatter, G. Celeux, &C. P. Robert (Chapman and Hall/CRC), 463Kuhn, M. A., Hillenbrand, L. A., Sills, A., Feigelson, E. D.,& Getman, K. V. 2019, ApJ, 870, 32,doi: 10.3847/1538-4357/aaef8cKutner, M. L., & Ulich, B. L. 1981, ApJ, 250, 341,doi: 10.1086/159380Lada, C. J., Lombardi, M., Roman-Zuniga, C., Forbrich, J.,& Alves, J. F. 2013, ApJ, 778, 133,doi: 10.1088/0004-637X/778/2/133Larson, R. B. 1981, MNRAS, 194, 809,doi: 10.1093/mnras/194.4.809Laugalys, V., Straiˇzys, V., Vrba, F. J., et al. 2006, BalticAstronomy, 15, 483Leung, H. W., & Bovy, J. 2019, MNRAS, 489, 2079,doi: 10.1093/mnras/stz2245Li, Z.-Y., & Nakamura, F. 2006, ApJL, 640, L187,doi: 10.1086/503419 tar and Gas Kinematics in NGC 7000/IC 5070 Kuhn et al. tar and Gas Kinematics in NGC 7000/IC 5070 Table 1.
Summary of Candidate NAP Members from the LiteratureReference Method Number of Ref.
Gaia
MembershipCandidates Notes Matches N.-Mem. Mem.(1) (2) (3) (4) (5) (6) (7)Merrill & Burwell (1949) H α a
60% 3 0Herbig (1958) H α
68 40% 6 21Welin (1973) H α b
89% 113 13Cohen & Kuhi (1979) H α c
57% 1 11Bally & Scoville (1980) IR 11 36% 3 1Marcy (1980) H α α
32 47% 1 14Comer´on & Pasquali (2005) Spec. 8 d
88% 5 2Laugalys et al. (2006) Phot./H α e
84% 340 21Witham et al. (2008) H α a
82% 7 25Straiˇzys & Laugalys (2008) CMD 5 60% 3 0Corbally et al. (2009) Spec. 34 f
62% 18 3Guieu et al. (2009) IRE 1,657 24% 125 272Rebull et al. (2011) IRE 1,329 g
29% 140 252Armond et al. (2011) H α /Clust. 54 22% 1 11Damiani et al. (2017) X-ray 721 19% 46 93H α
123 66% 22 59IRE 179 51% 58 34Combined sample All 3,473 35% 814 395
Note —References for published lists of candidate members of the NAP region. Column 2 is method for selection, includingH α emission, spectral classification, placement on the color-magnitude diagram (CMD), brightness in the infrared (IR) or IRexcess (IRE), spatial clustering, and X-ray emission. Column 3 is the number of candidates in the paper. Column 4 givesnotes about the reference. Column 5 is the percentage of sources with Gaia
DR2 counterparts that pass our quality criteria.Column 6 is the number of objects that we reject as members, and Column 7 is the number of objects whose membership wehave validated. The final row gives the statistics of the combined sample; the numbers of sources (Columns 3, 6, and 7) inpreceding rows do not sum to the values in this row because many objects have been repeatedly selected by multiple studies.aObjects were cataloged as H α emission stars in the referenced papers but later upgraded to potential NAP members by Rebullet al. (2011).bWe use the corrected list of 142 objects available from ftp://ftp.lowell.edu/pub/bas/starcats/welin.cyg, rather than the origi-nally published list of 141 objects.c The H α sample was enlarged by including stars spatially associated with known H α objects.dIn addition to the Bajamar Star, we consider all sources from Comer´on & Pasquali (2005) not classified as AGB stars to becandidate members.e We include all 430 stars in the direction of L935. However, a subset of 41 stars are flagged as having H α emission, 20 of whichare classified as non-members and 4 as members by our analysis.f We include all objects for which spectra are provided, including 19 stars with emission lines and 15 without. We classify 3emission-line stars as members and 9 as non-members, as well as 9 non-emission-line stars as non-members.gThis combines 1,286 YSO candidates identified in a Spitzer /MIPS-based search with 43 new IRAC-only candidates. Kuhn et al.
Table 2.
Astrometrically Validated Members
Gaia
DR2 α δ
Group(ICRS) (ICRS)(1) (2) (3) (3)2163138601938577664 20 51 19.8 +44 23 06 D2163137742945112960 20 51 20.6 +44 20 32 C2163138705017719296 20 51 21.3 +44 24 05 D2163136123738466688 20 51 12.0 +44 18 47 C2163137742945115136 20 51 22.6 +44 21 07 C2163148772421081728 20 51 22.8 +44 33 42 D2163149665774282368 20 51 23.6 +44 35 22 D2066862546309607552 20 51 26.4 +43 53 12 D2162947424350344448 20 51 24.4 +44 13 04 C2163135883219217280 20 51 24.6 +44 17 54 D
Note — Gaia
DR2 counterparts for the 395 previously pub-lished YSO candidates that passed our astrometric mem-bership criteria. Column 1: Source designation in the DR2catalog. Columns 2–3: Source positions are truncated,not rounded. We encourage readers to obtain the full-precision
Gaia astrometry for these sources directly fromthe
Gaia
Archive (https://gea.esac.esa.int/archive/). Col-umn 4: The kinematic group assignment for each star.(This table is available in its entirety in a machine-readableform in the online journal. A portion is shown here forguidance regarding its form and content.) tar and Gas Kinematics in NGC 7000/IC 5070 Table 3.
Astrometric Properties of Stellar GroupsGroup α δ N samp µ α (cid:63) , µ δ, (cid:36) r hm (ICRS) (ICRS) (stars) (mas yr − ) (mas yr − ) (mas) (pc)(1) (2) (3) (4) (5) (6) (7) (8)A 20 47 50 +43 47 27 21 − . ± . − . ± .
14 1 . ± .
02 1 .
1B 20 50 20 +44 37 58 19 − . ± . − . ± .
18 1 . ± .
03 1 .
8C 20 51 10 +44 19 25 47 − . ± . − . ± .
17 1 . ± .
01 1 .
1D 20 52 50 +44 22 13 235 − . ± . − . ± .
09 1 . ± .
01 1 .
9E 20 57 30 +43 46 36 59 − . ± . − . ± .
15 1 . ± .
02 1 .
3F 20 58 50 +44 09 43 14 − . ± . − . ± .
25 1 . ± .
02 2 . Note —Column 1: Name of the stellar group. Columns 2–3: Approximate coordinates of the group center. Column 4: Numberof
Gaia -validated members assigned to each group. Columns 5–6: Mean proper motions of the group. Column 7: Meanparallax of the group. Column 8: Characteristic radius for the stellar group, defined as the median distance of group members,in projection, from the group center. All values are in the
Gaia
DR2 system, with no correction for zero-point offset.
Table 4.
Clouds Properties from the CO Line MapDesig. α δ M CO v lsr distribution r cloud Optical Groupmean std. dev. mode FWHM(ICRS) (ICRS) ( M (cid:12) ) (km s − ) (km s − ) (km s − ) (km s − ) (pc)(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11)1 20 48 58 +43 48 10 330 0 . . . . . − . . − . . . † /Dark B,C3 20 51 12 +44 48 00 1100 − . . − . . . . . . . . . . . . . . . . . . − . . − . . . − . . − . . . . . . . . . . . . . . . . . . . . . . . . . . . . Note — Column 1: Cloud designation. Columns 2–3: Coordinates of the cloud center. Column 4: CO mass estimate of thecloud. Columns 5: Mean velocity of the cloud. Column 6: Standard deviation of the velocity distribution. Column 7: Modeof the velocity distribution. Column 8: ∆ v FWHM for the most-prominent mode. Column 9: Characteristic size of the cloud,defined as the radius that contains half the integrated emission. Column 10: Appearance of the cloud (or nebular material inthe direction of the cloud) in the optical image. Column 12: Group of stars associated with the cloud. † BRC = bright rim cloud. Kuhn et al.
Table 5.
Cloud Dynamical PropertiesCloud log ¯ ρ τ cross τ ff M M dyn α BM92 (g cm − ) (Myr) (Myr) (10 M (cid:12) )(1) (2) (3) (4) (5) (6) (7)1 − .
40 1 . . . . . − .
17 1 . . . . . − .
10 1 . . . . . − .
01 1 . . . . . − .
80 1 . . . . . − .
98 0 . . . . . − .
24 0 . . . . . − .
65 1 . . . . . − .
52 1 . . . . . − .
67 0 . . . . . − .
01 0 . . . . . − .
38 0 . . . . . − .
26 1 . . . . . Note — Column 1: Cloud designation. Column 2: Mean densityof the cloud, defined as ¯ ρ = M CO / (4 / πr . Column 3:Crossing timescale, defined as τ cross = r cloud / ∆ v . Column 4:Free-fall timescale defined using ¯ ρ from Column 3. Column 5:Mach number. Column 6: Dynamical mass of the cloud assuming α = 1 (Equation 9). Column 7: The Bertoldi & McKee (1992)estimate of the virial parameter, using M CO for the mass, r cloud for the radius, and σ nt for the velocity dispersion. tar and Gas Kinematics in NGC 7000/IC 5070 Table 6.
Internal Kinematic Properties of Stellar GroupsGroup Kendall’s τ Mean Motions Linear Expansion Model Parameters Total p x p y mean v out mean v az Gradient 1 Gradient 2 θ σ scatter σ (km s − ) (km s − ) (km s − pc − ) (km s − pc − ) ( ◦ ) (km s − ) (km s − )(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)A > > . ± . . ± . · · · · · · · · · · · · . ± .
2B 0.008 0.002 0 . ± . − . ± . · · · · · · · · · · · · . ± . < > . ± . . ± . . ± .
11 0 . ± .
17 81 ±
12 0 . ± . . ± . < < . ± . . ± . . ± .
03 0 . ± .
03 98 ± . ± . . ± . < > . ± . − . ± . . ± .
08 0 . ± .
07 87 ± . ± . . ± . > > − . ± . − . ± . · · · · · · · · · · · · . ± . Note —Column 1: Stellar group. Columns 2–3: Non-parametric tests for correlation between position and velocity in rightascension and declination. Columns 4–5: Mean outward velocity and mean azimuthal velocity as calculated in Section 7.2.Columns 6–9: Parameters from the linear model fit to the velocity from Section 7.3, including the velocity gradients along thedirection of maximum gradient (Column 6) and in the orthogonal direction (Column 7), the position angle of the anisotropy(Column 8), and the intrinsic velocity scatter left over after the velocity gradient is removed (Column 9). Column 10: Thetotal velocity dispersion. Kuhn et al.
Table 7.
New Member Candidates from
GaiaGaia
DR2 α δ
Group(ICRS) (ICRS)2067063447700277504 20 50 06.05 +44 17 48.9 C2166282204462982528 20 50 07.33 +45 49 22.4 distrib2067060007428506112 20 50 07.69 +44 08 30.1 D2163214738816798336 20 50 10.93 +44 51 22.0 D2065821858551580288 20 50 12.03 +41 38 44.1 distrib2163266832483920128 20 50 12.94 +45 32 43.1 distrib2066624635182982144 20 50 14.35 +42 15 43.0 distrib2067049428926311040 20 50 15.04 +43 58 11.8 D2166282964674128256 20 50 15.47 +45 54 02.3 distrib2067049532005529600 20 50 15.81 +43 58 59.9 A2067060625914494080 20 50 16.82 +44 11 53.5 B