A Kinematic Perspective on the Formation Process of the Stellar Groups in the Rosette Nebula

Our targets are young stars within a 1° × 1° region centered at R.A. = 06^h 31^m 5500, decl. = + 04° 300 (J2000). This SFR (b ∼ − 2072) is located in the Galactic plane, and therefore a large number of field interlopers are also observed in the same field of view. We intend to isolate member candidates first using multiple sets of archival data in this section, and then the final decision on membership will be made by using the recent astrometric data from the Gaia Early Data Release 3 (EDR3; Gaia Collaboration et al. 2021) in Section 3.

To identify probable member candidates in the survey region, we considered the intrinsic properties of young stars. O- and B-types stars spread over this SFR are probable member candidates because the positions of these stars are close to birthplaces within their short lifetime, particularly for O-type stars. We took a list of OB stars from the data bases of MK classifications (Mermilliod & Paunzen 2003; Reed 2003; Skiff 2009; Maíz Apellániz et al. 2013, SIMBAD) and combined them into a list after checking for duplicates. A total of 112 OB stars were selected as member candidates.

A large number of low-mass YSOs have a warm circumstellar disk. These disk-bearing YSOs appear bright in infrared passbands (Lada 1987; Gutermuth et al. 2008; Koenig et al. 2012). YSO candidates were identified using the published Spitzer and AllWISE data (Balog et al. 2007; Cutri et al. 2013); see Appendix appendix for details. We found 291 YSO candidates (10 Class I, 262 Class II, and 19 YSOs with a transitional disk) from the Spitzer data in total, while a total of 84 candidates (six Class I, 76 Class II, and two objects with a transitional disk) were identified from the AllWISE data. There are 37 YSO candidates in common between these data sets. Note that the survey region () of Spitzer (Balog et al. 2007) is smaller than ours. Our list contains 338 YSO candidates in total.

We used photometric data from the Gaia Data Release 2 (DR2; Gaia Collaboration et al. 2018) instead of the recent EDR3 (Gaia Collaboration et al. 2021) because the reddening law is empirically well-established (Wang & Chen 2019) for the former. We found counterparts for all the OB stars and for 320 YSO candidates in the Gaia DR2 (Gaia Collaboration et al. 2018). In addition, we searched for counterparts of X-ray sources detected in this SFR (Wang et al. 2008, 2009, 2010). Stars found within searching radii of 10 and 15 were classified as X-ray sources and X-ray source candidates, respectively. A total of 720 X-ray sources and 134 candidates were found in the Gaia data. Note that the lists of X-ray sources are complete down to 0.5–1 M_☉ (Wang et al. 2008, 2009, 2010).

Figure 1 displays the color–magnitude diagram (CMD) of stars in the Rosette Nebula. Most member candidates seem to occupy a specific locus in the diagram. The main-sequence turn-on is found at about 13 mag in G_RP. There are several B-type stars fainter than YSOs near the main-sequence turn-on. These stars may be background early-type stars. A large fraction of the X-ray sources may be YSO candidates, while some X-ray sources appear to fall on the main-sequence band. These X-ray emitting stars may be late-type field stars. All the YSO candidates with infrared excess emission were considered as member candidates because their colors and magnitudes could vary with different levels of internal extinction by material in disks or envelopes and accretion activities (Lee et al. 2020). Photometric errors increase with brightness, and therefore we limited our sample to 840 candidates brighter than 18 mag. The selection criteria of member candidates can be summarized as below:

• 1.  
  G_RP < 18 mag.
• 2.  
  O- or B-type stars brighter than 13.5 mag in G_RP.
• 3.  
  Class I YSOs.
• 4.  
  Class II YSOs.
• 5.  
  YSOs with a transitional disk.
• 6.  
  X-ray sources within the locus (see Figure 1).

Zoom In Zoom Out Reset image size

We performed spectroscopic observations of 316 YSO candidates with or without an X-ray emission and ionized gas on 2019 October 13, November 9, 13, and 18 using the high-resolution (R ∼ 34,000) multi-object spectrograph Hectochelle (Szentgyorgyi et al. 2011) on the 6.5 m telescope of the MMT observatory. The spectra of the YSO candidates were taken with the RV31 filter (5150–5300 Å) in 2 × 2 binning mode to achieve good signal-to-noise ratios. In addition, several tens of fibers were simultaneously assigned to blank sky to obtain sky spectra. The spectra of ionized gas were obtained from the fibers assigned to 238 positions on the Rosette Nebula (see Figure 1). The order separating filter OB25 (6475–6630 Å) was used to observe the forbidden line [N ii] λ6584. For calibration, dome flat and ThAr lamp spectra were also obtained just before and after the target observation. We present a summary of our observations in Table 1.

