Looking Deep into the Rosette Nebula's Heart: The (Sub)stellar Content of the Massive Young Cluster NGC 2244

Star clusters are the primary locations of star formation and the main suppliers of stellar and substellar objects in the universe (Lada & Lada 2003; Portegies Zwart et al. 2010). The youngest clusters are particularly important laboratories for understanding how star formation works. They give birth to objects spanning at least four orders of magnitude in mass, from the most massive stars reaching several tens of solar masses, down to substellar objects below the deuterium-burning limit at ∼0.012 M_⊙. Richly populated young clusters therefore provide an ideal testbed in which to probe a variety of yet unsolved questions related to the formation of stars, brown dwarfs (BDs), and clusters.

The distribution of stellar masses in young clusters, or the initial mass function (IMF), has been an important subject of many observational and theoretical efforts. On the high-mass side (≳1 M_⊙), it can be approximated by a power-law form dN/dM ∝ M^−α, with α = 2.35 (Salpeter 1955). Below ∼0.5 M_⊙, the mass distribution is significantly flatter (see, e.g., Luhman 2012 and references therein). Although the Salpeter slope is often considered universal (see, e.g., Bastian et al. 2010; Offner et al. 2014), there are some notable examples of young clusters where the mass function appears to be flatter, with α = 1.3–1.8 (Stolte et al. 2006; Harayama et al. 2008; Bayo et al. 2011; Peña Ramírez et al. 2012; Habibi et al. 2013; Andersen et al. 2017; Mužić et al. 2017). While the flattening of the slope could be an effect of crowding for the clusters that are several kiloparsecs away (NGC 3603, Arches, Westerlund 1) (Ascenso et al. 2009), this is certainly not the case for the more nearby ones (Collinder 69, σ Ori, RCW 38). Employing Bayesian statistics, Dib (2014) reports significant variations in IMF parameters of eight young Galactic clusters, including the slope of the IMF on the high-mass side. Recently, Dib et al. (2017) proposed a probabilistic formulation of the IMF, where the IMF is described by a tapered power law (De Marchi et al. 2010), and each parameter of this function is represented by a Gaussian probability distribution. The authors argue that the broad distribution of the IMF parameters might reflect the existence of equally broad distributions in the star-forming conditions (e.g., Svoboda et al. 2016).

On the low-mass side, the observed IMFs in young clusters and star-forming regions in general yield slopes α in the range 0.5–1, for masses below ∼1 M_⊙ and extending into the substellar regime (e.g., Luhman 2004a, 2007; Bayo et al. 2011; Peña Ramírez et al. 2012; Scholz et al. 2012; Lodieu 2013; Mužić et al. 2015, 2017). The number ratio between low-mass stars and BDs is typically found in the range 2–5. The range of values reflects different sources of uncertainties that contribute to the derivation of the IMF slope and the star-to-BD ratio (uncertainties in distances, ages, use of different evolutionary models, reddening laws, membership issues, small sample sizes; see, e.g., Scholz et al. 2013), rather than the variations between individual clusters. Although from the observations there is no convincing evidence for strong variations in the BD production efficiency in nearby star-forming regions, theoretical expectations are somewhat different. Most of the current BD formation theories in fact predict an increased production of substellar objects in dense environments or close to very massive stars. The factors that facilitate BD formation are high gas densities (Padoan & Nordlund 2002; Bonnell et al. 2008; Hennebelle & Chabrier 2009; Jones & Bate 2018), high stellar densities that favor BD formation through ejections (Bate 2012), and the presence of massive stars, where BDs can be formed through photoevaporation (Whitworth & Zinnecker 2004) or by fragmentation in their massive disks (Stamatellos et al. 2011; Vorobyov 2013). Of the mentioned scenarios, a few give quantitative predictions of the impact of environment on the efficiency of BD production. In the scenario by Bonnell et al. (2008), where new objects form by gravitational fragmentation of gas infalling into a stellar cluster, BDs are preferentially formed in regions with high stellar density. An increase in object density by an order of magnitude would result in an increase in the BD fraction by a factor of about two. The predicted effect is more subtle in the hydrodynamical simulations of molecular cloud collapse by Jones & Bate (2018). An increase in gas density by a factor of 100 leads to BD frequencies larger by ∼35%.

To explore potential environmental influences on BD formation, we have initiated a photometric study in a sample of young, massive clusters, characterized by high stellar densities and/or a substantial population of massive stars. In this sense, the environments we aim to cover are significantly different from those in nearby star-forming regions, with the exception of the Orion Nebula Cluster (ONC), the only site of massive star formation within 500 pc from the Sun. The latest results from the ONC show a bimodal form of an IMF, with a secondary peak in the substellar regime (Drass et al. 2016). The authors interpret the substellar peak as possibly formed by BDs ejected from multiple systems or circumstellar disks at an early stage of star formation. The spectroscopic confirmation of this result, however, is still pending. In the first paper resulting from our study of massive clusters, we studied the core of RCW 38, a young (1 Myr), embedded cluster located at a distance of 1.7 kpc (Mužić et al. 2017). RCW 38 is characterized both by high stellar densities (core stellar surface density Σ ∼ 2500 pc⁻², i.e., twice as dense as the ONC and more than ten times denser than NGC 1333; Mužić et al. 2017) and by a substantial population of massive OB stars (60 OB candidates in total, and a dozen in our small field of view of 0.5 × 0.5 pc²; Wolk et al. 2008; Winston et al. 2011). Given these characteristics, and the absence of more extreme environments that are close enough to allow the detection of BD candidates, this cluster was thought to be a good starting point in which to look for possible environmental differences in the substellar regime. Unlike what has been reported for the ONC, we found that RCW 38's IMF is consistent with those in nearby star-forming regions (α ∼ 0.7 for the range 0.02–0.5 M_⊙).

