A Gaia view on the star formation in the Monoceros OB1 and R1 associations

Stellar kinematics provides the key to understanding star formation process. In this respect, we present a kinematic study of the Monoceros OB1 (Mon OB1) and R1 (Mon R1) associations using the recent Gaia data and radial velocities of stars derived from high-resolution spectroscopy and the literature. A total of 728 members are selected using the criteria based on the intrinsic properties of young stars, parallaxes, and proper motions. The spatial distribution and kinematic properties of members show that these associations have distinct substructures. In Mon OB1, we find one northern group and two southern groups. Mon R1 is composed of three small stellar groups that are spatially and kinematically distinct. Some stars are found in a halo around these two associations. We detect patterns of expansion for most stellar groups in the associations. In addition, two stellar groups in Mon OB1 show the signature of rotation, which provides an important constraint on cluster formation. The star formation history of Mon OB1 is slightly revised. Star formation first occurred in the southern region and subsequently in the northern region. Recent star-forming events ignited deeper into the southern region, while some stars are escaping from Mon OB1, forming a halo. Mon R1 might have formed at the same epoch as the formation of the northern group in Mon OB1. Given that star formation is taking place on different scales along a large arc-like structure, Mon OB1 and Mon R1 may be the results of hierarchical star formation.


INTRODUCTION
Most young stars form in stellar systems such as stellar clusters and associations (Lada & Lada 2003;Porras et al. 2003). It is expected that only less than 10% of young stellar groups will remain gravitationally bound clusters (Lada & Lada 2003), so a significant portion of field stars may originate from the dissolution of such stellar * FNRS Senior Research Associate systems (Miller & Scalo 1978;Briceño et al. 2007). The formation of stellar systems is interconnected in space and time, which leads to form larger structures (Elmegreen et al. 2000;Gouliermis 2018). Therefore, they are superb laboratories to understand the star formation taking place on various spatial scales.
Stellar associations are, in general, composed of a single or multiple stellar clusters (or groups) and a distributed young stellar population (Blaauw 1964;Koenig et al. 2008). In addition, the internal structure is tightly associated with the kinematics of constituent stars (Lim et al. 2018(Lim et al. , 2019a(Lim et al. , 2020(Lim et al. , 2021. The morphological features and stellar kinematics provide hints in understanding the formation process of stellar associations. To obtain such inference, the high precision astrometry from the Gaia mission (Gaia Collaboration et al. 2016) is key as it allows us to select genuine members spread over a wide region and investigate their kinematics in detail, especially when combined with radial velocity (RV) data.
Extensive RV surveys have been performed for young stars in NGC 2264, the core of Mon OB1 (Fűrész et al. 2006;Tobin et al. 2015). The velocities of young stars follow the velocity field of the remaining molecular gas. A group of stars with RVs larger than the systemic velocity (∼ 5 km s −1 ) was found toward the O-type binary S Monocerotis (Skiff 2009). Tobin et al. (2015) claimed that this group might have formed on the far side of the remaining gas compressed by the strong wind from the massive star. They also reported the presence of another group of stars with systematically small RVs. Recently, the internal kinematics of this SFR has been investigated with the Gaia proper motion (PM) data (Kuhn et al. 2019;Buckner et al. 2020). As a result, a pattern of expansion was detected in this stellar group at the north of NGC2264.
Mon R1 is about 2 • west away from Mon OB1. It is composed of small SFRs surrounded by reflection nebulae such as IC 446, IC 447, NGC 2245, and NGC 2247. Infrared observations revealed several young stellar groups around Herbig Ae/Be stars in Mon R1 (Wang & Looney 2007;Gutermuth et al. 2009). Some Herbig-Haro objects were also discovered in those groups (Movsessian et al. 2021). The distance to Mon R1 was previously determined in the range of 660 pc to 715 pc (van den Bergh 1966; Monteiro et al. 2020;Dias et al. 2021; Movsessian et al. 2021). This association has thus been considered to be at the same distance as Mon OB1.
The molecular CO line observations of Mon R1 were carried out by Kutner et al. (1979). They found a semi ring-like structure, which is also evident in Planck continuum image at 550 µm (Bhadari et al. 2020). The overall velocity field is distinguishable from that of Mon OB1 (see also Oliver et al. 1996). In the direction of Mon R1, Kutner et al. (1979) identified two velocity components at −1 to 1 km s −1 and 3 to 5 km s −1 . The former component is physically associated with IC 446 and IC 447, while the latter one belongs to NGC 2245 and NGC 2247. Recently, Bhadari et al. (2020) considered the former component as a filament and discussed the star formation process along the filament in the context of the end-dominated collapse model (Bastien 1983;Pon et al. 2011). However, the kinematic properties of young stars in Mon R1 have not yet been studied in detail.
In this study, we aim to investigate not only the star formation process in Mon OB1 and Mon R1 but also the physical association between them by probing the kinematics of young stars. For this purpose, the recent Gaia Early Data Release 3 (EDR3; Gaia Collaboration et al. 2021) is used with RV data. We describe data and target selection in Section 2. The scheme of member selection is addressed in Section 3. In Section 4, we investigate the substructures in the two associations and the kinematic properties of young stars in the substructures. Star formation history is also inferred from a color-magnitude diagram (CMD). Star and cluster formation is discussed in Section 5. Our results are summarized in Section 6.
2. DATA 2.1. Selection of member candidates A 6 • × 6 • region centered at R.A. = 06 h 36 m 23. s 52, decl. = +10 • 04 55. 7 (J2000) was surveyed. In order to minimize the inclusion of field interlopers in the field of view, we first isolated member candidates based on the intrinsic properties of young stars as done in Lim et al. (2021).
Early-type (O-and B-type) stars are probable member candidates because of their short lifetime, particularly that of O-type stars. We compiled lists of earlytype stars taken from the data bases of MK classifications (Wenger et al. 2000;Reed 2003;Mermilliod & Paunzen 2003;Skiff 2009;Maíz Apellániz et al. 2013) and removed some duplicated stars. A total of 609 earlytype stars were selected as member candidates.
Low-mass young stellar objects (YSOs) with warm circumstellar disks appear bright at infrared wavelengths (Lada 1987). Hydrogen recombination lines are ob- Figure 1. Color-color diagrams of YSO candidates classified with the AllWISE data (Cutri et al. 2014). Red triangle, green square, and orange diamond represent Class I, Class II, YSO with a transitional disk, respectively. The color criteria for the YSO classification of Koenig & Leisawitz (2014) are shown by dashed lines. served in emission due to mass accretion (Muzerolle et al. 1998(Muzerolle et al. , 2003Fang et al. 2013). X-ray is emitted from their hot coronal region (Flaccomio et al. 2003;Caramazza et al. 2012;Rauw & Nazé 2016, etc). Based on these properties, the membership of young stars in NGC 2264 has been thoroughly evaluated by a series of multiwavelength studies (Sung et al. 1997;Park et al. 2000;Sung et al. 2004Sung et al. , 2008Sung et al. , 2009). We built a list of 992 member candidates from these studies.
The Wide-field Infrared Survey Explorer (WISE) has mapped whole sky in mid-infrared passbands (Wright et al. 2010). We used AllWISE catalogue (Cutri et al. 2014) to identify YSO candidates spread over the entire survey region. First, a number of spurious sources were rejected by adopting the criterion nm/m ≤ 0.2 in each passband (Koenig & Leisawitz 2014), where nm and m are the number of profile-fit flux measurements for a given source with signal-to-noise ratios (SNRs) larger than 3 and the total number of profile-fit flux measurements for the same source in a given passband, respectively. There could still be many spurious sources in W3 and W4 bands. We suppressed those sources following Koenig & Leisawitz (2014); where χ 2 is the reduced chi-square of profile-fit in a given passband.
Additional contaminants such as active galactic nuclei and star-forming galaxies were excluded by using the criteria of Koenig & Leisawitz (2014). In the end, we classified 75 Class I, 348 Class II, and three YSOs with transitional disks according to the scheme of Koenig & Leisawitz (2014) as shown in Figure 1. These YSO candidates were considered as member candidates.
We cross-matched the three lists that we built to create a master catalogue of member candidates. A total of 29 out of 609 early-type stars from the data bases of MK classification were found in the member candidate list of NGC 2264 (Sung et al. 1997;Park et al. 2000;Sung et al. 2004Sung et al. , 2008Sung et al. , 2009). There are 124 candidates in common between the YSO list from the AllWISE data and that of member candidates in NGC 2264. Therefore, a total of 302 sources are additional YSO candidates throughout our survey region, of which two are early-type stars, i.e. they are in the catalogue of 609 early-type stars. The total number of member candidates is 1872.
We searched for the counterparts of all member candidates in the catalogue of the Gaia EDR3 (Gaia Collaboration et al. 2021) within a radius of 3 . All the OB star candidates were found in the Gaia data. The Gaia data further contains 969 out of 992 member candidates belonging to NGC 2264 and 336 out of 426 YSO candidates identified with AllWISE data. A total of 1622 candidates brighter than 18 mag in G RP were used for member selection and analysis. Tobin et al. (2015) published the RV data of 695 stars located in NGC 2264. It is noted that the published RVs are actually line of sight velocities at the local standard of rest, not heliocentric RVs. We took the data and cross-matched them with the Gaia catalogue, leading to a total of 684 having Gaia counterparts.

Radial velocity measurements
We obtained the high-resolution spectra of 14 YSO candidates in Mon R1 using the Immersion GRating Infrared Spectrometer (IGRINS, R ∼ 45, 000 - Yuk et al. 2010;Park et al. 2014) attached to the 8.2-m Gemini South telescope on 2020 February 4, 5, 7, 9, 11, and 12. An ABBA nod technique was applied to our observations to subtract the sky background. Some A0V stars, such as HIP 30387, HIP 36796, HIP 33297, and HIP 28686, were observed as telluric standard stars. A large number of OH emission lines were observed from blank sky regions for wavelength calibration.
Data reduction was performed by using the IGRINS pipeline package version 2 (Lee et al. 2017) 1 . This pipeline sequentially executes aperture extraction, the subtraction of background emission, bad pixel correction, and wavelength calibration. The synthetic spectrum of Vega (Castelli & Kurucz 2004) was fit to those of the observed A0V stars. The telluric spectra were We plot only stars with parallaxes larger than three times their associated errors. Dashed lines in the left panel confine the probable members to distances between 500 pc and 1000 pc. The right panel shows the PM distribution of member candidates between 500 pc and 1000 pc. Red, blue, and black dots represent the genuine members of Mon OB1, Mon R1, and the halo, respectively, while gray dots denote probable nonmembers (see the main text for detail).
obtained from the spectra of the A0V stars divided by the best-fit synthetic spectrum of Vega. Target spectra were corrected by the telluric spectra. Synthetic stellar spectra in the wide temperature range of 3500 to 9000 K for the solar abundance were generated using SPECTRUM v2.76 (Gray & Corbally 1994) 2 based on a grid of the ODFNEW model atmospheres (Castelli & Kurucz 2004). The wavelength of the synthetic spectra in air was converted to that in vacuum by using the relation of Ciddor (1996). We derived the cross-correlation functions between the synthetic spectra and the observed spectra of the 14 YSO candidates with xcsao task in the RVSAO package . The velocities at the strongest correlation peaks were adopted as RVs. The task xcsao yields the RV uncertainty based on the r value as below (Tonry & Davis 1979): where h and σ a represent the amplitude of a crosscorrelation function and the root mean square value of its antisymmetric component, respectively. The RV uncertainty is then obtained from the relation 3w/8(1 + r) where w is the full width at half maximum of the peak of cross-correlation function . Some spectral orders showed poor cross-correlation functions because of the small number of lines. The r values tend to be lower than 6.0 for these orders. Therefore, we adopted the weighted-mean value and standard deviation of RVs measured from spectral orders with r larger than 6.0 as the final RV and RV error of a given YSO candidate, respectively. The inverse of squared uncertainty was used as the weight value. The RVs of YSO candidates were then converted to velocities in the local standard of rest frame using the IRAF/rvcorrect task.

Supplementary data
The Infrared Astronomical Satellite (IRAS) mission surveyed more than 95% of sky at 12, 25, 60, and 100 µm (Neugebauer et al. 1984). Later, the Improved Reprocessing of the IRAS Survey (IRIS) provided better quality of dust images over the sky (Miville-Deschênes & Lagache 2005). In addition, the AKARI satellite mapped almost all sky in the four far-infrared bands centered at 65, 90, 140, and 160µm (Doi et al. 2015).
These infrared maps help to investigate the distribution of interstellar material around Mon OB1 and Mon R1. We took the IRIS image at 100µm and the AKARI Far-Infrared Surveyor false-color image of our survey region (Miville-Deschênes & Lagache 2005;Doi et al. 2015) processed by the Aladin interactive sky atlas (Bonnarel et al. 2000;Boch & Fernique 2014).

MEMBER SELECTION
We may assume that young stars in an SFR have formed in the same molecular cloud. Therefore, they are almost at the same distance and share similar kinematic properties. Based on this conventional idea, we assessed the membership of the member candidates using the Gaia parallax and PM data (Gaia Collaboration et al. 2021). Systematic zero-point offsets that depend on magnitude, color, and position were found in the Gaia parallax (Lindegren et al. 2021). We corrected parallaxes for the zero-point offsets according to the recipe of Lindegren et al. (2021, https://gitlab.com/ icc-ub/public/gaiadr3 zeropoint).
Stars with parallaxes smaller than three times the associated errors were excluded in member selection. In addition, we did not use stars with negative parallaxes or close companions (duplication flag = 1), or poor astrometric solutions (RUWE > 1.4) were not used in analysis as well as member selection.
The left panel in Figure 2 displays the parallax distributions of member candidates. We considered only member candidates between 500 pc and 1000 pc given the previously determined distances of the two associations (∼ 700 pc van den Bergh 1966; Sung et al. 1997;Monteiro et al. 2020;Movsessian et al. 2021). We plot the PMs of members fulfilling this criterion in the right panel of Figure 2.
There are a few groups that have different PMs, on average. These groups may be related to Mon OB1 and Mon R1. We assigned the member candidates to each association based on the spatial distribution in Figure 3. There are the well-populated association Mon OB1 to the east and the loose association Mon R1 to the west. We considered the boundary of Mon OB1 as a red ellipse with a semi-major axis of 40 and an eccentricity of 0.7 centered at R.A. = 06 h 40 m 54. s 56, decl. = +09 • 39 43. 7 (J2000). The boundary of Mon R1 was assumed to be a blue circle with a radius of 35 centered at R.A. = 06 h 31 m 39. s 11, decl. = +10 • 05 55. 7 (J2000). A total of 653 and 48 stars were found in Mon OB1 and Mon R1, respectively, from these criteria.
An iterative process of removing PM outliers was performed to select a reliable set of members. We computed the weighted mean values and standard deviations of the PMs from the stars in each association. The inverse of squared error was adopted as the weight value. Then, stars with PMs within the standard deviations (1σ) from the mean PMs were used to determine better mean PMs and their standard deviations. These new statistical values were used as the initial values to select members. In each region, stars whose PMs are within four times the standard deviations (4σ) from the newly determined mean PMs were selected as members.  (Lindegren et al. 2021). We used stars with parallaxes larger than five times their associated errors. The black curves shows the best-fit Gaussian distributions.
The latter criteria was used to avoid eliminating possible member candidates. We redetermined the weighted mean PMs and standard deviations using the members. This procedure was repeated until the statistical values converge to constant values. The numbers of the final members in Mon OB1 and Mon R1 are 631 and 46 in total, and they are shown by red and blue dots in Figures 2 and 3, respectively. A total of 53 YSO candidates were also found between the two associations. All of them are sparsely distributed within an elliptical region (dashed curve in Figure 3). We refer to this low-stellar density region as halo and assume that there is no member outside the halo. The YSO members were selected in the same manner as above. As a result, we selected 51 YSO members in the halo. The halo YSO members are shown by black dots in Figures 2 and 3.
There are a number of early-type stars that do not belong to the two associations. These stars, on average, have larger PMs and PM dispersion than those of genuine members within the associations (most gray dots in the right panel of Figure 2). In addition, they are uniformly distributed over the surveyed region, except for the western part obscured by dark clouds. These facts implies that most of them may not be genuine mem-bers. Therefore, we did not additionally select earlytype members in the halo.
A total of 728 stars (631 in Mon OB1, 46 in Mon R1, and 51 in the halo) were finally selected as genuine association members. The list of members is presented in Table 1. We plot the distance distributions of these stars in Figure 4. The distance distribution of each group was fit by a Gaussian distribution. The distances to Mon OB1, Mon R1, and the halo were determined to be 704 ± 38 (s.d.) pc, 660 ± 35 (s.d.) pc, and 700 ± 47 (s.d.) pc, respectively, from the best-fit Gaussian distributions. These results are consistent with those of previous studies (Becker & Fenkart 1971;Feldbrugge & van Genderen 1991;Sung et al. 1997;van den Bergh 1966;Monteiro et al. 2020;Dias et al. 2021;Movsessian et al. 2021) within errors. The members of Mon R1 are systematically closer than Mon OB1 although there is only an 1σ level difference.

Substructures
A number of previous studies have probed substructures in many SFRs from the spatial distributions of stars. However, in the absence of kinematic information, it is unclear whether substructures are real physical systems. In this study, we search for substructures using the Gaia PM and RV data, as well as the spatial distribution of members.

Mon OB1
A deep photometric study of Sung et al. (2008) in optical passbands showed that NGC 2264 is composed of two active SFRs and a halo. The northern group is located around the massive O-type binary S Monocerotis, while the southern group is in the vicinity of the Cone Nebula. A halo surrounds these two SFRs. Later, Spitzer observations revealed that the southern group is composed of two subgroups of YSOs (Spokes and Cone(C), Teixeira et al. 2006;Sung et al. 2009). Additional smaller-scale substructures were identified in that southern group (Kuhn et al. 2014).
We searched for substructures in Mon OB1 from the correlations between PMs and positions of members. As a result, three stellar groups that are spatially and kinematically distinct were identified from the correlation between µ α cos δ (PM along R.A.) and declination in the upper left panel of Figure 5. We checked the spatial distributions of members within the three ellipses in the figure and found that members were populated in three specific regions, i.e. they do not show a random distribution across this association. We divided members into   (2) and (3) : The equatorial coordinates of members. Columns (4) and (5) : Absolute parallax and its standard error. Columns (6) and (7) : PM in the direction of right ascension and its standard error. Columns (8) and (9) (14) and (15) : Radial velocity at the local standard of rest and its error. Column (16) : Member type. 'E' represents O-or B-type stars obtained from the data bases of MK classification (Wenger et al. 2000;Reed 2003;Mermilliod & Paunzen 2003;Skiff 2009;Maíz Apellániz et al. 2013).
'Y' denotes young stellar objects or pre-main sequence members. Column (17)  A. and positions along declination. It is evident that there are three groups of stars that are spatially and kinematically distinct. Three ellipses were used to ascertain the spatial distributions of members within them. The solid lines represent the arbitrary boundaries of each group. The other panels display the spatial distributions of members belonging to each group. Red, blue, and orange dots represent the members of the S Mon group, the Cone group, and the THF15 group, respectively. The size of dots is proportional to the brightness of individual stars.
three groups by three solid lines as shown in Figure 5. These lines were adjusted through the visual inspection of their spatial distributions. Note that this division could be somewhat arbitrary because it is difficult to define the exact boundaries of each group. The northern group (hereafter S Mon group, red symbol in Figure 5) has a median PM (µ α cos δ, µ δ ) of (−1.649 mas yr −1 , −3.655 mas yr −1 ). We could not find additional clustering in the group. On the other hand, the southern group is separated into two groups. The group shown by blue symbol (the lower-left panel) has a median PM of (−2.455 mas yr −1 , −3.721 mas yr −1 ): we refer to this southern group as the Cone group according to the nomenclature of Sung et al. (2008). The other group (orange symbol in the lower-right panel) has a median PM of (−1.676 mas yr −1 , −3.685 mas yr −1 ). This group seems to correspond to the blueshifted population reported by Tobin et al. (2015), and therefore we refer to this group as THF15 from the names of the authors. There is no smaller-scale substructure within these two groups as well given the absence of additional enhancement in stellar surface density. The central coordinates of these three groups were obtained from the median positions of the relevant members. We summarize the basic properties of these groups in Table 2.   Table 2). The kinematic substructures can also be confirmed from the correlations between RVs and positions of members. Especially, there is a gradient of RVs with respect to declination (∼ 0.4 km s −1 pc −1 ). The median RVs of the S Mon group, Cone, and THF15 are about 5.7, 4.3, and 2.3 km s −1 , respectively. The members in THF15 have RVs systematically smaller than the others as reported by Tobin et al. (2015). Figure 7 displays the spatial distribution of stars in Mon R1. This stellar association is composed of three small stellar groups associated with reflection nebulae. The young open cluster Collinder 95 occupies the southern part of Mon R1. The reflection nebula IC 447 may have been formed by the early-type members of this cluster. We referred to this group as IC 447. A partially embedded group of stars in the vicinity of the reflection nebula IC 446 are found to the northwest of this association (hereafter IC 446). The other embedded group lies between the reflection nebulae NGC 2245 and NGC 2247 in the eastern region (hereafter N2245/47).  The boundaries of these groups were determined from the spatial distribution of members. Members below the declination of ∆decl. = −0. 5 were assigned to the members of IC 447. The members of IC 446 were confined to the northwestern stars (∆R.A. < −65. 0 and ∆decl. > 10. 0). The rest of stars were considered as the members of NGC 2245/47. We plot V R.A. , V decl. , and RVs of the Mon R1 members in Figure 8. The three stellar groups show somewhat complicated kinematic substructures. The members of N2245/47 and IC 446 show a large scatter in tangential velocities compared to the members of IC 447. Also, the V decl. of stars continuously vary with declination from IC 447 to IC 446 (∼ 0.5 km s −1 pc −1 ). The YSOs with RV measurements are found in only the two groups IC 447 and N2245/47, and their total number is somewhat limited to probe the global variation. Nevertheless, the RVs of the members show a tendency similar to the tangen- tial velocities. The RVs of IC 447 and N2245/47 follow the velocity fields of the remaining molecular gas. The former and latter correspond to the gas components in the RV ranges of −4 to 2 km s −1 and 2 to 10 km s −1 (Bhadari et al. 2020), respectively. It suggests that the complicated kinematics of stars seen in the tangential velocity distributions might have been inherited from that of their natal cloud. We determined the central coordinates of the three groups from the median positions of the associated members. Their basic properties are summarized in Table 2.

Internal kinematics in Mon OB1
We investigated the internal motions of stellar groups in Mon OB1. The upper left panel of Figure 9 yields the spatial distribution of stars along with the PM vectors (arrows) of three stellar groups relative to the systemic motion of Mon OB1, where the relative PM vec- tors are the median PM vectors of individual groups subtracted by the median PM of the entire system. The three groups do not have any significant motion along declination. The S Mon and THF15 groups are moving eastward with similar velocities, while the Cone group is moving westward at a larger velocity.
The relative PM vectors of individual members within a given group were obtained after subtracting their median PM. Figure  In order to quantitatively probe the internal motions of stars, we measured the vectorial angles of members as used in our previous studies (Lim et al. 2019a(Lim et al. , 2020(Lim et al. , 2021. The vectorial angle (Φ) is an angle between the position vector from the group center and the relative PM of a given star. A zero value means that a star is radially escaping from its host group. We present the Φ distribution of individual members belonging to each group in Figure 10. The Φ values of members in the S Mon group are clustered around 0 • , which is indicative of expansion. The Cone group also shows a pattern of expansion as many members of this group have Φ values around 0 • . On the other hand, the members of THF15 do not show clear outward motions given that there is no strong peak around Φ = 0 • .
Recent studies also detected expansion for the S Mon group (Kuhn et al. 2019;Buckner et al. 2020), but not for the southern groups. The southern stellar groups are seen as a single cluster in optical passbands (Sung et al. 2008), while two well-defined subgroups, the so-called Cone(C) and Spokes (Teixeira et al. 2006;Sung et al. 2009), are found around the embedded YSOs NGC 2264 IRS 1 and IRS 2 in infrared passbands. The fact that a large fraction of YSOs in the Cone(C) and Spokes groups were not found in optical CMD indicates that these two groups are deeply embedded (Sung & Bessell 2010). The presence of molecular clouds toward the southern region can directly be confirmed from Tobin et al. (2015). Therefore, we speculate that there are four stellar groups (Cone, THF15, Cone(C), and Spokes) along the line of sight at the south of Mon OB1.
If this is true, Kuhn et al. (2019) might have probed small parts of the Cone and THF15 groups, not actually the Cone(C) and Spokes groups. Then, the investigators could not find any pattern of expansion.

Internal kinematics in Mon R1
We investigated the kinematics of stellar groups in Mon R1 in the same way. The left panel of Figure 11 displays the motions of the three groups relative to the systemic motion of Mon R1. IC 447 is moving south, while N2245/47 and IC 446 are moving northeast and northwest, respectively. These groups are receding away from each other. The PM vectors of individual members relative to their host groups are shown in the right panel. Although it is difficult to define the center of each group due to the small number of stars, members tend to show outward motions. Figure 12 shows the Φ distributions. Stars in the three groups have Φ values around 0 • , indicating that the members of these groups are scattered radially outward.

Rotation
We searched for the signature of rotation for the three groups in Mon OB1 using the RVs of members as done by previous studies (Lane et al. 2009;Mackey et al. 2013;Lim et al. 2019b;Lim et al. 2021). A projected rotational axis passing through the central position of a given group was set at a position angle of 0 • (from north to south) in the projected sky plane. We computed the difference between the mean RVs of stars in the two areas separated by the axis. with an interval of 20 • in a counterclockwise direction (north to east). If a given cluster is rotating, the mean RV differences appear as a sinusoidal curve.
We found the signature of rotation for the members within a projected radius of 7 for the S Mon group and within 5 for the Cone group. Figure 13 exhibits the variations of the mean RV differences with respect to position angles. These observed variations were fit by the sinusoidal curve : where V rot , i, and Θ 0 represent the rotational velocity, inclination angle of a rotational axis along the line of sight, and phase, respectively. The amplitude of the best-fit sinusoidal curve corresponds to twice the projected rotational velocity (V rot sin i). The S Mon and Cone groups are rotating at 0.23 ± 0.02 km s −1 and 0.87 ± 0.03 km s −1 , respectively, if we assume an inclination angle of 90 • . The position angle of the projected rotational axis can be estimated from 270 • − Θ 0 . The projected rotational axes of the S Mon and Cone groups almost lie from north to south, but these groups are rotating in opposite directions. On the other hand, we could not find any signature of rotation for the THF15 group. The same method could not be applied to the stars in Mon R1 because of a small number of stars.  Figure 14 displays the CMD of stars in Mon OB1, Mon R1, and the halo. The G RP magnitudes of members were corrected by the distance moduli of their host groups (see Section 3). We superimposed four isochrones (gray curves) for 2, 3, 5, and 10 Myr (Dotter 2016;Choi et al. 2016) on the CMD. A mean total extinction (A V ) of 0.22 mag in visual magnitude (Sung et al. 1997) was applied to the isochrones. The faint part of the CMD is significantly affected by some factors such as large photometric errors, high internal extinction of disk-bearing stars, and variabilities Stauffer et al. 2014;Lim et al. 2015;Stauffer et al. 2016, etc). In addition, the systematic difference between the isochrones from the adopted evolutionary models and the observed CMD is found for late-type pre-main sequence stars. For these reasons, we only considered the bright part of the CMD as seen in Figure 14.

Ages of stellar groups
The magnitude (or luminosity) of the main sequence turn-on is sensitive to the age of a given coeval stellar system. The ages of members in Mon OB1 roughly range from 2 to 10 Myr. An age spread of 4 -5 Myr has been expected from the lithium abundances of pre-main sequence stars (Lim et al. 2016) and from the CMD analysis of Venuti et al. (2018). Hence, star formation in Mon OB1 might be sustained on a several Myr scale. The S Mon group (2 -3 Myr) seems to be younger than the Cone and THF15 groups ( 5 Myr) given that main sequence turn-on appears at higher luminosity. The members of Mon R1 have ages similar to those of the S Mon group. On the other hand, the CMD of the halo stars overlaps with that of the older populations of Mon OB1 (Cone and THF15).
The photometric errors of stars around the main sequence turn-on (G RP < 12 mag) are much smaller than 0.01 mag. Therefore, these are not the major source of uncertainties in age estimation. The error on distance (±40 pc) corresponds to ±0.1 mag in distance modulus. The contribution of this error to age estimation is also negligible. The other factor is differential reddening across the survey region.
The differential reddening over Mon OB1 is known to be small E(B − V ) = 0.07 ± 0.03 (Sung et al. 1997). Sung et al. (2008) found higher reddening values of E(B − V ) ∼ 0.2 for low-mass pre-main sequence stars in the region (see also Rebull et al. 2002). We plotted the isochrones reddened by this high reddening value (gray dashed curves) in Figure 14. But, the colors and magnitudes of main sequence stars are closer to those of the isochrones reddened by the mean reddening value of E(B − V ) = 0.07 (A V = 0.22 mag). This implies that the differential reddening is not high enough to affect the relative ages among the stellar groups in Mon OB1.
Most of the members of Mon R1 seen in Figure 14 are the bright members of IC 447. The reddening of these stars due to the intracluster medium may be small. Indeed, they were found in the cavity of the dusty cloud (see figure 1 of Bhadari et al. 2020). The B-type member with the bluest color in IC 447 has a color similar to those of members in Mon OB1. Therefore, the minimum reddening toward Mon R1 may be similar to the mean reddening of Mon OB1. There is a bright member of IC 446 in the CMD (G RP − DM = 1.24, G BP − G RP = 0.90). The age of the star is about 2 Myr, which is comparable to the ages of the IC 447 members. The ages of the bright stars in IC 446 and IC 447 may not be significantly altered ( 1 Myr) by differential reddening given the reddening vector.

A Large-scale Distribution of Stars and Gas
We plotted the IRIS image at 100µm and the AKARI false-color (blue : 65 µm, green : 90 µm, and red :140 µm) image of interstellar material over our survey region in Figure 15. Interestingly, the IRIS image reveals a large arc-like structure across the survey region. The members of Mon OB1, Mon R1, and the halo were superimposed on the AKARI image. It seems clear that star formation is actively taking place along the arc-like structure. Indeed, bright knots notably host Mon OB and Mon R1, but also G202.3+2.5 (Montillaud et al. 2019a,b) and some continuum sources that are bright at infrared and submillimeter wavelengths (Beichman et al. 1988;Di Francesco et al. 2008). Figure 15 also displays the PM vectors of members relative to the systemic motion of Mon OB1. A high fraction of the halo stars are found around Mon OB1 and tend to move outward from the association. Their Φ distribution has a strong peak at around 20 • . On the other hand, the members of Mon R1 are systematically moving toward south relative to Mon OB1.

Implication on cluster formation
The S Mon and Cone groups show patterns of expansion as seen in the other stellar clusters (Cantat-Gaudin et al. 2019;Kuhn et al. 2019;Lim et al. 2019aLim et al. , 2020Lim et al. , 2021. About 50% of members in the S Mon group are radially escaping from this group, but some members beyond 5 (Φ ∼ ±180 • , Figure 10) are still sinking into the group center. The fraction of these members may be less than 20% from the last two bins around Φ ∼ ±180 • in the histogram of Figure 10. The Cone group shows a similar pattern, but less clear than that of the S Mon group. Such a trend was also found in the Orion Nebula Cluster (Platais et al. 2020). On the other hand, the young open clusters IC 1805 and NGC 2244 only show patterns of expansion without signature of collapse (Lim et al. 2020(Lim et al. , 2021. The different internal kinematics among these clusters may result from the different initial conditions of cluster formation and evolution time. Many theoretical studies have tried to explain the expansion of stellar clusters as the result of their dynamical evolution after rapid gas expulsion (Tutukov 1978;Hills 1980;Lada et al. 1984;Kroupa et al. 2001;Banerjee & Kroupa 2013, 2015. Our previous study (Lim et al. 2020) (Beichman et al. 1988;Di Francesco et al. 2008) in the northern knot were marked by purple dots. The colors of the other symbols are the same as in Figures 5 and 7. formation and evolution of these two groups. In addition, the rapid gas ejection and stellar feedback could affect the structure and dynamics of these groups (Gavagnin et al. 2017).
The S Mon and Cone groups also show the signature of rotation. Some young stellar clusters were also found to be rotating, e.g., R136 (Hénault-Brunet et al. 2012), Trumper 15 (Kuhn et al. 2019), and NGC 2244 (Lim et al. 2021). These results provide important constraints on the cluster formation process, such as the monolithic collapse of rotating clouds or the hierarchical assembly of subclusters (Corsaro et al. 2017;Mapelli 2017). A large fraction of molecular clouds in external galaxies were found to be rotating (Rosolowsky et al. 2003;Tasker 2011;Braine et al. 2020). Rotating clusters can naturally form in such rotating clouds. However, the groups in Mon OB1 have different directions of rotation, which cannot be explained by only the monolithic collapse of a single molecular cloud.
Theoretically, collisions between molecular clouds can result in a larger cloud with retrograde rotation with respect to the galactic rotation on a large spatial scale (Dobbs et al. 2011). If such collisions could occur between gaseous and stellar clumps on small spatial scales (several pc), then the angular momentum vector of the natal cloud could be changed (Mapelli 2017). In addition, the observational evidence for an on-going merger of stellar clusters was reported (Sabbi et al. 2012). Therefore, our findings suggest that at least one of the groups in Mon OB1 might have been formed via the hierarchical merging process.
On the other hand, the stellar groups in Mon R1 do not seem to form as self-gravitating systems given the weak central concentration of stars. As these groups expand and move away from each other, their members will eventually disperse to become field star population in the Galactic disk (Miller & Scalo 1978;Briceño et al. 2007).

Star formation on different spatial scales
The projected distance and the line of sight distance between Mon OB1 and Mon R1 are about 20 and 40 pc, respectively, which is the typical size of giant molecular clouds. Members in these two associations have velocities in almost the same range (see Figures 6 and 8). As a comparison point, note that the Orion A cloud shows an RV variation larger than 10 km s −1 from north to south (Yun et al. 2021). In addition, molecular gas constituting the arc-like structure as seen in Figure 15 was found in the same velocity channel (0 < V [km/s] < 15 from figure 2-b of Oliver et al. 1996). Hence, the two associations and the other small SFRs related to bright knots have formed in the same molecular cloud. However, it is unclear how their formation is related to each other.
Based on the age distributions of pre-main sequence stars in Mon OB1, Sung & Bessell (2010) proposed an outside-in star formation history for Mon OB1. Star for- mation has been initiated from the halo and propagated via the S Mon and Cone groups to the Cone(C) and Spokes groups around the embedded YSOs NGC 2264 IRS 1 and IRS 2 (see also Venuti et al. 2018). However, this scenario should be slightly modified, according to our results. Indeed, star formation initiated in the southern region (Cone and THF15) and then occurred in the northern region (S Mon). Recent star-forming activity has ignited in the molecular gas behind the southern region, which may have resulted in the formation of the Cone(C) and Spokes groups. The star formation history of Mon OB1 and Mon R1 is illustrated in a simple cartoon (Figure 16).
A number of halo stars were found around Mon OB1, and they seem to be escaping from the association (Figure 15). Most of them might thus have formed in Mon OB1 as the first generation of stars. The fact that the halo stars have almost the same ages as those of the older populations in Mon OB1 supports this argument. Therefore, the eastern halo might have been formed by stars escaping from the association rather than the star formation in-situ proposed by Sung & Bessell (2010).
The systemic motion of Mon R1 may provide clues to the formation of this association. If this association had been formed in the compressed clouds by feedback from the O-type binary S Monocerotis which is located to the east of Mon R1, its members would move westward. However, they do not show such a systematic motion toward west ( Figure 15). Furthermore, there is no significant age difference between the S Mon group and Mon R1. The stellar groups in Mon R1 might thus have spontaneously formed at about the same epoch as the formation of the S Mon group in Mon OB1.
Herschel images at submillimeter wavelengths reveal the networks of filaments in many molecular clouds (André 2015). The presence of filamentary structures in molecular clouds has been accepted as a ubiquitous feature. Turbulence could play a key role in the for-mation of filamentary structures (Larson 1981;Padoan et al. 2001). Recently, there is increasing evidence that magnetic fields significantly contribute to the formation of filaments on a small scales (Wang et al. 2019;Doi et al. 2020, etc). Cores and protostars form after the gravitational fragmentation of filaments (André et al. 2010). Filament hubs with high densities are the sites of stellar cluster formation (Schneider et al. 2012;Galván-Madrid et al. 2013;Treviño-Morales et al. 2019). The relics of such filaments have been found toward the Vela OB2 association and the Orion region (Jerabkova et al. 2019;Beccari et al. 2020;Pang et al. 2021). In our survey region, SFRs including Mon OB1 and Mon R1 are distributed along a large arc-like structure in a hierarchical way. The formation of structures on different spatial scales by the actions of turbulence, gravity, and magnetic field may be the essential formation process of Mon OB1 and Mon R1.

SUMMARY
In this study, we investigated the spatial distributions and kinematic properties of young stellar population in the two stellar associations Mon OB1 and Mon R1 to understand star formation process and their physical association.
We first isolated member candidates in a 6 • ×6 • survey region using the published data sets. Then, a total of 728 members were finally selected from the criteria based on the Gaia parallaxes and PMs. The spatial distributions of these stars show substructures that are kinematically distinct. Mon OB1 contains three optically visible stellar groups, the S Mon, Cone, and THF15 groups. We also suggested the possibility that there are two embedded groups (Spokes and Cone(C)) behind the optically visible groups. Mon R1 hosts the open cluster IC 447 and two partially embedded groups (N2245/47 and IC 446). In addition, some stars were found in the halo region.
The stellar groups, except for THF15 show patterns of expansion as seen in many associations. In addition, the signature of rotation was detected for the S Mon and Cone groups. Interestingly, these groups are rotating in opposite directions, which could be a trace if clouds having merged in the past.
We analyzed the CMD of members to infer the star formation history in the survey region. The members of Mon OB1 have ages ranging 2 Myr to 5 Myr. The ages of the Mon R1 members (2 -3 Myr) are similar to those of the younger population in Mon OB1, while the halo stars ( 5 Myr) have ages similar to those of the older population. Furthermore, the motions in Mon R1 are not pointing away from Mon OB1. This suggests that Mon OB1 and Mon R1 might have formed independently in a giant molecular cloud.
In addition, Mon OB1 and Mon R1 belong to a large scale arc-like structure comparable to the size of typical giant molecular clouds. Actually, more star formation activities on small scales are found along this large structure, forming a hierarchy: isolated stars, clusters, and associations. Hence, these two associations might have formed within the same cloud in a hierarchical way. In addition, the expansion of stellar groups plays a crucial role in the formation of the halo population. This paper has made use of data obtained under the K-GMT Science Program (PID: GEMINI-KR-2020A-003 and Gemini program number: GS-2020A-Q-239) funded through Korean GMT Project operated by Korea Astronomy and Space Science Institute (KASI) and from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has also made use of the SIMBAD database,operated at CDS, Strasbourg, France. Based on observations obtained at the international Gemini Observatory, a program of NSF's NOIRLab [ include additional acknowledgment here, see section below ], which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (KASI) (Republic of Korea). This work used the Immersion Grating Infrared Spectrometer (IGRINS) that was developed under a collaboration between the University of Texas at Austin and the KASI with the financial support of the US National Science Foundation under grants AST-1229522 and AST-1702267, of the University of Texas at Austin, and of the Korean GMT Project of KASI. This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (Grant No : NRF-2019R1C1C1005224 and 2022R1C1C2004102). Y.N. acknowledges support from the Fonds National de la Recherche Scientifique (Belgium), the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Programme (contracts linked to XMM and Gaia).