Table 1. Summary of Observations

	Date	Setup	Filter	Na	Exposure time	Binning	Seeing
	2019 Oct 13	1	RV31	111	35 minutes × 3	2 × 2	09
	2019 Nov 13	2	RV31	110	40 minutes × 3	2 × 2	14
YSOs					38 minutes × 1
		3	RV31	95	37 minutes × 1	2 × 2	13
					35 minutes × 1
	2019 Nov 18	3	RV31	95	35 minutes × 2	2 × 2	13
Ionized gas	2019 Nov 9	4	OB25	238	15 minutes × 5	1 × 1	09

Note.

^a N represents the number of observed targets.

Download table as:  ASCIITypeset image

We preprocessed the raw mosaic frames using the IRAF ⁶ /MSCRED packages in a standard way. One-dimensional spectra were then extracted from the reduced frames using the dofiber task in the IRAF/SPECRED package. Target spectra were flattened using dome flat spectra. The solutions for the wavelength calibration obtained from ThAr spectra were applied to the target and sky spectra.

Our observations were performed under bright sky condition, and therefore the scattered light was unevenly illuminated over the field of view (1° in diameter), resulting in a spatial variation of sky levels. For a given observing setup, we constructed a map of the sky levels and combined sky spectra into a master sky spectrum with an improved signal-to-noise ratio. The sky level at a target position was inferred by interpolating the position of the target to the map, and its sky spectrum was obtained from the master sky spectrum scaled to the inferred sky level. For each target spectrum, we subtracted the corresponding scaled master sky spectra and then combined sky-subtracted spectra into a single spectrum for the same target observed on the same night. Finally, the target spectra were normalized by using continuum levels traced from a cubic spline interpolation.

Since the spectra of low-mass stars contain a large number of metallic lines, their RVs can be measured more precisely than those of high-mass stars. The RVs of the YSO candidates were measured by applying a cross-correlation technique to the observed spectra. Assuming the solar chemical abundance, a total of 31 synthetic spectra in a wide temperature range of 3800–9880 K were generated from the MOOG code and Kurucz ODFNEW model stellar atmosphere (Sneden 1973; Castelli & Kurucz 2004). We derived the cross-correlation function (CCF) between the observed spectra and synthetic ones using the xcsao task in the RVSAO package (Kurtz & Mink 1998) and adopted the central velocities at the strongest CCF peaks as RVs. The errors on RVs were estimated from the equation , where w, h, and σ_a represent the full widths at half maximum of CCFs, their amplitudes, and the rms from antisymmetric components, respectively (Tonry & Davis 1979; Kurtz & Mink 1998). The measured RVs were converted to velocities in the local standard of rest frame using the IRAF/RVCORRECT task.

The spectra of 95 YSO candidates (Setup 3 in Table 1) were obtained on two different nights. The RVs of these stars were measured for each night. The median difference between two RV measurements for 52 common stars is about 0.3 km s⁻¹, which is smaller than the mean of measurement errors (∼1.5 km s⁻¹). The weighted means of two measurements were adopted as the RVs of given stars, where the inverse of the squared error was used as the weight value. We measured the RVs of 224 out of 316 YSO candidates in total. The spectra of the other 92 stars have insufficient signals to derive CCFs, or some of them show only emission lines.

We used the forbidden line [N ii] λ6584 to probe the kinematics of ionized gas. In general, the critical density of this line is higher than the typical electron densities of the Galactic H ii regions (Copetti et al. 2000), and photons emitted from the singly ionized nitrogen atoms are not absorbed along the line of sight. Therefore, this emission line can be used to trace the structure and kinematics of ionized gas distributed across an H ii region. The RVs of the ionized gas were measured from the line center of the best-fit Gaussian profile. We fit multiple Gaussian profiles to some complex line profiles and measured the RVs of the multiple components along the line of sight.

We used radio data to investigate the physical association between stellar groups and remaining molecular gas. The ¹²CO and ¹³CO (J = 1 − 0) line maps were obtained from Heyer et al. (2006). These radio maps covering both the Rosette Nebula and the eastern molecular cloud are useful for investigating the velocity fields of remaining molecular gas at the boundary of the H ii region. In addition, the excitation temperature T_ex and column densities can be estimated by assuming of local thermodynamic equilibrium.

Since the ¹²CO line is, in general, optically thick, T_ex can be derived from the following equation by Pineda et al. (2008):

where is the brightness temperature at the peak intensity of ¹²CO. The column density of ¹³CO was estimated from the integrated intensity along the line of sight (in units of K km s⁻¹) assuming an optically thin case (Pineda et al. 2008). On the other hand, the column density of ¹²CO was estimated using the relation of Scoville et al. (1986). However, this relation considers an optically thin case. The effects of optical depth were corrected by adopting the correction factor addressed in Schnee et al. (2007).