In this second paper of the series, we present the observation of the young cluster NGC 2244, located at the center of Rosette Nebula, one of the most prominent features of the Mon OB2 Cloud. The cluster contains more than 70 massive OB stars, which are presumed to be responsible for the evacuation of the central part of the nebula (Román-Zúñiga & Lada 2008). Estimates of the distance to NGC 2244 vary from 1.4 to 1.7 kpc (Ogura & Ishida 1981; Perez et al. 1987; Hensberge et al. 2000; Park & Sung 2002; Lombardi et al. 2011; Martins et al. 2012; Bell et al. 2013; Kharchenko et al. 2013), whereas most authors agree on ∼2 Myr for the age of the cluster (Perez et al. 1987; Hensberge et al. 2000; Park & Sung 2002; Bell et al. 2013).

The stellar content of NGC 2244 has been extensively studied from X-ray to mid-infrared (MIR) wavelengths. The most sensitive X-ray data to date are the Chandra observations published in Wang et al. (2008). The authors identify more than 900 X-ray sources in a ∼17' × 17' field, of which 77% have optical or near-infrared (NIR) counterparts and are considered young cluster members. The X-ray-selected population of Wang et al. (2008) is nearly complete down to 0.5 M_⊙. The deepest published NIR data on NGC 2244 are those taken with Flamingos at the 2.1 m telescope at the Kitt Peak National Observatory, as part of the imaging survey of the Rosette complex (complete down to J ∼ 18.5 mag; Román-Zúñiga et al. 2008). The analysis of several clusters in the complex reveals NGC 2244 to be the largest one. That paper, however, deals with the entire Rosette complex, and does not present further specific analysis of the NGC 2244 cluster population. The early attempts to derive the IMF of the cluster reported a shallow slope α ∼ 1.7−1.8 (Massey et al. 1995; Park & Sung 2002) for stars with m > 3 M_⊙, but at the same time Park & Sung (2002) cautioned about the incompleteness at the intermediate and low masses in their sample. More recently, Wang et al. (2008) derive the slope α ∼ 2.1 for the mass range 0.5–30 M_⊙, close to the Salpeter slope.

The Rosette Nebula has also been studied within the MYStIX project (e.g., Kuhn et al. 2014, 2015), who identify 255 probable members in the regions Rosette D and E, overlapping with the position of the NGC 2244 cluster. They estimate the peak surface density of Σ ≈ 300 stars pc⁻² in the region Rosette D, centered close to the star HD 46150, located also in the center of our observed field. In the circular area with a radius of 0.5 pc around the peak position, the average surface density is ∼100 stars pc⁻², assuming a distance of 1.4 kpc. This value is similar to those found in nearby star-forming regions.

This paper is structured as follows. Section 2 contains the details of the observations and data reduction. The data analysis, including the point-spread function (PSF) fitting, photometric calibration, and completeness analysis, is presented in Section 3. In Section 3.2, we derive the distance to the cluster. In Section 4, we discuss the cluster membership, derive the masses and their distribution, estimate BD frequencies for NGC 2244, and discuss the stellar densities in various clusters and star-forming regions. The results are discussed in Section 5. Summary and conclusions are given in Section 6.

The photometry was performed using the DaophotII/Allstar PSF-fitting algorithm (Stetson 1987), with a Gaussian PSF that varies quadratically with position in the frame (VARIABLE PSF parameter set to 2). We excluded all the sources with the Daophot sharpness parameter $| {sh}| \gt 0.7$, which helps to get rid of galaxies, clumps, knots, and spurious detections caused by ghosts in our images. The photometric zero-points and color terms were calculated from a comparison with the data from the United Kingdom InfraRed Telescope Infrared Deep Sky Survey (UKIDSS; Lawrence et al. 2007) Galactic Plane Survey (Lucas et al. 2008) Data Release 10. This catalog was chosen because of its photometric depth, which allows a significant overlap with the F2 source lists. The photometric uncertainties provided by the UKIDSS pipeline only include random photon noise error, and are therefore underestimated. We use a relation provided by Hodgkin et al. (2009) to calculate more realistic photometric errors that also include the systematic calibration error component. Stars with counts surpassing the linear regime of the F2 detector were excluded from the comparison. The transition to the nonlinear regime occurs at J ∼ 15 mag, H ∼ 16 mag, and K ∼ 15 mag. To avoid potential issues due different spatial resolutions of the F2 and UKIDSS data, we have discarded those sources with a neighbor at a distance <10 and brightness contrast ΔH < 2 mag. Furthermore, we have rejected objects with uncertainties larger than 0.1 mag in the UKIDSS data or in the F2 instrumental magnitudes. Finally, the photometric zero-points and color terms were calculated using the following equations:

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}J & = & {J}_{{\rm{instr}}}+{{\rm{ZP}}}_{1}+{c}_{1}\ast (J-K)\\ H & = & {H}_{{\rm{instr}}}+{{\rm{ZP}}}_{2}+{c}_{2}\ast (J-K)\\ K & = & {{Ks}}_{{\rm{instr}}}+{{\rm{ZP}}}_{3}+{c}_{3}\ast (J-K).\end{array}\end{array}\end{eqnarray} \tag{ 1 }$$

For the cluster field, the zero-points and the color terms are ZP₁ = 25.07 ± 0.02, ZP₂ = 25.25 ± 0.02, ZP₃ = 24.68 ± 0.03, c₁ = 0.06 ± 0.02, c₂ = 0.04 ± 0.02, and c₃ = −0.08 ± 0.03. For the control field we obtain ZP₁ = 25.00 ± 0.04, ZP₂ = 25.28 ± 0.04, ZP₃ = 24.64 ± 0.05, c₁ = 0.05 ± 0.04, c₂ = −0.01 ± 0.04, and c₃ = −0.06 ± 0.05. The photometric uncertainties were calculated by combining the uncertainties of the zero-points and color terms, and the measurement uncertainties supplied by Daophot, and are shown in Figure 2 for the cluster field. The median uncertainties of the F2 photometry are 0.04 mag in J- and H-bands and 0.05 mag in the K-band. The photometry of the nonlinear stars was replaced by the UKIDSS measurements, with the exception of the two brightest stars in the field which are saturated even in UKIDSS. Their photometry was taken from the Two Micron All Sky Survey (2MASS, Skrutskie et al. 2006), and transformed to the UKIDSS photometric system using the relations given in Hodgkin et al. (2009). Finally, we removed all the sources identified as galaxies in the UKIDSS catalog, verifying visually that they are not resolved into multiple sources in F2 images (eight galaxies removed). The final JHK catalog contains 793 objects (Table 2).