The integrated intensity maps of molecular gas were also constructed from the ¹²CO (J = 1 − 0) cube data. Note that the ¹²CO lines between velocities of −5 to 25 km s⁻¹ were integrated. In addition, position–velocity diagrams were obtained by integrating the data cube along R.A. and decl.

Figure 5 displays the radial surface density profile fit to the model profiles of Elson et al. (1987) in a logarithmic scale. The surface density of NGC 2244 decreases with the radial distance, but is influenced by NGC 2237 beyond the radius of (see also the left panel of Figure 4). Therefore, we adopted this boundary as the apparent radius of NGC 2244 (r_cl). This radius encompasses the majority of cluster members. The fundamental parameters and kinematics of this cluster were investigated within this radius.

Zoom In Zoom Out Reset image size

4.1.1. Age

In order to derive the fundamental parameters (reddening, age, mass, and half-mass–radius) of this cluster, we fit the isochrones from MESA models considering the effects of stellar rotation (Choi et al. 2016; Dotter 2016) to the CMD. The distance modulus of 10.73 mag (equivalent to 1.4 kpc) was applied to the unreddened isochrones (A_V = 0 mag). For the early-type members with colors bluer than G_BP − G_RP = 0.6, the reddening of the individual stars was obtained by fitting them to the distance-corrected isochrones along the reddening vector of Wang & Chen (2019) [A_GRP/E(G_BP − G_RP) = 1.429]. The mean reddening is then estimated to be 〈E(G_BP − G_RP)〉 = 0.733 ± 0.078.

Figure 6 shows the CMD of the members in NGC 2244. The isochrones for four different ages (2, 3, 5, and 10 Myr) were superimposed on the CMD after correction for the mean reddening. Some B-type stars may be already evolving into the main sequence. The magnitude of the main-sequence turn-on is sensitive to the age of the cluster. The isochrone for 2 Myr well-matches the main-sequence turn-on. We adopt 2 Myr as the age of NGC 2244, and this is in good agreement with the previous estimates (Hensberge et al. 2000; Park & Sung 2002). However, the entire shape of the CMD seems to be fit by isochrones in a wide range of ages (ΔAge > 7 Myr). There could be systematic uncertainties in the calibration adopted in the theoretical models (Choi et al. 2016; Dotter 2016). Observationally, photometric errors, imperfect reddening correction, binaries, and variabilities of YSOs are possible sources of the observed scatter in the CMD (Lim et al. 2015). In addition, the luminosity spread can be caused by multiple star formation events.

Zoom In Zoom Out Reset image size

4.1.2. The Initial Mass Function and Total Stellar Mass

The mass function is important to estimate the total cluster mass, which allows us to infer the virial state of NGC 2244 at a current epoch. To derive the mass function of this cluster, we first obtained the masses of individual stars from the CMD by means of the isochrones from the MESA evolutionary models (Choi et al. 2016; Dotter 2016). The masses of early-type main-sequence members (G_BP − G_RP < 0.6) were obtained by interpolating their magnitude to the mass–luminosity relation constructed from the isochrones for 2, 5, and 10 Myr, where the individual reddening values were applied to the isochrones. For the other members, we adopted a mean reddening because the individual reddening values of these stars cannot be obtained from the current data. The colors of stars were interpolated to the grid of the isochrones with different ages (0.1–40 Myr) in order to construct the mass–luminosity relations at the given colors. The masses of individual members were then obtained from the interpolation of their magnitude to the mass–luminosity relations.

To compute the uncertainty in mass estimation, the photometric errors were added to or taken out from the observed magnitudes and colors. We then found the masses associated to these increased or decreased magnitudes and colors in the same way as done before. The mass uncertainties were adopted from half the difference between the upper and lower values. A typical uncertainty due to photometric errors is about 0.1M_⊙. Differential reddening can also affect the mass estimation. The uncertainty in mass is then about 0.3 M_⊙, on average, if we considered the differential reddening of 0.078.

We derived the present-day mass function of this cluster by counting the number of stars in given logarithmic mass bins (). A bin size of 0.2 was used to count the number of stars with mass smaller than 10 M_⊙, while a larger bin of 1.08 was used for higher-mass stars because of their paucity. Then, the mass function was normalized by the associated bin sizes. To avoid the binning effects, we repeated the same procedure by shifting the bin size by half. The Poisson noise was adopted as the errors on the mass function.

Figure 7 displays the present-day mass function. The difference between the initial and present stellar masses is negligible in the first few Myr. According to the MESA model (Choi et al. 2016; Dotter 2016), the most massive star in NGC 2244 has lost 4M_⊙ for 2 Myr. Hence, the present-day mass function can be considered as the initial mass function. We constrained the slope (Γ) of the mass function down to the turnover mass (∼1 M). The slope is estimated to be Γ = −1.4 ± 0.1 from a least-squares fitting method for stars with masses larger than 1 M_⊙. This is in good agreement with those obtained from other SFRs in the solar neighborhood (Salpeter 1955; Kroupa 2001).

Zoom In Zoom Out Reset image size

However, it does not guarantee that our member selection is complete down to the turnover mass. Since we selected the most likely members using strict selection criteria, many lower-mass members with large errors in parallax and PM may have been excluded. We estimated the upper limits of the mass function including member candidates fainter than 14 mag in G_RP band in the same way as above (open circles in Figure 7). The upper limits follow a trend similar to the original one for masses larger than 1 M_⊙ and become larger for the lower-mass regime. The number of stars obtained from the original mass function is about 30% lower than that from the upper limits at around the turnover mass. However, this may not be the actual completeness limit of our sample, because some field interlopers may be included in the sample of member candidates. If we adopt the upper limits of the mass function, then the slope of the mass function is estimated to be Γ = −1.6 ± 0.2, which is consistent with the value found above within errors.

The sum of the stellar masses yields a total mass of 533 M_⊙. However, this result does not consider a number of stars with very low mass (<1 M_⊙). In order to infer their population, we scaled the Kroupa initial mass function to ours (the original mass function) using the mean difference between these mass functions. This scaled Kroupa initial mass function was then integrated in the mass range of 0.08 M_⊙ to 1 M_⊙. The number of cluster members down to the hydrogen-burning limit is expected to be about , where the upper and lower limits were inferred from the integration of the scaled Kroupa initial mass functions adjusted by the uncertainty of the scaling factors (i.e., the standard deviation) in the low-mass regime (<1 M_⊙). The total stellar mass of NGC 2244 is estimated to be about 879 M_⊙, which is slightly larger than those previously obtained by other studies (770M_⊙—Pérez 1991; 625M_⊙—Bonatto & Bica 2009).

We also estimated a projected half-mass radius (r_h) of (equivalent to 2.4 pc) from our limited sample. A recent study derived the initial mass function of stars in a very small central region, down to the subsolar mass regime (Mužić et al. 2019). The mass function appeared similar to the Kroupa initial mass function, which implies that there is no sign of dynamical mass segregation. Hence, r_h obtained in this study may not be significantly altered even if the radial mass distribution of stars in the full mass range is considered.

4.1.3. Internal Kinematics

The PM vectors of the members relative to the systemic motion of NGC 2244 are shown in Figure 4. Stars in the central cluster show outward motions, which were also detected in Kuhn et al. (2019). To probe the direction of their motion in detail, we calculated the vectorial angle (Φ) between the position vectors from the cluster center and the relative PM vectors of stars (see also Lim et al. 2019, 2020). A Φ close to 0° means that a given star is radially moving away from the center of the cluster, while a Φ of 180° indicates that the star is moving toward the cluster center.

Figure 8 displays the Φ distribution of stars in NGC 2244. The cluster members have Φ values in a wide range. We present the histograms of Φ values in the right panel of the figure. There is a peak at 0°, which indicates that some stars are radially escaping from NGC 2244. Given that this feature seems weaker within r_h (blue open squares in the figure), the majority of these escaping stars may be distributed beyond r_h.

Zoom In Zoom Out Reset image size

In order to search for the signature of cluster rotation, the RV distribution of stars was investigated as addressed in many previous studies (Lane et al. 2009; Mackey et al. 2013). We first considered a projected rotational axis across the cluster center from north (0°) to south (180°). The difference between the mean RVs of stars in the two regions separated by the projected rotational axis was computed. The same procedure was repeated for various position angles of the projected rotational axis with an interval of 20°.

We found a variation of the mean differences of RVs with the position angles of the projected rotational axis, as shown in the left panel of Figure 9. This variation was fit to the sinusoidal curve as below:

where V_rot, i, and ϕ₀ represent the rotational velocity, inclination of a rotational axis, and phase, respectively. The amplitude of the sinusoidal curve corresponds to twice the projected rotational velocity (). The position angle of the projected rotational axis can be inferred from 270° − ϕ₀. The projected rotational axis is tilted by ∼20° (±14°) from north to east. We investigated a variation of using stars within various radii from to . The middle panel of Figure 9 displays the variation of with the projected radial distances. In the inner region (<65), of stars increases up to 0.67 km s⁻¹ with the radial distance, and drops to ∼0.20 km s⁻¹ beyond this radius. It is interesting that this cluster rotation is evident around r_h. A clear pattern of cluster rotation is also seen in the position–velocity diagram (the right panel of Figure 9).

Zoom In Zoom Out Reset image size

The signature of cluster rotation was also searched for using PMs. The one-dimensional tangential velocities (V_t) of stars along R.A. and decl. were calculated by adopting a distance of 1.4 kpc. Since the diameter of the cluster is very small compared to the distance, the error on V_t by the extent of the cluster is less than 1% under the assumption that this cluster has a spherical shape. The rotational velocity on the celestial sphere is given by the equation below:

where r, d_R.A., V_R.A, d_decl., and V_decl. represent the radial distance, the position and V_t in R.A., and the position and V_t in decl., respectively. The nonzero mean velocity would yield the rotational velocity of this cluster. However, the mean velocity was calculated to be 0.03 ± 1.32 km s⁻¹, meaning that this cluster does not have significant azimuthal motion. This result is consistent with that found by Kuhn et al. (2019). Since the cluster rotation was only detected in RV, the rotational axis of NGC 2244 may be tilted by almost 90° from the line of sight.

4.1.4. Dynamical State

We computed the virial velocity dispersion of NGC 2244 by using the equation (Parker & Wright 2016) below:

where G, M_total, R, and η represent the gravitational constant, enclosed total mass, radius, and the structure parameter, respectively. Note that we can consider the total cluster mass of 879 M_⊙ as the enclosed mass, given that gas has already been removed around NGC 2244 according to ¹²CO (J = 1 − 0) and ¹³CO (J = 1 − 0) observations (Heyer et al. 2006). The r_h of 2.4 pc was used in Equation (4). The structure parameter η has a value between 1 and 12 depending on the geometry of the cluster. Here, η of 10 was adopted because the observed surface density profile (γ = 3.5 ± 0.3; Figure 5) is close to that (γ = 4) of Plummer (1911; see also Figure 4 of Portegies Zwart et al. 2010). The expected velocity dispersion at a virial state is about 0.6 km s⁻¹ (for reference, the virial velocity dispersion for other clusters is in the range of 0.4–1.6 km s⁻¹; see Kuhn et al. 2019). The virial velocity dispersion is not very sensitively varied with the total mass, because it is proportional to . The virial velocity dispersion is almost constant even if we adopted the upper limit of the total cluster mass.

Figure 10 displays the distributions of V_R.A, V_decl., and RV (the upper panels) with their error distributions (the lower panels). We derived the systemic velocities and velocity dispersions from the best-fit Gaussian distributions (Table 3). The typical errors on the measurements are estimated from the weighted-mean errors, where the associated error distributions were used as the weight values for given errors. We obtained the intrinsic velocity dispersions of 0.8 and 0.8 km s⁻¹ along R.A. and decl., respectively, from quadratic subtraction between the observed velocity dispersions and their typical errors. The effect of system rotation on the observed RV dispersion can be simulated by using the ideal Gaussian velocity distribution based on the underlying rotation with position angles. As a result, the observed velocity dispersion can be inflated by about 14% compared to the intrinsic velocity dispersion. Since our measurement of the velocity dispersion along the line of sight, after quadratic subtraction by the typical error of 1.6 km s⁻¹, is about 0.8 km s⁻¹, the intrinsic RV dispersion would be 0.7 km s⁻¹. The mean value from the intrinsic velocity dispersions along R.A., decl., and the line of sight is 0.8 ± 0.1 km s⁻¹.

Zoom In Zoom Out Reset image size

Table 3. Velocity Dispersions

	Vsys	σobs	σerr	σint
	(km s−1)	(km s−1)	(km s−1)	(km s−1)
VR.A	−11.4	1.0	0.6	0.8
Vdecl.	2.0	1.0	0.6	0.8
RV	12.8	1.8	1.6	0.7 a

Notes. Velocity dispersions were obtained for members within r_h . Columns (2)–(4) represent the systemic velocities, the observed velocity dispersions, typical errors, and the intrinsic velocity dispersions along R.A, decl., and the line of sight, respectively.

^a The contribution of the cluster rotation (14%) has been corrected for the RV dispersion.

Download table as:  ASCIITypeset image

The initial size of this cluster might have been smaller than the current one. In addition, a large amount of the natal cloud mass probably remained around the cluster. The virial velocity dispersion might have thus been larger than what we derived. While the velocities of cluster members at a supervirial state are not significantly altered in the first several Myr (Schoettler et al. 2019), the virial velocity dispersion could decrease with the evolution of the cluster because of the cluster expansion and gas expulsion. At the current epoch, the observed velocity dispersion is comparable to times the virial velocity dispersion, which implies that the total kinetic energy is almost in balance with the total potential energy. This cluster may be gravitationally unbound. The weak pattern of expansion within r_h supports this argument.

The relaxation time of this cluster was estimated from the equation of Binney & Tremaine (1987). The crossing time of stars for r_h is about 2.9 Myr, adopting the mean velocity dispersion of 0.8 km s⁻¹. Given the total number of stars (1510), the dynamical relaxation time was estimated to be about 76 Myr. The age of NGC 2244 is much younger than this timescale, and therefore dynamical mass segregation may not be found in this cluster. Figure 11 displays the two-dimensional V_t () distribution with respect to the radial distance and stellar mass. Since dynamical mass segregation is the result of energy equipartition (Binney & Tremaine 1987), the signature of this dynamical process can be probed from a correlation between stellar mass and velocities, i.e., the low-mass stars have velocities higher than those of high-mass stars. However, there is no evidence for the energy equipartition among cluster members, confirming the result of Chen et al. (2007). Furthermore, the most massive O-type stars are widely distributed across the cluster region.