Zoom In Zoom Out Reset image size

The completeness of the photometry was assessed with an artificial star test. Artificial stars were added to the original image using Daophot's function ADDSTAR. The resulting image was then used as the input for a routine identical to the one used to obtain the photometry for the cluster and the control field (Section 3.1), using the same PSF solution. We inserted only 50 stars at a time to avoid crowding, and repeated the procedure 40 times at each magnitude (in steps of 0.2 mag) to improve statistics. The ratio between the number of recovered and inserted artificial stars is plotted against magnitude in Figure 3 for both the cluster (black lines) and the control field (blue lines). The 90% completeness limits are J = 21.7 mag, H = 20.9 mag, and K = 20.2 mag for NGC 2244, and J = 21.7 mag, H = 21.1 mag, and K = 20.1 mag for the control field.

Zoom In Zoom Out Reset image size

The distance to NGC 2244 was previously estimated to be between 1.4 and 1.7 kpc (Ogura & Ishida 1981; Perez et al. 1987; Hensberge et al. 2000; Dias et al. 2002; Park & Sung 2002; Lombardi et al. 2011; Martins et al. 2012; Bell et al. 2013; Kharchenko et al. 2013). The recently published Gaia Data Release 2 (DR2, Gaia Collaboration et al. 2018) provides precise positions, proper motions, parallaxes, and G-band magnitudes for more than 1.3 billion sources. We queried the Gaia DR2 catalog within 10' radius from the brightest star in the field (HD 46150; 06:31:55.52, +04:56:34.3). This area contains the region observed with F2, but covers a much larger portion of the cluster. The resulting list contains 3498 objects. We then cross-matched the Gaia DR2 result with the catalog of 718 clusters members from Meng et al. (2017), constructed based on X-ray detection, infrared excess, proper motion, optical color selection, and a cut-off in the clustrocentric distance. The total number of matched sources with a valid five-point astrometric solution is 528.

The proper motions of the Gaia sources in the field are shown in the top panels of Figure 4. There is a distinct overabundance of sources close to the position μ_α, μ_δ ≈ (−1.5, 0.2) mas yr⁻¹, corresponding to the location of the cluster members. An enlarged version of this plot is shown in the upper right panel, with the blue open circles marking the matched members from Meng et al. (2017). The mean proper motion of the matched members is μ_α = −1.62 ± 0.65 mas yr⁻¹ and μ_δ = 0.15 ± 0.70 mas yr⁻¹.

Zoom In Zoom Out Reset image size

The color–magnitude diagram (CMD) of the same set of sources constructed by matching the Gaia DR2 to the UKIDSS survey shows the probable member sequence clearly standing out to the right of the bulk of the field stars. We can therefore use the proper motions and the position in the CMD as criteria to assign cluster membership. To consider a source a member in this exercise, we require it to be located inside the red ellipse in the proper motion space, to the right of or above the dividing line in the CMD, and to be classified as a star in UKIDSS (class = −1). The semimajor and semiminor axes of the red ellipse are equal to 1σ widths of the proper motion distributions in R.A. and decl. The transition from the pre-main sequence to the main sequence at the age of NGC 2244 occurs roughly at J = 12 (see Figures 16 and 8), motivating the horizontal red line. After applying the cuts mentioned above, we keep 345 out of the initial 3498 sources. These 345 sources are shown with red open circles in the parallax plot in the lower right panel of Figure 4. 242 of the 345 selected sources belong also to the members list of Meng et al. (2017). We note that our criteria may cause us to miss some members (mainly due to the relatively narrow allowed proper motion space, especially considering individual proper motion uncertainties), but the idea here is to minimize the number of interlopers that could bias our distance analysis, rather than to do a sophisticated membership analysis at this point.

We estimate the distance to the cluster through a maximum likelihood procedure (see, e.g., Cantat-Gaudin et al. 2018), maximizing the following quantity:

$$\begin{eqnarray}&&{ \mathcal L }\propto \displaystyle \prod _{i}P({\varpi }_{i}| d,{\sigma }_{{\varpi }_{i}})=\displaystyle \prod _{i}\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{{\varpi }_{i}}^{2}}}\exp \left(-\displaystyle \frac{{\left({\varpi }_{i}-\tfrac{1}{d}\right)}^{2}}{2{\sigma }_{{\varpi }_{i}^{2}}}\right),\end{eqnarray} \tag{ 2 }$$

where $P({\varpi }_{i}| d,{\sigma }_{{\varpi }_{i}})$ is the probability of measuring a value ϖ_i for the parallax of star i, if its true distance is d and its measurement uncertainty is ${\sigma }_{{\varpi }_{i}}$. Here we assume that there are no correlations between individual parallax measurements, and that the likelihood of the cluster distance can be represented as the product of the individual likelihoods of all its members. We neglect any cluster spread, i.e., all the stars are assumed to be located at the same distance, which should have a negligible effect on our distance determination, given the relatively large distance to the cluster. In this way we obtain a distance to NGC 2244 of d = 1.59 ± 0.02 kpc. If we restrict ourselves only to the members identified by Meng et al. (2017) that also pass our cuts (242 objects), we obtain a very similar value of d = 1.57 ± 0.02 kpc. Our distance estimate is in agreement with the most recent value ${1.55}_{-0.09}^{+0.1}$ kpc from Kuhn et al. (2019), derived also from the Gaia DR2 data and assuming the systematic parallax uncertainty of 0.04 mas.

It has been shown that Gaia DR2 parallaxes suffer systematic uncertainties, with a global zero-point offset of −0.029 mas (Arenou et al. 2018; Lindegren et al. 2018; Luri et al. 2018). However, applying this correction to small areas is not advisable, since the local variations can be significantly different (Arenou et al. 2018). Moreover, a low number of observed QSOs prevents a calculation of the local bias close to the Galactic plane. Arenou et al. (2018) report a median systematic offset in parallax of ∼−0.065 mas for the cluster samples of Kharchenko et al. (2013) and Dias et al. (2002), with a dependence on cluster distance and stellar colors. Judging from the right panels of Figure 16 of Arenou et al. (2018), an appropriate offset value for NGC 2244 is ∼−0.08 mas, i.e., the Gaia distance derived here might be overestimated by ∼180 pc, or 11%. Considering this potential systematics, it might seem that the distance constrained by Gaia does not provide a significant improvement with respect to previous measurement; however, we have to keep in mind that this is only the second data release, and that by the end of the mission, the systematics are expected to be at a microarcsecond level even at distances of a few kiloparsecs. To understand the effect this uncertainty might have on the results, we will report them for two values of the distance, 1.6 and 1.4 kpc.

As a side note, in the bottom left panel of Figure 4, there is a group of relatively bright sources not recognized as members by Meng et al. (2017), located to the right of and above the red lines (roughly J < 13 and G − J > 1). Only about 10% of these sources pass all our membership cuts. The parallaxes of those that do not, signal that they are mostly located behind the cluster, and their position in absolute CMDs is consistent with a giant branch. About a quarter of these objects have been observed by the LAMOST survey,¹⁰ and have determinations of surface gravity (log g) and effective temperature (T_eff). About 75% of those have T_eff = 4500–5000 K and log g < 3, i.e., they are background giants.

Since this is a long and relatively complicated section, here we present a short summary of the methodology, giving readers a chance to skip some sections that might be too technical or detailed.

• 1.  
  The goal of this section is to derive the IMF and the star-to-BD number ratio for NGC 2244 (Sections 4.4 and 4.5). The results are presented for three different isochrones and two distances (1400 and 1600 pc), in order to estimate the effect that different assumptions have on the results.
• 2.  
  Membership determination (Section 4.2). The first step in the analysis is to separate cluster stars from the contaminants sharing the same line of sight. For the sources brighter than J ∼ 18, this can be done using the color, proper motion, and parallax cuts on the Gaia DR2 and Pan-STARRS data (Section 4.2.1). However, our NIR catalog extends more than 3 mag deeper, and in this regime we employ a statistical determination of membership, comparing the CMD of the population observed in direction of the cluster to that of the control field (Section 4.2.2). In this procedure, there are several parameters that are varied (cell sizes and positions, difference in extinction between the cluster and the control field), and each combination of these parameters can give a slightly different solution. We keep all the possible solutions (324 member lists) and use them later to evaluate the uncertainties that this approach introduces to the IMF and star-to-BD ratio.
• 3.  
  Isochrones (Appendix A). Stellar parameters are derived in comparison to theoretical models. The basic products of the stellar models are the bolometric luminosity and effective temperature, which are converted into magnitudes and colors, by applying bolometric corrections and T_eff–color relations. The conversion can be done by using the atmosphere models or empirical relations. Here we decide to test both approaches, and we derive our results using three different isochrones.
• 4.  
  Derivation of stellar parameters (Section 4.3). We compare three different methods to derive effective temperature, mass, and extinctions. The first method is based on fitting the spectral energy distribution (SED) over a large wavelength range (optical to mid-infrared, where available; Section 4.3.1). The other two methods explore the derivation of parameters from the three NIR colors only (JHK), in one case assuming that the observed colors are photospheric (Section 4.3.2), and in the other case assuming that some intrinsic excess at the NIR wavelengths is present (Section 4.3.3). We find that the results are consistent with each other, validating the use of only NIR magnitudes and colors in mass determination at the faint end, where this is our only option (Section 4.3.4).

4.2.1. Gaia DR2 and Pan-STARRS1

For the brighter portion of sources in our field of view, data from Gaia DR2 and Pan-STARRS1 (PS1; Chambers et al. 2016) are available. 199 sources from our JHK catalog have a counterpart with a five-parameter astrometric solution in Gaia DR2 (down to J ≈ 17.5−18) within the matching radius of 05. The average distance between the matched entries in the two catalogs is 003, justifying our choice of the relatively small matching radius. We consider a source to be member candidate if it satisfies the following criteria, as depicted in Figure 5: (1) it is located above and to the right of the red lines in the J, (G − J) diagram (open gray circles in top panels); (2) its 1σ proper motion errors intersect the circular area with a radius of 1.5 mas yr⁻¹, centered at the average value of the proper motions derived in Section 3.2 (red circle in the top right panel); (3) its 1σ errors in parallax intersect one of the red lines plotted in the bottom left panel (marking 1.4 and 1.7 kpc distance, see Section 3.2). In total 72 of 199 sources pass our cuts and are considered member candidates. They are marked with orange diamonds in the bottom right panel of Figure 5.

Zoom In Zoom Out Reset image size