Zoom In Zoom Out Reset image size

The physical association between stellar groups and remaining gas can be a key clue to understanding the formation of substructures. Toward this aim, the distributions of stars and gas were investigated in position–velocity space. We constructed a three-dimensional position–position–velocity diagram of ionized gas as made in our previous studies (Lim et al. 2018, 2019). A data cube composed of 90 × 90 × 90 regular volume cells was created, and then a Delauney triangulation technique was used to interpolate the intensities and velocities of the forbidden line [N ii] λ6584 into the individual cells. From this data cube, we also obtained position–velocity diagrams by the sum of the counts along R.A. and decl., respectively. The positions and RVs of stars are overplotted in these diagrams.

Figure 12 shows the distribution of stars and ionized gas in position–velocity space. A cavity of the H ii region is clearly seen in the central region occupied by NGC 2244 (red). In the position–velocity diagrams, the other stellar groups (blue) are found along the H ii bubble, and they have velocities roughly similar to those of the ionized gas. This indicates that these groups and remaining gas belong to the same physical system. The RVs of ionized gas range from 0 to 30 km s⁻¹, and therefore the Rosette Nebula is expanding at velocities up to 17.2 km s⁻¹ from NGC 2244 (12.8 km s⁻¹; Table 3).

Zoom In Zoom Out Reset image size

In order to find the physical association between stellar groups and gas, it is necessary to probe them in detail in smaller position–velocity space. We investigated the three stellar groups that are likely associated with adjacent gas structures (the left panel of Figure 4). An amount of molecular gas still remains in the vicinity of these groups. Figure 13 exhibits the distribution of stars on the integrated intensity maps and position–velocity diagrams of molecular gas. The RVs of ionized gas measured from the [N ii] λ6584 were also superimposed on the figure. The ionized and molecular gas appear to have similar RVs.

Zoom In Zoom Out Reset image size

Most stars in subregion (a) may be members of the neighboring cluster NGC 2237 (Li 2005; Wang et al. 2010). There are two gas components in the position–velocity diagrams, which is indicative of expanding shell components. The ionized gas also traces both components in RV. The heads of the gas pillars seen in the left panel of Figure 4 correspond to the molecular gas clumps at (ΔR.A., Δdecl.) ∼(,5). The RVs of these clumps range 0–5 km s⁻¹. Given the systemic velocity of NGC 2244 (12.8 km s⁻¹), they may be part of the near side of the H ii bubble. Some gas structures associated with this bubble extend southwest. This component may be expanding at velocities up to −12.8 km s⁻¹, although this is a projected expanding velocity. The other gas component is found at (ΔR.A., Δdecl.) ∼(,1) and has RVs in a range of 15–20 km s⁻¹. This component may be part of the bubble receding away from NGC 2244 at projected velocities up to 7.2 km s⁻¹. NGC 2237 is located in the vicinity of the heads of the gas pillars. Most stars in this group tend to systematically move toward the south. Stars with RV measurements have RVs in range 15–21 km, except for one star. Therefore, this group seems kinematically associated with the latter gas component, rather than the gas pillars.

There is a gas pillar in subregion (b). This gas structure is moving at velocities of about 10 km s⁻¹ along the line of sight, which is smaller than the systemic velocity of NGC 2244 (12.8 km s⁻¹). It implies that this is part of the near side H ii bubble approaching observers. There is a small group of stars in the vicinity of the gas pillar. These stars have RVs similar to that of the gas pillar, although the number of stars with RV measurements is only three out of six. In addition, their PM vectors are systematically oriented toward the southeast (i.e., away from the central cluster).

The remaining gas in subregion (c) is part of the northwestern ridge of the Rosette Molecular Cloud. This molecular gas has RVs slightly larger than the systemic velocity of the central cluster, and therefore it is receding away from the cluster in the line of sight. There are more than 20 members of PL02 in the subregion. Their RVs range from 7 to 19 km s⁻¹. It is unclear whether or not PL02 is physically associated with this gas, due to the spread of the RVs. Given that the PM vectors of the PL02 members are small and random, the velocity field of this group may be closer to that of NGC 2244 than that of the molecular gas.

The formation of stellar clusters has been explained by two scenarios: monolithic formation (Kroupa et al. 2001; Banerjee & Kroupa 2013, 2015) and hierarchical mergers of subclusters (Bonnell et al. 2003; Bate 2009). Cluster rotation can also provide additional implications with regard to the cluster formation (Corsaro et al. 2017; Mapelli 2017). Here, we discuss the formation of the central cluster NGC 2244 in the two-theories framework.