PS1 contains neither proper motions nor parallaxes, but its optical photometry is deeper than Gaia, and can further help in selection of candidate members. As shown in the left and middle panels of Figure 6, the sources are separated in two bands when combining optical with NIR photometry.¹¹ Orange diamonds mark the member candidates selected by Gaia, and gray crosses those rejected in the same procedure. Following the sequence of the Gaia member candidates, we draw dividing lines in r, (r − J) and i, (i − J) CMDs (solid red lines). We select the sources not classified by Gaia and located to the right of the red lines in both diagrams as member candidates (blue open circles). 47 out of 362 F2–PS1 matches (matching radius 05) with valid r or i magnitudes are selected in this way. Judging from the Gaia color versus astrometric selection, ∼15 of these sources might still be interlopers. The selected sources are identified in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  Member Candidates from Gaia DR2 and PS1, along with the Best-fit Parameters from the SED

ID	α (J2000)	δ (J2000)	Teff (K)	AV (mag)	log(g (cm s−2))	Selection
2	06:31:50.58	04:57:29.3	3100	1.9	3.75	PS1
3	06:31:50.64	04:58:35.8	6000	4.5	4.25	Gaia DR2
5	06:31:50.85	04:57:38.0	3100	2.9	3.75	PS1
6	06:31:50.86	04:57:20.0	2900	0.2	3.50	PS1
18	06:31:52.28	04:58:33.9	3000	1.2	3.50	PS1

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The selection using Gaia and PS1 is nearly complete down to J ∼ 18.5. About 6% of the F2 sources with J < 18.5 (21/337) cannot be classified as members or non-members in this way, because no counterpart with valid r or i magnitude is found in PS1, or no astrometry is given in Gaia. Most of these sources are located very close to the two brightest stars in our FoV, which might have affected their measurements in Gaia or PS1.

We matched the sources selected in this section with the member list of Meng et al. (2017), who select as member candidates all sources detected in both X-rays and Spitzer/IRAC, or those that show excess at wavelengths of 8–24 μm. The depth of their survey seems to be fairly complete down to the 2MASS completeness limit at J = 15.8 mag, and we restrict the comparison down to that value. Of the 61 Gaia/PS1 selected candidates brighter than J = 15.8, 47 are found in the member list of Meng et al. Of the 14 sources that are not in their list, 10 are not in the X-ray/NIR combined catalog of Wang et al. (2008), which was used for selection in the former work. The remaining four objects were observed by Spitzer (Balog et al. 2007), but they do not show an excess at these wavelengths, which is possibly the reason why they were rejected by Meng et al. (2017). On the other hand, of the 59 objects rejected by Gaia/PS1, 20 are found in the member list of Meng et al. Of these 20, six were rejected based on either colors or proper motions, and the remainder based on parallax. It might be worth noting that the relative error of the parallax for these sources is relatively large (median ∼0.3), but only for four of them do the lower and upper bounds of the distance confidence interval from Gaia DR2 (Bailer-Jones et al. 2018) agree with the distance to NGC 2244.

4.2.2. Flamingos-2

Our F2 JHK photometry extends more than three magnitudes in J deeper from the objects detected in PS1. In the magnitude range covered by Gaia/PS1 the probable member sequence separates nicely from the background even in the NIR CMD, but below J ∼ 17 the cluster members occupy the same space in the CMD as the foreground and background field stars. To estimate the relative contribution of the contaminants to the stellar sample observed in the direction of the cluster, we compare the cluster CMD with that of the control field (see Figure 7). The method is similar to the one we applied to the cluster RCW 38 (Mužić et al. 2017) but includes a few modifications. The main steps of the method are as follows.

• 1.  
  The CMD is subdivided into grid cells with a step size (Δcol, Δmag) in the color and magnitude axes, respectively.
• 2.  
  For each CMD cell, we calculate the expected number densities of stars in the cluster field and the control field. The number density of field stars has in principle to be scaled to account for different on-sky areas covered by the images, although in our case this number is close to one.
• 3.  
  A number of objects corresponding to the expected field object population is then randomly removed from corresponding cells of the cluster CMD.

Zoom In Zoom Out Reset image size

To calculate the number densities of stars in each cell, we use the method outlined in Bonatto & Bica (2007). Each star is represented by a Gaussian probability distribution of magnitude and color, with the widths determined by the uncertainties. The expected number density of stars in each cell is calculated as the sum of each individual star's probability of being in a corresponding cell. In this way we take into account the fact that stars with large uncertainties have a non-negligible probability of populating more than one CMD cell.

The choice of the cell size is an important parameter to take into account, and results from a compromise between having a sufficient number of stars in each cell and preserving the morphology of different CMD evolutionary sequences. We have tested the cell size of Δcol = (0.4, 0.5, 0.6) mag and Δmag = (0.4, 0.6, 0.8) mag. Furthermore, for each combination of the color and magnitude cell sizes, we repeat the procedure with the grid shifted by ±1/3 of the cell width in each dimension. This results in 81 different configurations.

The control field population seems to have, on average, a slightly higher extinction than the cluster field. This can be appreciated by looking at the upper sequence in the CMD, where there is a slight shift between the sequence of the cluster on the left and the control field. We estimate that the difference is somewhere in the range 0.5–1 mag. The control field sequence is therefore dereddened at the beginning of the procedure by ΔA_V = (0.5, 1.0) mag. Finally, we applied the procedure to both the J, (J − H) and J, (J − K) CMDs. In total, this procedure results in 81 × 2 × 2 = 324 different member lists, which can then be used to calculate statistical errors of the IMF and the star-to-BD number ratio, which are introduced by the decontamination procedure.

We estimate that (61 ± 2)% of the sources observed toward the cluster are actually field contaminants. The number of cluster stars surviving the procedure is 310 ± 18 (the uncertainty is the standard deviation of the 324 outcomes of the above procedure).

Our procedure successfully recovers the member sequence at J ≲ 17 mag, which appears separated from the field objects located to the left (see Figure 8). Green diamonds in this plot show the candidate members from the Gaia and PS1 selection.

Zoom In Zoom Out Reset image size

In this section we explore different approaches to derive stellar effective temperatures, masses, and extinction (A_V) of candidate members of NGC 2244 from multi-band photometry. For the candidate members selected with the help of Gaia and Pan-STARRS1, we can take advantage of the existence of photometry over a large wavelength range, and derive these parameters by fitting their SEDs to the synthetic photometry obtained from atmosphere models. This, however, only works for the bright part of our catalog. For the objects with J ≳ 18.5 where we only have JHK photometry available, the two-parameter SED fitting becomes degenerate. Instead, to simultaneously derive the extinction and the masses of the cluster members, we compare the JHK photometry to model isochrones. The first approach here is to assume that the observed colors come from the reddened stellar photosphere. However, as argued by Cieza et al. (2005), classical T-Tauri stars (CTTSs) can present an excess already in the J-band, and therefore a simple dereddening of the photometry to the model isochrones can overestimate the extinction, and consequently also the stellar luminosity and mass. About 70% of the stellar members of NGC 2244 host dusty disks or envelopes (Meng et al. 2017), which can contribute some excess flux to their near-infrared colors. For this reason we also develop a method that takes into account this possible intrinsic excess when deriving the extinction toward individual stars in the cluster.

Throughout this section, we assume an age of 2 Myr (when employing the isochrones; Perez et al. 1987; Hensberge et al. 2000; Park & Sung 2002; Bell et al. 2013), and adopt the extinction law from Cardelli et al. (1989) with the standard value of the total-to-selective extinction ratio R_V = 3.1 (Fernandes et al. 2012). For a detailed description of the isochrones used in Sections 4.3.2 and 4.3.3, see Appendix A. We derive the parameters assuming the cluster distance of 1.6 kpc derived from the Gaia data, and use this for the comparison of the different methods. At the end, we also derive the masses for the distance of 1.4 kpc, which should be a lower limit on the distance judged from the Gaia systematics.

4.3.1. Optical to Near-infrared SED Fitting

We construct the SEDs of the Gaia/PS1 candidate members using the available optical and NIR photometry mentioned in the previous sections. We also use the Spitzer/IRAC photometry from Balog et al. (2007) where available, but only for the sources not showing excess at mid-infrared wavelengths (15 sources). Whether a source exhibits an excess was determined with the help of the tool VOSA (Bayo et al. 2008), which detects mid-infrared excess based on the slope of the infrared SED. The validity of this selection was also visually inspected. We retrieved the BT-Settl atmosphere models (Baraffe et al. 2015) for log g values of between 3.5 and 4.5, and T_eff between 2000 and 25,000 K¹² and interpolate the models to the intermediate log g values of 3.75 and 4.25 (young stars with an age of up to a few million years and masses below 10 M_⊙ should have log g between 3.5 and 4.25). The metallicity is set to the solar value. The extinction was varied between 0 and 10 in steps of 0.1 mag. We find the model best matching the data by minimizing the χ² parameter calculated as

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{\left({F}_{i\_\mathrm{obs}}-D\times {F}_{i\_\mathrm{model}}\right)}^{2}}{{\sigma }_{i\_\mathrm{obs}}^{2}},\end{eqnarray} \tag{ 3 }$$

where F_obs and F_model are observed and model fluxes, respectively, and σ_obs are observed flux uncertainties. D is the dilution factor, D = (R/d)², with R being the stellar radius and d the cluster distance (fixed to 1600 pc). The value of R is fixed for each combination of T_eff and log g and was extracted from the models. This is motivated by the relation $R=\sqrt{G\times m/g}$, where m stands for the object's mass, which is directly related to T_eff at a particular value of log g. We visually inspected all the fits, and rejected two of them due to a clear discrepancy between the observed and the best-fit model fluxes. In Figure 9 we show a subset of the object SEDs, together with the best-fit model synthetic photometry (red). In the leftmost panels we show two objects with optical to near-infrared photometry, in the middle panels those with no excess in Spitzer photometry, and in the rightmost panels those showing the excess (the excess points were not used to fit the SED). The best-fit parameters are given in Table 3.

Zoom In Zoom Out Reset image size

Finally, the masses were obtained by interpolating the T_eff–mass relation from the 2 Myr models to the best-fit T_eff values.

4.3.2. From NIR Colors, without Excess

4.3.3. From NIR Colors, with Excess

4.3.4. Comparison of the Three Methods

In Figure 10, we compare the three mass/T_eff methods described in Sections 4.3.1 (SED fitting; method 1 for the rest of the section), 4.3.2 (JHK without excess; method 2), and 4.3.3 (JHK with excess; method 3). The top panels compare the results of method 3 with those of method 1, and the bottom panels method 3 with method 2. In the left and middle panels, we show the comparison of T_eff and masses, where in the case of methods 2 and 3 we show only the results using isochrone 1. The T_eff error bars for the SED fitting are set to 100 K, reflecting the spacing of the fitting grid. For the methods based on JHK photometry only (methods 2 and 3), we have a distribution of T_eff and masses for each object, which does not necessarily follow the normal distribution, and in most cases is not symmetric. The T_eff and mass for the plot were calculated as the weighted average, with the weights provided by the distribution function. The error bars represent 95% confidence limits. We note that the large T_eff error bars for the massive stars in methods 2 and 3 arise from the T_eff/A_V degeneracy at the main-sequence transition region of the isochrones. Basically, there is more than one combination of T_eff and A_V that can be used to deredden a star's photometry onto the isochrone, and the span of T_eff values in this region is large (∼14,000 K).

Zoom In Zoom Out Reset image size

The right panels of Figure 10 show the IMFs derived using the masses obtained by different methods. The IMFs are derived in the same way as explained in the following section. For the upper plot the input object list contains the Gaia/PS1 candidate members, while in the lower plot the IMF is averaged over 324 lists resulting from the control-field-based membership analysis (Section 4.2.2).