Stellar clusters may form along the dense regions in a turbulent molecular cloud because star formation efficiencies can be high in this environment (Kruijssen 2012). Recent simulations of monolithic cluster formation have been capable of reproducing the internal structures and kinematics of clusters (Kroupa et al. 2001; Banerjee & Kroupa 2013, 2015; Gavagnin et al. 2017). Rapid gas expulsion and stellar feedback can significantly change the structure and kinematics of clusters (Tutukov 1978; Hills 1980; Lada et al. 1984; Kroupa et al. 2001; Banerjee & Kroupa 2013, 2015; Gavagnin et al. 2017). However, it is also possible to explain the formation of stellar clusters even without consideration of these effects. Our previous study (Lim et al. 2020) simulated the dynamical evolution of a model cluster formed in an extremely subvirial state. This model cluster collapses in the first several million years and expands. Subsequently, its structural features and kinematics were compared with those of the young open cluster IC 1805 (3.5 Myr—Sung et al. 2017) composed of a virialized core and an expanding halo. The model cluster had properties that compared well with the observed structure and kinematic properties of IC 1805.

Many other clusters in OB associations also tend to show expansion (Cantat-Gaudin et al. 2019; Kuhn et al. 2019), but it seems difficult to detect cluster rotation by using only PMs, as rotating clusters are rare in the sample of Kuhn et al. (2019). However, the rotation of NGC 2244 was found via RVs in this study. Therefore, the spectroscopic observations are crucial to perform a more complete search of cluster rotation. Indeed, the significant rotation of the young massive cluster R136 in the Large Magellanic Cloud was detected from spectroscopic observations (Hénault-Brunet et al. 2012). The rotational energy accounts for more than 20% of total kinetic energy of this cluster. The detection of cluster rotation implies that clusters can form in rotating molecular clouds or from hierarchical mergers of small stellar groups with different motions.

Systematic survey of molecular clouds in external galaxies found that molecular clouds are rotating (Rosolowsky et al. 2003; Tasker 2011; Braine et al. 2020). More than 50% of the observed molecular clouds have prograde rotation with respect to the galactic rotation. The angular momenta of these clouds may be acquired from the differential rotation of disk and self-gravity. The clouds with retrograde rotation can form via cloud–cloud collisions (Dobbs et al. 2011). The rotation of stellar clusters can thus be inherited from their natal clouds. However, it is not yet understood how efficiently the angular momenta of their natal clouds are redistributed to individual stars and clusters.

The hierarchical merging of gaseous and stellar clumps in a turbulent cloud can also lead to the formation of rotating clusters (Mapelli 2017). There are some possible candidates. The core of 30 Doradus in the Large Magellanic Cloud has an elongated shape, and it is composed of two stellar populations with different ages (Sabbi et al. 2012). The dense group is known to be the young massive cluster R136, which presents rotation (Hénault-Brunet et al. 2012), while a diffuse stellar group is located at the northeast of this cluster. It is believed that the elongated structure and age difference indicate a recent or an ongoing merger of these two groups (Sabbi et al. 2012). Similarly, the young massive cluster Westerlund 1 has an elliptical shape (Gennaro et al. 2011) and hosts two different age groups (Lim et al. 2013). However, the rotation of this cluster has not been reported yet.

NGC 2244 could monolithically form in a rotating cloud or filament hub, after which it might undergo cold collapse depending on its initial virial state. Some stars are now escaping from the cluster, as seen in Figure 8. Gas expulsion may also play a role in dispersing the cluster members. In addition, there are several stars in the north and south beyond r_cl, and these stars are also moving away from the cluster (see Figure 4). We traced their positions back using PMs. As a result, about 85% (11/13) of these stars were in the cluster (<r_cl) 2 Myr ago. Therefore, the spatial distribution of these stars can also be explained by the dynamical evolution of this cluster after the monolithic formation.

The observed features in IC 1805 could be explained by the monolithic collapse and rebound (Lim et al. 2020). In the same context, we compare the pattern of expansion in NGC 2244 with those in IC 1805. While IC 1805 shows a strong expansion pattern in its outer region (Lim et al. 2020), NGC 2244 has, in comparison, a rather weak tendency toward expansion (Figure 8). One possible explanation would be that the initial rotation supported the system against the subvirial collapse. Dynamical mass segregation is a natural consequence of the violent relaxation during the cold collapse (Allison et al. 2009). However, the absence of notable mass segregation in NGC 2244 also implies that the cluster did not have a strong cold collapse in the past.