• 1.  
  Method 3 versus method 1  (top panels of Figure 10): The values of T_eff obtained by the two methods are in reasonable agreement below ∼6000 K, where most of the objects are located. Above this temperature there is a set of objects with significantly higher T_eff from the JHK method (albeit with very large error bars). These are the objects that occupy the space in the CMD around the pre-main-sequence transition (masses roughly in the 2–6 M_⊙ range, see Figure 8). Their T_eff and mass distribution tends to be bimodal, resulting in the larger error bars, and values typically higher than those obtained by the SED fitting. There are also a handful of sources just above the smaller blue box in the left panel of Figure 10. SED fitting for these sources prefers higher T_eff and A_V than the JHK method. There is no obvious explanation for this discrepancy, but we note that most of them (4/5) are either at the faint end in Gaia and have very large parallax errors, or were selected in PS1 (color selection only). They might simply be contaminants. For the objects located inside the blue box of the top left plot of Figure 10, we obtain the average T_eff difference (method 3 – method 1) and its standard deviation ΔT_eff = −80 ± 310 K. The resulting IMFs are very similar. The difference in the bin at around 3 M_⊙ comes from the five objects discussed above, where SED fitting yields higher masses. Removing those objects, the red point becomes consistent with the black one. To conclude, the overall good agreement between the parameters obtained by the two methods, and in particular by the similar IMFs, is reassuring, and validates the methods based on JHK colors only. This is important, because it represents the only option to derive T_eff and mass for the fainter sources.
• 2.  
  Method 3 versus method 2  (bottom panels of Figure 10): As expected, the method ignoring the NIR excess predicts systematically higher T_eff and masses than the method allowing for the infrared excess. The average T_eff difference (method 3 – method 2) and its standard deviation ΔT_eff = −80 ± 150 K. The results agree well within the errors, and the impact on the resulting IMF is minimal, and only at the lowest-mass bins. Clearly, ignoring the NIR excesses in the determination of mass does not seems to be a crucial point for studying the IMF, but anyway for the remainder of the paper we decide to use the masses derived from method 3.

In Figure 11 we show histograms of the extinction values derived by the three methods. The top panel is comparing the A_V values from methods 2 and 3 for the objects from the full JHK sample. Note that the A_V derivation relies on dereddening the photometry to the 2 Myr isochrone shifted to the distance of the cluster. The middle and bottom panels show the histograms of the three methods, but only for the objects in the Gaia/PS1 sample of candidate members. The average extinction of the cluster is 1.7 ± 1.1 mag from the SED fitting (method 1), 1.3 ± 1.1 mag from JHK photometry without excess (method 2), and 1.4 ± 0.9 mag from JHK photometry with excess (method 3). The listed uncertainties are standard deviations of the three distributions. The extinctions derived by the three methods are in general consistent with each other, and we see no large intracluster spatially variable spread of extinctions.

Zoom In Zoom Out Reset image size

In this section we derive the IMF of NGC 2244, for masses in the range from 0.02 to 7 M_⊙. Stellar masses are binned with logarithmic bin sizes of Δlog(m/M_⊙) = 0.2. As can be appreciated from Figure 8 our lower mass limit is well above the 90% completeness limit at the average cluster extinction (A_V = 1.5 mag), and should also contain the complete dereddened control field population to ensure the lowest-mass bins are not artificially overpopulated.

In Section 4.2 we have derived 324 lists containing statistically cleaned probable members. Each object was assigned a mass distribution using the method described in Section 4.3.3. For each of the 324 member lists, we run a Monte Carlo simulation to derive the values and uncertainties of dN/dM for each bin. The dN/dM values of the resulting 324 IMF versions are combined (as a weighted average) to derive the final cluster IMF. The Monte Carlo simulation is very similar to the one previously applied in the cluster RCW 38 (Mužić et al. 2017). The mass of each star is drawn from the distribution derived in Section 4.3.3. This is performed 100 times, and for each of the 100 realizations we do 100 bootstraps, i.e., random samplings with replacement. In other words, starting from a sample with N members, in each bootstrap we draw a new sample of N members, allowing some members of the initial sample to be drawn multiple times. This results in 10⁴ mass distributions, which are used to derive dN/dM and the corresponding uncertainties.

In the upper panels of Figures 12 and 13 we show the IMF in the power-law form (dN/dM ∝ M^−α), derived from three different sets of isochrones and for distances of 1600 pc and 1400 pc, respectively. In all cases the IMF is well represented by two power laws, with a break at ∼0.4 M_⊙. We fit separate power laws to the mass ranges 0.4–7 M_⊙ and 0.02–0.4 M_⊙, as well as to the one encompassing only the substellar regime (∼0.02–0.1 M_⊙). The least-squares fits are shown as gray lines in Figures 12 and 13, and the resulting slopes α are given in Table 4. For each distance, the α values derived from the three isochrones agree within the uncertainties. The largest variation is seen in the 0.02–0.1 M_⊙ range, where the fit is more sensitive to single-point variations given the small number of points available for it. We can also observe that the uncertainty in distance of the cluster does not significantly affect our results, which again agree within the uncertainties. For a comparison, we overplot the Kroupa segmented power-law IMF, normalized to the total number of objects in the cluster (orange line).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To estimate the effects of the choice of the bin size and location, we repeat the IMF calculation and power-law fits (within the same mass limits) for two additional bin sizes, Δlog(m/M_⊙) = 0.15 and 0.25. We also repeat the calculation with the same bin size as before (Δlog(m/M_⊙) = 0.2) but shift the bin locations by half the bin size. For each run, we combine the obtained slopes from the three isochrones as a weighted average, and we show the results in Table 5 of Appendix B. The first line in the table is the default one and is shown to facilitate the comparison. There is no major effect on the slope of the IMF in any of the considered mass ranges.

In the lower panels of Figures 12 and 13 we plot the IMF in the log-normal form (dN/dlog M), and overplot the Kroupa and Chabrier (Chabrier 2005) mass functions, both normalized to match the total number of objects in the cluster. Unlike the power-law tail at masses >1 M_⊙, the log-normal part of the Chabrier mass function does not seem to represent the NGC 2244 mass function very well. Possible implications of this result will be discussed in Section 5.