The western half of the cluster has a negative net RV relative to the systemic velocity as it rotates (Figure 9), while NGC 2237 has a positive net RV (Figures 12 and 13). The discrepancy between the directions of the cluster rotation and stellar motions in the subregion cannot be simply explained by the monolithic cluster formation in a rotating cloud. On the other hand, star formation is still taking place along the filaments in the Rosette Molecular Cloud, forming several groups of stars (Poulton et al. 2008; Schneider et al. 2010; Cambrésy et al. 2013). In particular, a clustering of the stellar groups PL04, PL05, PL06, and REFL08 is noteworthy (Cambrésy et al. 2013; see also Cluster E in Poulton et al. 2008). This circumstance is a favorable condition for merging of stellar groups. Similarly, NGC 2244 has a substructure toward the southeast (Figure 4). The structural feature and rotation of this cluster may be the result of hierarchical mergers of stellar groups.

The Herschel observation reveals that the Rosette Molecular Cloud is being affected by the ultraviolet irradiation from the massive stars in NGC 2244 (Schneider et al. 2010). A gradient of dust temperature was found toward the southeast. The age differences among individual stellar groups indicate the presence of an age sequence among YSOs from the H ii bubble toward the Rosette Molecular Cloud. Schneider et al. (2010) argued that the formation of the stellar groups may be triggered by the radiative feedback from massive stars in NGC 2244.

However, this argument was later refuted by Cambrésy et al. (2013). They measured the slope of spectral energy distributions of individual YSOs and better estimated the ages of the stellar groups. The stellar groups close to the Rosette Nebula are systematically younger than NGC 2244. On the other hand, the stellar groups PL06 and PL07 are far from the Rosette Nebula, but their ages are comparable to that of NGC 2244. As a conclusion, their age distribution could not be fully understood in the context of the feedback-driven star formation (Elmegreen & Lada 1977). Star formation in the Rosette Molecular Cloud seems to take place in two different ways (Phelps & Lada 1997). Several previous studies have also reached similar conclusions regarding the formation of these stellar groups (Poulton et al. 2008; Román-Zúñiga. et al. 2008; Ybarra et al. 2013).

To understand the formation of the three groups that we defined at the border of the Rosette Nebula (the left panel of Figure 4), it is necessary to consider both models of hierarchical star formation along filaments (Bonnell et al. 2011) and feedback-driven star formation (Elmegreen & Lada 1977) as suggested by Phelps & Lada (1997). Figure 14 displays the distributions of T_ex and column density along the remaining molecular gas in these subregions. The heads of gas pillars in subregion (a) have T_ex and column densities lower than those in the other subregions. The column density traced by the optically thin ¹³CO line compared to ¹²CO line shows no agglomeration deeper inside the pillar heads.

Zoom In Zoom Out Reset image size

A large amount of molecular gas still remains in subregions (b) and (c). T_ex in subregions (b) and (c) appears higher than 25 K, which is consistent with dust temperature (Schneider et al. 2010). The molecular gas in subregion (b) condenses along the midplane, and a local condensation is found at the ridge of the molecular gas in subregion (c). These results indicate that molecular gas is compressed along the ridge of the Rosette Molecular Cloud because the clouds are exposed to the ultraviolet irradiation from the O-type stars in NGC 2244. Hence, subregions (b) and (c) might be the possible sites of feedback-driven star formation, compared to subregion (a). If PL02 is physically associated with the adjacent molecular gas, its age being younger than NGC 2244 (Cambrésy et al. 2013) may support this argument.

Here, we discuss the formation of the three groups based on their kinematics as shown in Section 4.2. If stars form in the remaining gas compressed by the expanding H ii bubble, it is expected that they are receding away from ionizing sources in NGC 2244 and have similar velocities to that of the gas. The stars in subregion (b) seem to follow this expectation. These stars are moving away from the cluster (toward the southeast) and have RVs similar to that of the adjacent gas pillar. In addition, this group contains a low-mass Class I YSO at an early evolutionary stage of protostars. This group may be slightly younger than the other groups, and hence its formation might be triggered by the feedback from massive stars in NGC 2244. We note that there is another Class I object in NGC 2244, but this may be a slightly old Herbig Ae/Be star candidate with a thick disk (Figure 6).

NGC 2237 in subregion (a) does not show any physical association between stars and the gas pillars (Figure 13). Also, this group is moving toward the south instead of the west (the direction of radial expansion of the H ii bubble). This confirms that region (a) is probably not a site of feedback-driven star formation. Although subregion (c) was considered as a possible site of feedback-driven star formation in Phelps & Lada (1997), the kinematic properties of PL02 and adjacent clouds do not support the claim. If PL02 was formed by the compression of the clouds, the members of PL02 should show PMs systematically receding away from NGC 2244. However, they have almost random motions, rather than any systematic motion toward the southeast as seen in subgroup (b). In addition, the physical association between PL02 and the remaining clouds is unclear. Hence, the formation of this group seems not to be related to feedback from the massive stars in the cluster. These two groups, NGC 2237 and PL02, might have been independently formed through the hierarchical fragmentation of a molecular cloud.