For the masses above 0.4 M_⊙, α = 2.12 ± 0.08 is close to the Salpeter slope (α = 2.35; Salpeter 1955), and identical to the one derived for X-ray-selected members with masses >0.5 M_⊙ (Wang et al. 2008). Below 0.4 M_⊙, a single power law with α = 1.03 ± 0.02 describes well both the low-mass stellar and substellar regimes. We can compare this with the values of other clusters and star-forming regions summarized in Table 4 of Mužić et al. (2017). The slope of the low-mass IMF in NGC 2244 is steeper than the one in RCW 38 (located at 1.7 kpc), where we find α = 0.7 ± 0.1 for the mass range 0.02–0.5 M_⊙. The slopes found in nearby star-forming regions at distances of up to 400 pc, typically lie in the range 0.5–1.0. The slope we derive for NGC 2244 at masses 0.02–0.4 M_⊙ is at the high end of this range, but still in agreement with it, taking into account the statistical uncertainties. Furthermore, there are a number of other sources of uncertainty that are more difficult to take into account, such as the errors in age, systematics associated with the evolutionary models, extinction laws, and finally the choice of the mass range included in the fit. The latter is demonstrated by the power-law fits in the substellar regime, which are somewhat shallower than those encompassing also very low-mass stars.

To estimate the ratio between stars and BDs, we use the masses derived in Section 4.3.3 and apply the same method as we did in deriving the IMF, generating 10⁴ mass distributions for each of the 324 member lists. The stellar–substellar mass boundary is set to 0.075 M_⊙, the value at the solar metallicity. We calculate the ratio using two low-mass limits on the BD side, 0.02 and 0.03 M_⊙. On the stellar side, we set an upper limit at 1 M_⊙ as commonly found in the literature, but also calculate the ratio by taking into account all stellar candidates. The results are given in Table 4 for the three different isochrones used in this work. The first column contains the results for the mass ranges 0.03–0.075 M_⊙ for the BDs and 0.075–1 M_⊙ for the stars, which are the limits most commonly used in the literature and the most suitable for comparison to other works. The star-to-BD number ratio in this mass range for the core of NGC 2244 is 3 ± 0.3. We find that the star-to-BD ratio in NGC 2244 is slightly larger than that in RCW 38 (2.0 ± 0.6; Mužić et al. 2017); however, the values are still in agreement within their uncertainties. Furthermore, the star-to-BD ratio in NGC 1333 was found to be 1.9–2.4, and that in IC 348 was 2.9–4 (Scholz et al. 2013). For Cha-I and Lupus 3, Mužić et al. (2015) derive the ratio of 2.5–6.0, although the analysis based on the completeness levels of the spectroscopic follow-up suggests that these ratios might in reality be on the lower side of the quoted span. These numbers are also in agreement with the star-to-BD ratios found in the ONC (${3.3}_{-0.7}^{+0.8}$; Andersen et al. 2008) and NGC 2024 (${3.8}_{-1.5}^{+2.1}$; Andersen et al. 2008).

In Appendix C, we derive stellar surface densities for a number of young clusters and star-forming regions, and in Figure 15 we show the dependence of star-to-BD ratio on this parameter. Different regions are represented by polygons of different colors. The height of each polygon represents either the ±1σ range around the mean value or the range in star-to-BD ratio if given in this form. The width of the polygons has no physical meaning. The filled polygons represent the regions with few or no massive stars, while the dashed ones mark the regions with substantial OB star population. The star-to-BD ratio does not seem to depend on either the stellar density or the presence of OB stars.

To compare our results with expectations from the standard forms of mass function, we perform a simulation in which we draw 310 stars in the mass range 0.02–10 M_⊙ (our estimate of the average number of members and mass limits in the data) from the underlying Kroupa or Chabrier mass functions. We then estimate the star-to-BD ratio counting the sources in the mass ranges 0.03–0.075 M_⊙ and 0.075–1 M_⊙. Repeating this process 10⁴ times, we obtain distributions of star-to-BD number ratio expected from the two standard IMFs. The results are shown in Figure 14, where the gray histograms show the distribution obtained for NGC 2244 from the 324 member lists (total 324 × 10⁴ data points). The dashed blue and solid orange lines show distributions expected from the underlying Chabrier and Kroupa IMFs, normalized to the peak value of the histogram. The six panels of Figure 14 show the star-to-BD ratios derived for the three isochrones and two distances (1600 and 1400 pc). As expected from the IMF, the Kroupa mass function provides a better representation of the NGC 2244 star-to-BD ratio than the Chabrier one, which predicts higher values of the ratio than what we find in NGC 2244.

Zoom In Zoom Out Reset image size

The low-mass end IMF slope in RCW 38 is shallower than that of NGC 2244, while the star-to-BD ratio we derived for RCW 38 is lower than the one derived for NGC 2244. This is the opposite of what is expected, since a shallower IMF slope should in principle result in fewer BDs with respect to stars (i.e., higher star-to-BD number ratio). The reason for this is the way we treated the incomplete mass bins in the survey of RCW 38. The data in the lowest-mass bins were corrected for incompleteness in the IMF calculation, according to the average K-band magnitude of the sources belonging to each bin. However, when calculating the star-to-BD ratio, all BDs are counted in a single bin and therefore this approach could not be applied. In this case we added a "missing" object for each source if its K magnitude was in the incomplete magnitude range. The mass function we derived for RCW 38 is consistent with a slightly higher star-to-BD ratio (∼4). The value is still within the range expected from other star-forming regions, but we note that the way the incompleteness is treated is clearly another important source of uncertainty. In this work, however, we do not encounter this issue since the analysis was restricted to the levels where our photometry is more than 90% complete.