© The Authors 2023

1. Introduction

Brown dwarfs (BDs), objects with masses below 0.075 M_⊙, populate the mass spectrum between the realms of stars and planets, being mainly characterized by the lack of hydrogen burning in their cores. As is the case for stars, most BDs are born in clusters and associations, making those perfect laboratories where one can study the formation of stars and BDs. One of the most studied observable outcomes of the star formation process is the initial mass function (IMF), which represents the mass distribution at birth. Since the pilot work of Salpeter (1955), the high mass side of the IMF has been found to be reproduced by a single power law with the form dN/dM ∝ M^−α, with α ∼ 2.35, except for some notable exceptions (Bayo et al. 2011; Stolte et al. 2006; Mužić et al. 2017; Hosek et al. 2019).

The IMF flattens below 1 M_⊙ with a turnover at 0.1–0.7 M_⊙ (Chabrier 2003). The IMF on the low-mass side can also be well reproduced by a single power law. In nearby star-forming regions (SFRs), the low-mass side α has been found to range between 0.5–1 (Luhman 2007; Bayo et al. 2011; Peña Ramírez et al. 2012; Alves de Oliveira et al. 2012; Lodieu 2013; Scholz et al. 2013; Mužić et al. 2015; Suárez et al. 2019). At the same time, it has been found that two to five stars are born for every BD in a single SFR (Mužić et al. 2019). The wide range of values found reflects the overall contribution of the uncertainties of all parameters involved in the derivation of the slope of the IMF and the star-to-BD number ratio, rather than variations between clusters. Based on this result and the comparison of other properties of BDs and stars (e.g., kinematic and spatial distributions, protoplanetary disk properties, and multiplicity), it is generally accepted that at least most of the BDs with higher masses form similarly to stars (Luhman 2012), from the turbulent fragmentation of the molecular cloud (Padoan & Nordlund 2004). However, due to the intrinsic faintness of BDs, most of the environments where the substellar content has been studied in detail are located nearby and are characterized by relatively loose groups of low-mass stars, with no or very few massive stars. The current BD formation theories predict that the production of BDs could be affected by the environment where they are born. The production of BDs could be increased in environments with high gas and/or stellar densities (Bonnell et al. 2008; Jones & Bate 2018), which may also favor BD formation through ejections (Bate 2012). The presence of massive stars can also increase the efficiency of BD formation through photoevaporation of the outer layers of a forming protostar (Whitworth & Zinnecker 2004) or by fragmentation in their massive disks (Stamatellos et al. 2011; Vorobyov 2013).

In order to probe environments whose properties in terms of stellar density and number of OB stars are more extreme than that found in nearby SFRs, we have a started a program to study a sample of massive young clusters. In Mužić et al. (2017), we studied the core of RCW 38, a young (1 Myr), embedded cluster located at a distance of 1.7 kpc. RCW 38 is characterized both by high stellar densities (with a core stellar surface density Σ ∼ 2500 pc⁻² that is twice as dense as the ONC and more than ten times denser than NGC 1333) and by a substantial population of massive OB stars (∼60 OB stars, Wolk et al. 2008; Winston et al. 2011). We found the IMF of the cluster to be consistent with nearby SFRs (α ∼ 0.7 for the range 0.02–0.5 M_⊙). In Mužić et al. (2019), we studied the core of NGC 2244, a young (2 Myr, Hensberge et al. 2000; Bell et al. 2013; Wareing et al. 2018) cluster located at a distance of 1.5 kpc. It is characterized by a low stellar density (core stellar surface density Σ ∼ 30 pc⁻², similar to the sparse nearby SFRs) and a rich population of OB stars (> 70 massive stars, Román-Zúñiga & Lada 2008). These OB stars are presumed to be responsible for the evacuation of the gas in the central region of the cluster (Román-Zúñiga & Lada 2008), and may have also eroded the protoplanetary disks around stars in their vicinity (Balog et al. 2007). We found the IMF to be on the high end of BD production while still consistent with nearby SFRs (α ∼ 1 for the range 0.02–0.5 M_⊙). The membership analysis in Mužić et al. (2019) was carried out in a statistical way, so, in order to confirm those results, we need to perform a spectroscopic follow-up. Spectroscopy is a fundamental step in IMF studies, where the membership of cluster candidate members can be confirmed, removing one of the sources for uncertainty in the derivation of the IMF. Thanks to the low extinction and distance to NGC 2244, spectroscopy of candidate members is feasible into the substellar regime with current facilities.

In this work, we present a spectroscopic follow-up of 85 NGC 2244 faint candidate members (∼0.03–0.4 M_⊙). In Sect. 2 we present the spectroscopic observations and data reduction, a new reduction of Spitzer mid-infrared data, and a compilation of all the known stellar members of the core region of the cluster. In Sect. 3 we derive the spectroscopic properties of the candidate members with low-mass signatures and evaluate the youth of the cluster candidate members using gravity-sensitive spectral indices. We also estimate the infrared excess and derive the masses of all the members of the cluster. In Sect. 4 we derive the IMF and star-to-BD number ratio and evaluate the disk fraction into the BD regime. In Sect. 5 we discuss the effect that massive stars may have had in the formation of BDs in NGC 2244.

2. Data set

2.1. Near-infrared spectroscopy

Mužić et al. (2019) analyzed the IMF of the core region of NGC 2244 into the BD regime using deep photometric observations. These observations were taken with the Flamingos-2 near-infrared (NIR) camera at the Gemini-South telescope in the JHK bands on a field of ∼2.4 pc² (black dots in Fig. 1, we refer to these observations as the Flamingos-2 catalog hereafter). We use the same ID notation to refer to the candidate members of NGC 2244¹. The 90% completeness limit of the Flamingos-2 catalog is J = 21.7 mag. Here we present the spectroscopic follow-up of 135 sources too faint to have proper motions from Gaia (J > 16.5 mag, see Mužić et al. 2022) down to J ∼ 19.8 mag (limited by the feasibility of the observations).

Fig. 1. Color magnitude diagram (CMD) of the core of NGC 2244. We show the Flamingos-2 catalog (Mužić et al. 2019) with black dots, objects observed with KMOS are represented by orange empty diamonds, KMOS members as solid orange circles, KMOS possible members as blue solid circles, VIMOS observed objects as empty gray diamonds, and members from VIMOS, Mužić et al. (2022) and Meng et al. (2017) as gray circles. The black solid line represents the 2 Myr isochrone shifted to the distance of 1500 pc, and the dashed line represents the same isochrone reddened by the mean extinction toward NGC 2244 (i.e., AV = 1.5 mag, Mužić et al. 2019). The isochrone shown comes from the PARSEC stellar models (Bressan et al. 2012) down to 1 M⊙, and BT-Settl below 1 M⊙ (Baraffe et al. 2015), that is isochrone 1 in Mužić et al. (2019). The right panel represents the inset of the CMD to the region populated by the KMOS observations. The black arrow represents the effect of AV = 1.5 mag.

2.1.1. Observations

The observations were carried out using the K-band Multi Object Spectrograph (KMOS, Sharples et al. 2013), located at the Very Large Telescope (VLT), under program 0102.C-0682(A) (PI: K. Muzic). KMOS performs integral field spectroscopy and has 24 arms, each with a squared field of view (FoV) of 2.8″ × 2.8″. The total FoV of the instrument is 7.2′ in diameter. KMOS was operated in service mode using a nod to sky ABA configuration, where A is the target and B is the sky. Between each exposure a small dithering was performed. In order to optimize the integration time, we divided the observations in two groups according to their brightness: bright (J < 18.5 mag) and faint (J > 18.5 mag). For the bright set of observations, 10 exposures with 90 s detector integration time (DIT) each were obtained, whereas for the faint set, 20 exposures with DIT = 120 s were taken. The spectra were obtained using the HK filter (∼1.5–2.4 μm), which has a mean spectral resolving power of ∼1800. The sources were randomly selected from the complete Flamingos-2 catalog with J = 16.5 − 19.5 mag to maximize the number of arms allocated in each observation block. This process was optimized using the KMOS ARM Allocator (KARMA) software. In this process we also allocated one brighter source (#770) to a free arm. Together with the target spectra, all the calibration frames needed for a successful data reduction were also obtained: dark, flat, arc-lamp, standard star and telluric standard stars.

2.1.2. Data reduction

The data reduction was performed using the KMOS pipeline (SPARK, Davies et al. 2013) in the ESO Reflex automated data reduction environment (Freudling et al. 2013). The pipeline produces calibrated 3D cubes, by performing flat field correction, wavelength calibration, sky subtraction and telluric correction (using Molecfit, Smette et al. 2015). However, we noticed that the quality of the subtraction of the OH sky lines (telluric emission features) in the data products was poor, especially for the cases with lower signal-to-noise ratio (S/N). We found SkyCorr (Noll et al. 2014) to improve the OH sky line residuals in the H-band and the signal at the end of the K-band and in the ∼2–2.05 μm region where there is a telluric H₂O absorption band. Therefore, we used the instrument’s pipeline to perform the entire data reduction except for sky subtraction. Since the objects were not discernible in the data cubes without sky subtraction, we used the pipeline sky subtracted data cubes to define the spectrum extraction mask. We fit a 2D Gaussian to each individual exposure of the pipeline sky-subtracted data cubes. The 2D Gaussian was used to define a mask with the weights of the extraction at each spaxel. The mask was applied to the non-sky subtracted data cubes, as well as to the corresponding dithered sky exposures. Then, we performed the sky subtraction using SkyCorr, the input being the extracted object and sky 1D spectra. We found that SkyCorr performed better when the H and K-band spectra are treated separately, removing the region in between the bands (1.8–1.95 μm) that is dominated by telluric absorption. Finally, we combined all the exposures of each object with a median approach. The noise of the spectra at each wavelength was derived as:

$$\begin{array}{c}\hfill N=\frac{\sigma }{\sqrt{n}},\end{array}$$(1)

where N is the noise, σ is the standard deviation of all the exposures at each wavelength, and n is the number of exposures.

We inspected whether the signal obtained in each spectrum is enough to perform any analysis. An initial visual inspection of the spectra confirmed that 50 out of the original 135 spectra had a low S/N (S/N < 4 in H-band) which prevented any meaningful analysis, and were therefore excluded from the data set. Based on the brightness of the remaining sources and the distance and average extinction of NGC 2244 (Mužić et al. 2022), the masses of the observed candidate members should roughly correspond to 0.03–0.4 M_⊙.

2.2. Mid-infrared photometry

We performed a new reduction of the available Spitzer data of the NGC 2244 region. We retrieved all the images obtained with IRAC (Fazio et al. 2004) on board the Spitzer Space Telescope (Werner et al. 2004) from the Spitzer Heritage Archive within a radius of 1.5° centered on NGC 2244. Table 1 gives the list of program ID, principal investigator and associated exposure times and number of frames (i.e., the number of BCD). Our reduction began from the S18.25.0 pipeline-processed artifact Corrected Basic Calibrated Data (CBCD). We then combined these into mosaics using the recommended version 18 of MOPEX (MOsaicker and Point source EXtractor) provided by the Spitzer Science Center using the standard parameters (see the MOPEX User’s Guide for details on the data reduction). We chose to process the short 0.4 s and long 10.4 s frames separately and to produce a mosaic for each data set. This allowed us to nicely cover the largest possible dynamical range. The presence of a bright extended nebulosity compromises the detection of sources and measurement of their photometry as the background estimations implemented in SExtractor (Bertin & Arnouts 1996) and APEX (the photometry package that is part of MOPEX) are not optimised to deal with such variable extended emission. We therefore applied the nebulosity filter described in Bertin et al. (in prep.) to the pipeline produced mosaics. Sources were then detected using SExtractor, and their aperture photometry was then extracted using APEX. Apertures of 1$\stackrel{″}{.}$2, 1$\stackrel{″}{.}$8, and 2$\stackrel{″}{.}$4 radii and a sky annulus between 14$\stackrel{″}{.}$4 and 24$\stackrel{″}{.}$0 were used to measure the photometry in the four IRAC bands. Finally, the recommended aperture corrections were applied to the fluxes.

Finally, we cross-matched the extracted IRAC catalog with the Flamingos2 catalog. We visually checked that the correspondence between the NIR and mid-infrared sources is correct. In the few cases were the IRAC photometry did not correspond to a single object in the Flamingos-2 catalog we rejected their IRAC photometry.

2.3. Previous members

We compiled all the previously reported candidate members encompassing the core of NGC 224, that is the region covered by the Flamingos-2 catalog. Mužić et al. (2022) built a catalog of probable members of NGC 2244. The members were obtained by applying the Probabilistic Random Forest algorithm using as input optical to NIR photometry as well as proper motions from Gaia EDR3. 85 of these candidate members are in the Flamingos-2 catalog. Mužić et al. (2022) also presented optical spectroscopy of 501 candidate members of NGC 2244. The objects were observed using VIMOS/VLT. The spectral type (SpT) and visual extinction (A_V) were derived for the 327 objects that presented late SpT spectral shape by comparison with spectral templates. Strong fringing in the VIMOS spectra longward of 8000 Å prevented an analysis of gravity-sensitive absorption lines in that range from being performed, leaving Hα as the only signature of youth. The equivalent width of Hα was measured and compared with accretion thresholds defined in the literature (Barrado y Navascués & Martín 2003; Fang et al. 2009). If a source had Hα emission above these thresholds, this emission is assumed to be associated with accretion, which is a signature of youth. But if the Hα emission is below this threshold, then this emission is also consistent with chromospheric activity and hence not a definite signature of youth. We found 11 candidate members with accretion-related VIMOS Hα emission, 9 of which are in the Mužić et al. (2022) member list. The remaining two candidate members are below or close to the magnitude completeness limit of Mužić et al. (2022), which is most probably why they were not selected as members there.

Meng et al. (2017) presented a catalog of probable members of NGC 2244 based on X-ray data (Wang et al. 2008) and infrared excess from Spitzer (Balog et al. 2007). There are 84 sources from this member list in the Flamingos-2 catalog, 23 of which are not present either in Mužić et al. (2022) or have accretion-related Hα emission in their VIMOS spectra. Most of these 23 sources are below or close to the magnitude completeness limit of Mužić et al. (2022) as well. Eight of these sources have a VIMOS spectrum, seven of which present a spectroscopic shape typical of low-mass stars and BDs. The remaining source (#734) is brighter and its spectrum would not show low-mass signatures. We visually examined the colors, proper motions and parallaxes of the 23 sources and only found three to be inconsistent with membership to the cluster (#146, #623, #754). The remaining 20 sources are considered candidate members of the cluster. Overall, we find 107 candidate members of the cluster, all of them are shown as gray circles in the color-magnitude diagram in Fig. 1.

3. Analysis

In this section we present the membership assessment (Sects. 3.1.4 and 3.3) and mass derivation of all the NGC 2244 candidate members. We also corrected the census of cluster members for the fact that we did not observe with KMOS all the sources in the faint regime (see Sect. 3.5).

3.1. Near-infrared spectroscopy

In this section we present the analysis of the KMOS NIR spectra. First, we visually inspect all the spectra looking for well-known features present in the NIR spectrum of low-mass stars and BDs. We derive the SpT and extinction of the sources that present these features by comparison with spectral templates. We validate the SpT measurement by comparing it with the SpT coming from spectral indices. Finally, we assess the membership to the cluster of these candidate members using three gravity-sensitive spectral indices. These indices indicate whether the studied objects have a low surface gravity and are therefore young.

3.1.1. Visual inspection

Based on the expected mass range of the observed candidates (∼0.03–0.4 M_⊙), the members within the observed candidates should roughly have SpT ≥ M2. The HK bands spectra of cool dwarfs present distinct characteristics that make them visually identifiable. We therefore perform a visual analysis of the KMOS spectra looking for features characteristic of late SpT objects. The HK bands of cool dwarfs are dominated by broad H₂O absorption bands at both red and blue edges of the H and K bands. These features start to dominate at SpT ∼ M4 and their strength increase with decreasing SpT (Cushing et al. 2005). The NIR spectrum of cool dwarfs is also shaped by numerous absorption features coming from different species (e.g., CO, NaI, KI, CaI), but they are not well resolved in the KMOS spectra due to their noisy nature. Therefore we are only able to visually inspect the broad shape of the H and K bands. Through visual inspection we find 42 (49% of the targets) to have late SpT spectroscopic shape. The rest of the objects are not discussed further since they present a flat spectrum not in accordance with the cool dwarf shape discussed above. Low-mass members would present a flat spectrum, and therefore be missed in this step, if they are at an earlier evolutionary class (i.e., Class I). However, NGC 2244 has a very low fraction of Class I sources (Balog et al. 2007; Poulton et al. 2008), so the chance of having observed a Class I source is very low. The contaminants will then most probably be background reddened high-mass stars and extra galactic sources.

3.1.2. Comparison with spectral templates

For the 42 late SpT candidate members, we derive their SpT and extinction simultaneously by comparison with spectral templates. We use the young (< 10 Myr) spectral templates from Luhman et al. (2017). The templates are defined for half-integer SpTs from M0 down to L0, complemented with L2, L4 and L7 templates. Each object is directly compared with all the spectral templates for a wide range of extinction values: A_V = 0–12 mag with a step of 0.2 mag. We use the extinction law from Cardelli et al. (1989) with R_V = 3.1 (Fernandes et al. 2012) to model the effect of extinction. The object spectra were resampled to match the spectral template wavelength grid. The comparison is made for the H and K bands neglecting the telluric region in between them (1.8–1.95 μm). All the spectra are normalized at 1.66 μm prior to the analysis, and to each comparison we add a grid of wavelengths at which the template is normalized in order to help find the best fit. The goodness of fit of each comparison is evaluated using the reduced χ² minimization:

$$\begin{array}{c}\hfill {\chi }^2=\frac{1}{N-m}{\displaystyle \sum _{i=1}^N\frac{{(O_i-T_i)}^2}{{\sigma }^2}}\end{array}$$(2)

where O is the object spectrum, T the template spectrum, σ is the noise of the observed spectrum, N the number of data points, and m the number of fitted parameters (m = 3). An example of the output distribution of χ² is shown in Fig. 2 (trimmed at L0 for visualization purposes), where darker regions represent lower χ². We also show in orange the best fit SpT and extinction, and in green the SpT derived using SpT indices evaluated for the entire extinction grid (see Sect. 3.1.3). Cool dwarfs spectra typically show some degree of degeneracy between the spectral type and the extinction (Luhman et al. 2017), which is also evident in this case.

Fig. 2. Probability map of SpT and AV from comparison with spectral templates. Darker regions represent lower values of χ2 from the comparison with the young templates for object #40. We show the best fit SpT and AV in orange. The green line represents the SpT derived from a combination of various spectral indices (see Sect. 3.1.3) for the entire AV grid. At the right and bottom panels we also show the χ2 histogram of SpT and AV.

For each object we associate the SpT and A_V of the template and extinction configuration that minimizes the χ² in the template fitting process. We also performed a visual inspection of all the fits. We noticed that in the case of object #439, the second best fit (M8) reproduced better the blue edges of both the H and K bands than the best fit (M5.5). The spectrum of this object presents a lower S/N and that makes our methodology less precise. We therefore adopted M8 as the SpT of object #439. In Fig. 3 we show the KMOS spectra together with the best fit template de-reddened by the best fit extinction except for #439 as discussed above. The spectra are color-coded by their final membership assessment (see Sect. 3.3): members (orange), possible members (blue) and nonmembers (gray). The distribution of the retrieved SpT and A_V values is represented in Fig. 4. The objects in our sample show SpTs between M2.5 and M8.5. There is a drop in low-mass sources with SpT M7, which could be a bias of the observations, or a real feature associated with a known dip in the luminosity function of SFRs between M7 and M8 (Dobbie et al. 2002). This feature was explained to be possibly caused by the formation of dust in the atmospheres at these temperatures. Most objects present an extinction below 4 mag, and the median extinction is 2.3 mag for the entire low-mass sample, dropping to 2 mag if only the bona fide members are considered. These values are slightly larger than the previously derived mean extinction of the cluster (∼1.5 mag, Mužić et al. 2019).

Fig. 3. All the KMOS spectra of objects consistent with low-mass (sub)stellar nature, together with the best fit young templates (black) de-reddened by the best fit extinction. The selected SpT and AV (in magnitude units) are also shown. The color-coding represents the final membership of the targets: orange are members, blue are possible members and gray are nonmembers.

Fig. 4. Histogram of SpT and AV of the objects observed with KMOS derived by direct comparison with spectral templates. We also show the histograms of the members (orange) and possible members (blue).

3.1.3. Spectral type indices

SpT indices are calibrated to provide a straightforward measurement of the SpT. Their usage is specially interesting when the wavelength range available is constrained to one or two spectral bands, as it is the case here. Using indices for SpT derivation has one main drawback, which is that they generally are extinction dependent, and therefore extinction needs to be derived by a different method. We de-redden the spectra using the extinction value derived in Sect. 3.1.2. We derive the SpT with the method described in Almendros-Abad et al. (2022) using the indices that fall in the KMOS spectra wavelength range: WK (Weights et al. 2009), H2O-2 (Slesnick et al. 2004), sH2O^K (Testi et al. 2001) and TLI-K (Almendros-Abad et al. 2022). In Fig. 5 we compare the SpT derived here with the value derived in Sect. 3.1.2. We observe that both measurements are in excellent agreement, generally lying within 1 subtype.

Fig. 5. Comparison of the spectral templates best fit SpT (Sect. 3.1.2) with the SpT coming from spectral indices (Sect. 3.3). The objects are color-coded based on the final membership discussed in Sect. 3.1.4: members (orange circles), possible members (blue squares) and nonmembers (black crosses). The solid line represents the 1:1 SpT relationship, and the dashed lines represent a ±1 SpT deviation from the 1:1 relationship.

3.1.4. Youth evaluation

For a late SpT candidate to be confirmed as a member of a young cluster it needs to present some kind of youth indicator in its spectrum. Young objects are still contracting and present a larger radii compared with old objects of the same temperature. As a consequence, young objects possess lower surface gravity than their older counterparts. The NIR spectra of cool dwarfs present several gravity-sensitive features that allow a distinction between field, giant and young low-mass objects (Gorlova et al. 2003; Allers & Liu 2013; Almendros-Abad et al. 2022). However, the only feature that can be inspected in our KMOS spectra is the broad band shape of the H-band. Low surface gravity objects present a more triangular and sharp shape of the H-band than field objects (Lucas et al. 2001). This feature has been explained by a decrease in H₂ collision induced absorption with pressure (Kirkpatrick et al. 2006). Using machine learning models, Almendros-Abad et al. (2022) showed that the H-band peak is indeed the most important feature in the NIR spectra when it comes to surface gravity classification.

There are two gravity-sensitive spectral indices defined to model the change of the H-band peak with surface gravity that fall in the KMOS wavelength range: the HPI (Scholz et al. 2012) and TLI-g (Almendros-Abad et al. 2022) indices. Although originally defined for derivation of SpT, Almendros-Abad et al. (2022) showed that the HPI index also presents a gravity-sensitive behavior. The TLI-g index has been shown to be able perform an almost complete separation between young and field cool dwarfs for a wide range of SpT (M0-L3). Additionally, we looked for other SpT indices evaluated in Almendros-Abad et al. (2022; see their Figs. C1, C2), that present a good overall separation of the gravity classes and also fall in the KMOS wavelength range. We find the WH index (Weights et al. 2009) to have a pronounced separation between the gravity classes from very early SpTs (see Fig. 6). We use these three gravity-sensitive indices to perform the membership classification. For the three indices we derive the decision boundary that best represents the separation between the “young” and “field” sources using the data set presented in Almendros-Abad et al. (2022; see Appendix A for further details).

Fig. 6. Gravity-sensitive indices used in Sect. 3.1.4 to evaluate the youth of the KMOS sources compared with the SpT. Dots represent the data set from Almendros-Abad et al. (2022) color-coded according to the age class of each object (young in orange, mid-gravity in blue and field in gray). The black solid line represents the decision boundary between young and field objects (see Appendix A). The objects with KMOS spectra are color- and symbol-coded as in Fig. 5. We also show with a black arrow the effect that an extinction of AV = 5 mag would have on the values of each index.

In Fig. 6 we show the three gravity-sensitive indices as a function of SpT for the KMOS sample and the Almendros-Abad et al. (2022) data set together with the derived decision functions. We perform the membership classification as follows: if an object is classified as low-gravity in at least two indices, it is labeled as member (Y); if it is classified as low gravity only in one index, it is labeled as a possible member (?); it is labeled as a nonmember otherwise (N, see Table 2). With this methodology we aim to keep all the objects that show any hint of having a low surface gravity. In Fig. 6 the objects with KMOS spectra are color-coded based on this surface gravity criteria: young members as orange circles, possible members as blue squares and nonmembers objects as empty black crosses. We observe that four of the sources do not display enough S/N in the KMOS spectrum to rely on the measurement of the spectral indices (#77, #405, #439, #567). We classified these four sources as possible members. Overall, we find 25 bona fide members and 10 possible members. The distribution of the extinction of the possible members is not different to that of the bona fide members, all having A_V < 4 mag (see Fig. 4). There is therefore no clear evidence based on the reddening to distinguish them from foreground contaminants.

3.2. Infrared excess

Young objects with protoplanetary disks present excess emission at infrared wavelengths coming from the inner part of the disk that is heated by the central object. This feature can be mimicked by other types of objects (see e.g., Gutermuth et al. 2009), but if the (sub)stellar nature of the source is confirmed, then the infrared excess is a definite proof of the object’s youth. While this feature is a very strong indicator of youth, it will not uncover full populations of young objects, since not all young objects have disks (Haisch et al. 2001).

We now classify all the candidate members of NGC 2244 based on their infrared excess. We define the excess in the four IRAC channels from the K_S – IRAC colors. We first derive the K_S – IRAC photospheric colors as a function of SpT for a compilation of known diskless members of Cha-I (Luhman et al. 2008; Esplin et al. 2017), Upper Sco (Esplin et al. 2018), ρ-Ophiuchus (Esplin & Luhman 2020), Taurus (Esplin & Luhman 2019) and Corona Australis (Esplin & Luhman 2022) SFRs (see Appendix B). We observe that in the four IRAC colors an excess of 0.25 mag separates most (> 97%) of the diskless sources (see Appendix B). We therefore define 0.25 mag as the boundary to classify an object as having excess in a certain color.

We explore the infrared excess of all the NGC 2244 candidate members that have a SpT measurement, that is the candidate members with KMOS or VIMOS spectra. Mužić et al. (2022) derived the T_eff together with extinction from fitting the spectral energy distribution (SED) using the Virtual Observatory SED Analyzer (VOSA; Bayo et al. 2008). We convert the T_eff to SpT using the Pecaut & Mamajek (2013)T_eff–SpT scale (see Sect. 3.4 for justification of the choice of scale). We therefore also evaluate the infrared excess of the Mužić et al. (2022) candidate members. We classify an object as having excess when it presents excess in at least two colors. We require two colors for the excess classification in order to discard objects that are marginally detected in a single band. With this criteria we may be missing sources with excess only at longer wavelengths (IRAC3 or IRAC4) that are associated with disks at later stages of their evolution such as transitional disks. Additionally, the color excess is affected by the determination of the SpT, which is specially important at later SpTs (> M0) where the photospheric colors have a steep relationship with SpT. Therefore, we also visually inspect the SED of all the sources classified as having infrared excess and those that present excess only in IRAC3 or IRAC4. For each source we compare their SED with the BT-Settl atmospheric model (Allard et al. 2012) of the T_eff and A_V used in Sect. 3.4. We find 6 candidate members to have excess only in IRAC4 (#33, #49, #95, #197, #726, #763), and we classify them as having excess. We also find 4 candidates with color excess to be diskless when inspecting their SEDs (#205, #237, #737, #739). In Fig. 7 we show the color excess as a function of SpT for all the members of the cluster in the four IRAC bands. The locus for sources with full disks and those that are diskless in the four colors are represented with orange and blue contours respectively. The NGC 2244 candidate members that are classified with infrared excess are highlighted with a big black empty circle.

Fig. 7. KS – IRAC color excesses for the candidate members of NGC 2244 with SpT or Teff measurement. Orange circles represent KMOS bona fide members, blue squares KMOS possible members and gray diamonds represent VIMOS objects and Mužić et al. (2022) members. Big black circles represent the objects that we found to have infrared excess. The locus of sources from the literature with full disks is represented with orange contours, and the locus of sources without disks is represented with blue contours (see Appendix B). Black dashed line represents the 0.25 mag excess boundary.

3.3. Final census of members

In this section we perform the final assessment of the census of members of the core region of NGC 2244. So far, we have found 107 candidate members from the literature (see Sect. 2.3) and 25 members and ten possible members from the analysis of the KMOS spectra (see Sect. 3.1.4). Additionally, we find KMOS possible member #77 to present an infrared excess. This object only has IRAC1 and IRAC2 photometry, but it presents strong excess in both bands. We also find seven of the candidate members with VIMOS spectra that did not present accretion-related Hα emission to present infrared excess. Two of the members that have KMOS spectra are also members in both Mužić et al. (2022) and Meng et al. (2017) catalogs (#387, #770), and another two members appear only in the Meng et al. (2017) catalog (#18, #152). KMOS member #18 also has VIMOS spectrum giving the same SpT and extinction (within errors) as in the KMOS spectrum (see Sect. 3.1). The VIMOS spectrum of object #18 does not present accretion-related Hα emission and the object does not present infrared excess either.

We consider a candidate member to be a bona fide member if any of these conditions are met: (1) is a member in Mužić et al. (2022), (2) has a VIMOS spectrum with accretion-related Hα emission and/or infrared excess, (3) has a KMOS spectrum with low-gravity classification and/or infrared excess. In total, we obtain 117 bona fide members. We also consider a candidate member to be a possible member if any of these conditions are met: (1) is a member in Meng et al. (2017) and did not meet any of the bona fide conditions, (2) has a KMOS spectrum with a possible member classification based on the gravity-sensitive indices (see Sect. 3.1.4) and does not present infrared excess. In total, we obtain 22 possible members, one of which (#208) meets both of the possible member conditions. In Table 3 we show how many of these bona fide and possible members come from each member list.

3.4. Derivation of stellar parameters

We derive the stellar parameters (T_eff and mass) of all the candidate members in the core region of NGC 2244 (see Sect. 3.2). For the KMOS members and possible members we derive T_eff together with log g by direct comparison with the BT-Settl models (Allard et al. 2012). We de-redden all the spectra using the A_V derived in Sect. 3.1.2. We fix the metallicity at the solar value leaving only T_eff and log g as free parameters. The model T_eff grid we use is 1500 – 4200 K with a step of 100 K, encompassing the highest and lowest possible T_eff expected in our sample based on their SpT. We constrain the model log g grid to 3.5 and 4 dex, since these are the values expected at the age of the cluster and the mass range of the candidate members based on the Baraffe et al. (2015) theoretical isochrones. The comparison is performed following the same steps of Sect. 3.1.2, but in this case the free parameters are T_eff, log g and the normalization wavelength. We adopt the T_eff from the best fit model.

In Fig. 8 we show the relationship between SpT and T_eff for our sample. We also show the young T_eff–SpT scale from Luhman et al. (2003) with a red solid line, the young scale from Herczeg & Hillenbrand (2014) with a black dotted line and the field scale from Pecaut & Mamajek (2013)² with a gray dashed line. We observe that while the three scales agree well with our results, the scale from Pecaut & Mamajek (2013) agrees better with our results of the model fitting in the M6–M8 range (see Appendix C for further discussion). For the members from Mužić et al. (2022) we take the T_eff and A_V values given there coming from SED fitting. For the members and possible members with VIMOS spectroscopy, we use the T_eff–SpT scale of Pecaut & Mamajek (2013) to convert the SpT into T_eff, and take the A_V given in Mužić et al. (2022) obtained from a template fitting process.

Fig. 8. Teff derived in Sect. 3.4 compared with the SpT derived in Sect. 3.1.2 for the bona fide members (orange circles) and possible members (blue squares). The young Teff scale from Luhman et al. (2003, L03, red solid line), the field Teff scale from Pecaut & Mamajek (2013, PM13, gray dashed line), and the young Teff scale from Herczeg & Hillenbrand (2014, HH14, black dotted line) are also shown.

Using A_V, a distance to the sources of 1500 pc and the photometry from the Flamingos-2 catalog we obtain the absolute J-band magnitude (M_J). Combining M_J with T_eff we construct the Hertzsprung-Russell (HR) diagram of the entire data set (see Fig. 9). We overplot the Baraffe et al. (2015) theoretical isochrones for ages between 0.5 and 100 Myr. We observe that most of the sources lie above the 2 Myr isochrone. Some recent studies suggest that the cluster may be younger than previously thought (between 1 and 2 Myr; Wareing et al. 2018; Mužić et al. 2022; Román-Zúñiga et al. 2023), although these results are still in agreement with an age of 2 Myr.

Fig. 9. Hertzsprung-Russell diagram of all the members of NGC 2244. The color-coding is the same as in Fig. 1. The isochrones (black dashed lines, shown for 0.5, 2, 5, 10 and 100 Myr) and the lines of constant mass (green dashed lines) shown are from Baraffe et al. (2015).

For each object, we associate the age and mass of the closest point of the isochrones interpolated to a finer grid of 1000 by 1000 points in age and mass respectively. The uncertainty of the derived parameters was estimated with a Monte Carlo (MC) simulation by varying T_eff and M_J by their respective uncertainties. For T_eff we assumed an uncertainty of 100 K (model grid step). For M_J we include the photometric errors as well as an uncertainty of 0.25 mag in A_V, added in quadrature. The distance uncertainty is not included in this calculation. The derived mass and age distributions for each object do not typically follow a normal distribution. Therefore, we save the derived distributions, which are later used to draw random samples for the derivation of the IMF. Objects #205 (KMOS bona fide member), #368 (KMOS possible member) and #119 (member of Mužić et al. 2022 with VIMOS spectrum) are subluminous. Both #205 and #119 present infrared excess in the Spitzer data, which could mean that their disk is seen edge-on (Crapsi et al. 2008).

There are ten possible members coming from Meng et al. (2017) without VIMOS or KMOS spectroscopy for which we only have JHK_S photometry. We derive their mass and A_V by de-reddening the objects to the 2 Myr (Baraffe et al. 2015) isochrone (shifted to the distance of NGC 2244, 1500 pc) in the J − K_S/J diagram. In order to obtain the uncertainties of these measurements we perform a MC simulation by varying J and J − K_S by their respective uncertainties. In Table D.1 we provide the properties (coordinates, SpT, A_V and T_eff) of all the members and candidate members in the core of NGC 2244 we have found. This table also includes the origin of their membership candidacy to the cluster, the infrared excess key derived in Sect. 3.2 and the final membership assessment (see Sect. 3.3).

3.5. Completing the census

The KMOS member sample only represents a fraction of all the cluster members in its brightness range. In order to have a complete census we have to apply the obtained confirmation rate (number of confirmed members over number of observed objects) to the objects without spectroscopy in this brightness range. These “new members” together with the bona fide members represent the underlying full population of the cluster. We define the KMOS brightness range to be J = 16.5 − 19 mag. J = 16.5 mag is the completeness limit of the Mužić et al. (2022) study, and J = 19 mag roughly corresponds to 0.03 M_⊙ based on the distance, age and mean extinction of the cluster. However, below J = 18.5 mag (∼0.045 M_⊙) the quality of the KMOS spectra starts to worsen and we only find one bona-fide member (see right panel of Fig. 1). We may therefore have an incomplete census in the J = 18.5 − 19 mag range. There are 143 objects in the Flamingos-2 catalog within the defined brightness range that were not observed with KMOS and are not members in any of the other catalogs.

The confirmation rate of the KMOS sample is found to be 36 ± 2%, which is in good agreement with the contamination rate of ∼60% found in Mužić et al. (2019). However, we observe that the confirmation rate changes significantly at different brightness levels, being ≈70% at J < 17.5 and ≈30% at J > 17.5 mag. We therefore follow a similar approach to that used in Mužić et al. (2019) for the statistical cluster membership (see their Sect. 4.2.2). The CMD is divided in several grid cells with a certain Δcol and Δmag. In each cell grid we calculate the confirmation rate by integrating the contribution of each KMOS object to that cell and then obtaining the fraction of members. Each KMOS object is represented by a gaussian probability distribution, where the width is equal to its uncertainty. The measured in-cell confirmation rates are then applied to the objects in the Flamingos-2 catalog without spectra contained in each cell. Following Mužić et al. (2019) methodology, we vary both the size and position of the grid cells. We have tested a cell size of Δcol = (0.3,0.4) mag and Δmag = (0.5,0.6) 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. Additionally, we perform the analysis on both the J − K_S/J and J − H/J CMDs, resulting in 72 different grid cell configurations. Since this membership extrapolation method is statistical, we obtain 1000 different member lists by bootstrapping the 72 grid cell configurations.

We derive the mass and A_V of all the sources in the KMOS brightness range without spectroscopy. We use the same methodology used in Sect. 3.4 for the Meng et al. (2017) possible members de-reddening to the 2 Myr Baraffe et al. (2015) isochrone in the J − K_S/J diagram. The mass completeness of the combined census is given by the magnitude limit imposed in the membership statistical extrapolation, which is roughly equivalent to 0.03 M_⊙. However, as mentioned previously, this limit may be overestimated. Therefore, we set the completeness limit of the spectroscopic observations to J = 18.5 mag, that is 0.045 M_⊙. While it is not straightforward to assess the completeness limit of a spectroscopic follow up, we opted for setting this rather conservative limit to avoid artifacts due to incompleteness at lower masses.

4. Results

4.1. Initial mass function

The IMF is obtained using the methodology presented in Mužić et al. (2022). The mass of each member is represented by a distribution derived from the different MC simulations. Each of these distributions is smoothed by a Gaussian kernel density estimator (KDE), and used to extract a random sample of N₁ values per object, constructing thus N₁ realisations of the mass distributions for the entire sample. For each of these realisations, we perform N₂ bootstraps, resulting in N₁ × N₂ mass distributions. Each distribution is then binned onto the same grid, and the final value and its uncertainty for each bin is calculated as mean and standard deviation of N₁ × N₂ values. We set N₁ = N₂ = 100. The masses are binned to logarithmic bin sizes of Δlog(m/M_⊙) = 0.2 for masses between 0.03 and 1 M_⊙. The final IMF is derived from the mean and standard deviation of the IMF obtained using the 1000 different member lists. Furthermore, we repeat this procedure for different treatments of the possible members. We weighted the possible members with three different probabilities, that is P = 0, 1, and 0.5. We considered these three scenarios separately for the KMOS and the Meng et al. (2017) possible members, exploring nine different scenarios for the treatment of the possible members.

In order to minimize the dependence of the IMF slopes derived on the choice of the bin size and location, we repeat the IMF calculation 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 location of the bins by half the bin size. In total we derive 36 different IMFs. In Fig. 10 we show the IMF obtained when including the KMOS and Meng et al. (2017) possible members with 0.5 probability of being members in each member list with bin sizes of Δlog(m/M_⊙) = 0.15. In the top panel of Fig. 10 we observe that the mass bin located below the mass completeness limit at 0.045 M_⊙ (red dashed line), is well below the power-law fits we performed and the Kroupa mass function. This confirms the unreliable nature of the membership methodology below this limit. The treatment of the KMOS possible members affects all mass bins below 0.3 M_⊙, specially the last one, where possible members dominate. The treatment of the Meng et al. (2017) possible members affect the mass bins between 0.2 and 0.8 M_⊙.

Fig. 10. IMF of the central region of NGC 2244 considering both KMOS and Meng et al. (2017) possible members with P = 0.5 (see text), represented in equal-size logarithmic bins with size Δlog(m/M⊙) = 0.15. Two different IMF representations are shown: dN/dm (top panel) and dN/dlog(m) (bottom panel). The green dashed and the pink dotted lines in the top panel represent the power-law fits performed between 0.045–0.2 M⊙ and 0.045–0.4 M⊙ respectively. The orange solid line (both panels) represents the Kroupa segmented power-law mass function (Kroupa 2001), the blue solid line (bottom panel) shows the Chabrier mass function (Chabrier 2005) and the gray dash-dotted line (bottom panel) represents a log-normal function with a characteristic mass (mc) of 0.13 M⊙ and standard deviation (σ) of 0.45. The three representations of the IMF are normalized to match the total number of objects in the cluster. The black vertical dashed line at 0.03 M⊙, and the red line at 0.045 M⊙ represent the KMOS survey limit and the true mass completeness of our sample respectively.

In the bottom panel of Fig. 10, we plot the IMF in the log-normal form (dN/d log M) together with the Kroupa (Kroupa 2001, orange line) and Chabrier (Chabrier 2005, blue line) mass functions. Both mass functions are normalized to have the same number of objects as the IMF derived for NGC 2244. We observe that the peak mass of the IMF is lower and the spread smaller than predicted by the Chabrier mass function, whereas the Kroupa mass function reproduces the overall shape of the IMF. We find the log-normal form of the IMF to be well-reproduced by a single log-normal function (gray dash-dotted line) down to the completeness limit of the observations. The shown log-normal function has a characteristic mass (peak of the mass distribution) of 0.13 M_⊙ and a standard deviation of 0.45. This characteristic mass is lower than that found for NGC 2244 in Damian et al. (2021, ∼0.23–0.28 M_⊙). However, the completeness of the observations of NGC 2244 in Damian et al. (2021) is ∼0.07 M_⊙. It is also lower than that found in nearby SFRs (see Suárez et al. 2019, and references therein). A possible reason for this difference is that we have not consider objects with masses larger than 1 M_⊙, which could potentially shift the characteristic mass of the IMF of NGC 2244 toward larger values.

We fit the 36 IMFs with power laws (dN/dM ∝ M^−α) at three different mass ranges: 0.045–0.2 M_⊙, 0.045–0.4 M_⊙, 0.045–1 M_⊙. The fit performed up to 1 M_⊙ produces an α that is larger than any previously reported α value (1.41 ± 0.15). However, we observe that this fit is well above the lowest mass bins, which means that the α slopes derived from this fit are overestimated. We do not discuss these results further since they are not representative of the low-mass population of the cluster. The other two fits are able to reproduce the shape of the IMF in their mass ranges and are shown in the top panel of Fig. 10 as green dashed and pink dotted lines respectively. In Table 4 we present the α slopes obtained in the two mass ranges for the nine treatments of the possible members. Here, the results for the different mass bins have been combined together. We also provide the mean and standard deviation of the α for the two mass ranges explored.

The α slope we found below 0.2 M_⊙ (α = 0.7 ± 0.12) agrees well with the previous result found for NGC 2244 in Mužić et al. (2019) for the 0.02–0.1 M_⊙ mass range (α = 0.75 ± 0.13). While the result we obtain below 0.4 M_⊙ (α = 1.06 ± 0.16) agrees with the results for the 0.02–0.4 M_⊙ mass range in Mužić et al. (2019, α = 1.03 ± 0.03). The 0.045–0.4 M_⊙ slope is steeper than the slope found in Mužić et al. (2017) for RCW 38 (α = 0.7 ± 0.1 for the mass range 0.02–0.5 M_⊙). In nearby SFRs at distances of up to 400 pc, the α slope typically ranges between 0.5–1 (Mužić et al. 2017; Gennaro & Robberto 2020; Suárez et al. 2019, and references therein). The 0.045–0.4 M_⊙ slopes we derived are on the high end of production of low-mass stars and BDs, but they are still in agreement taking into account the uncertainties.

An environmental property that may affect the slope of the IMF is metallicity (Kroupa 2002; Yan et al. 2020; Li et al. 2023). Several studies have found a bottom-light IMF in regions with very low metallicity (Reylé & Robin 2001; Geha et al. 2013; Li et al. 2023). We cross-matched our members with measurements coming from the LAMOST (Luo et al. 2015) and Gaia-ESO (Randich et al. 2022) surveys. We obtain a metalliticy of [Fe/H] = −0.08 ± 0.11 and [Fe/H] = −0.07 ± 0.2 for LAMOST and Gaia-ESO respectively. These results are in agreement with the recent results derived for the entire Rosette nebula from APOGEE of [Fe/H] = −0.1 ± 0.04 (Román-Zúñiga et al. 2023) and with a solar metallicity.

4.2. Star/BD ratio

We estimate the ratio between the number of stars and BDs (star-to-BD number ratio) using the same methodology as we did in deriving the IMF: for each of the 1000 member lists we extract N₁ mass values from the gaussian KDE of the mass distribution, then for each of the N₁ realizations, we perform N₂ bootstraps. We derive the star-to-BD number ratio and its uncertainty as the mean and standard deviation of the 1000 member lists. We place the star-BD boundary at 0.075 M_⊙ (Chabrier et al. 2023). We derive the ratio between 0.045 M_⊙ up to 1 M_⊙, and we do it for the nine different treatments of the possible members studied. We calculate the star-to-BD number ratio in this mass range because it is the range used in most measurements of the ratio in the literature. We find a mean star-to-BD number ratio of 3.6 ± 0.9. However, most studies in the literature calculate the star-to-BD number ratio down to 0.03 M_⊙. Therefore, we derived the star-to-BD number ratio down to 0.045 M_⊙, and in the 0.03–0.045 M_⊙ mass range, we added the contribution from the IMF power-law fits (given in Table 4). We find a mean star-to-BD number ratio of 2.8 ± 0.5 and 2.2 ± 0.3 when correcting for completeness using the IMF fit < 0.2 M_⊙ and < 0.4 M_⊙, respectively. We present these results in Table 5. These values are in agreement with the mean value found in Mužić et al. (2019, 3  ±  0.3) for the core of NGC 2244. They are also in agreement to the star-to-BD number found in RCW 38 (3.0  ±  1; Mužić et al. 2017) and the values typically found in nearby star forming regions (2–5, e.g., Scholz et al. 2013; Mužić et al. 2015; Andersen et al. 2008), but at the same time, they are on the high end of BD production.

We take the stellar densities derived in Mužić et al. (2019) for several star forming regions and their star-to-BD number ratios from the literature and compare them with the values derived in this work. In Fig. 11 we show the values of the star-to-BD number ratio for different star forming regions as a function of the stellar surface density. Regions with high number of OB stars are shown with a hatched pattern. We represent the star-to-BD number ratio obtained when correcting for completeness using the < 0.4 M_⊙ IMF power-law fit with a black hatched polygon. The star-to-BD number ratio does not seem to depend on either the stellar density or the presence of OB stars. However, we observe that NGC 2244 hosts one of the (if not the) richest population of BDs ever found in a SFR.

Fig. 11. Dependence of the star-to-BD number ratio on cluster surface density. 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 number 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 result shown for NGC 2244 is the one calculated correcting for completeness with the IMF power-law fit performed between 0.045 and 0.4 M⊙.

4.3. Distance to OB stars

Here we explore the distribution of projected distances between the OB stars and the stellar and substellar population of NGC 2244. We get the OB stars of the cluster from Román-Zúñiga & Lada (2008) and assume that all members and OB stars are located 1500 pc from us. We then calculate the minimum projected distance of all the members to the OB stars in the cluster. In Fig. 12 we show the projected distance to OB stars cumulative distribution function (CDF) of stars (orange solid line) and BDs (cyan solid line). We observe that BDs are typically closer to OB stars than stars. We performed an Anderson–Darling test to study the null hypothesis that the two samples are drawn from the same underlying population, obtaining a probability of 0.7% that they are drawn from the same distribution, which confirms that the two distributions are different. This result holds when only the bona-fide members are used in the calculation (i.e., not including the possible members), and also when calculating the distance to the O-type stars (i.e., leaving the B-type stars out). In order to test that this result is not an observational bias in the sense that most of the faint sources with good KMOS spectra were by chance in the fields located closer to OB stars, in Fig. 12 we also show the projected distance CDF of the KMOS faint nonmembers (J > 17 mag) and the observed KMOS sources that did not have enough S/N for the analysis. We observe that these sources present a CDF much more similar to that of stars. The Anderson-Darling test confirms that these distributions are indeed very similar, finding a probability of 29.2% and 17.2% that stars and the faint nonmembers or the sources with insufficient S/N are drawn from the same distribution, respectively.

Fig. 12. Cumulative distribution function of the projected distance to the closest OB star of the stars (orange line) and BDs (cyan line) members of the core of NGC 2244. We also show the distribution of nonmembers with J > 17 mag (black line) as well as the KMOS sources with very low S/N that were not considered in any part of the analysis (gray line).

In Fig. 13 we show the spatial distribution of all the members of the core region of the cluster and that of the OB stars coming from Román-Zúñiga & Lada (2008). We observe that the stars are located all over the studied area, whereas the BDs seem to be preferentially located in the central region of the field. The overall asymmetry of both stellar and substellar members seen in Fig. 13 is a consequence of the footprint of the Flamingos-2 observations (Mužić et al. 2019). This asymmetry should not affect our results, since it should have the same effect in the distribution of both stars and BDs. Mass segregation could also be affecting the spatial distribution of the cluster members. Several authors have explored mass segregation in NGC 2244 and found contradictory results (Chen et al. 2007; Wang et al. 2008; Mužić et al. 2022). Mužić et al. (2022) found that some degree of mass segregation may be present in NGC 2244, but it does not seem to be very pronounced.

Fig. 13. Sky distribution of all the members of the core region of NGC 2244. Orange circles represent stars (< 1.2 M⊙), blue circles represent BDs and yellow stars represent the OB stars in the region coming from Preibisch et al. (2008). On the background we show the KS-band image of the center of NGC 2244 taken with Flamingos-2/Gemini-S (Mužić et al. 2019).

4.4. Disk fraction

Due to the particular environment of the cluster (high number of OB stars while keeping nearby SFR-like stellar density), we can probe whether OB stars have had any impact in the evolution of protoplanetary disks in this cluster. We explored the fraction of members of NGC 2244 that present infrared excess (i.e., that have disks) evaluated in Sect. 3.2. We compute the disk fraction as the fraction between the number of sources that present disk excess and the total number of sources with IRAC data in at least two channels (analogous to: Class IIs over Class IIs and Class IIIs). We find a total disk fraction of 39 ± 9% (error coming from binomial probability distribution with Clopper–Pearson method). We find the disk fraction of stars (40  ±  9%) to be higher than that of BDs (31  ±  22%), although the errors in the disk fractions are very large to confirm any real difference. We observe no difference in the distribution of the projected distance to the closest OB star between the disk-bearing and diskless sources, nor do we find any difference when separating the sample in stars and BDs. We find the same result when limiting the analysis to the O-type stars.

The Spitzer data are available in at least two bands for most sources (> 96% of all the candidate members). However, the completeness of the IRAC3 and IRAC4 bands decreases significantly below J = 17 mag. The 50% completeness is located at J = 18.5 mag in IRAC3 (∼0.045 M_⊙), and at J = 17 mag in IRAC4 (∼0.1 M_⊙). This implies that we may be missing some BDs that only present disk excess in IRAC3 or IRAC4 and may account for the difference in the disk fraction we have found between stars and BDs. The same issue may be affecting the derived total disk fraction. Given the infrared excess detection probabilities provided in Table B.1, the possibility to not detect the infrared excess of a source with full disk but no detection at longer wavelengths than IRAC2 is below 30%. Based on this, and taking into account the sources that only present Spitzer detections in the IRAC1 and IRAC2 bands, the disk fraction may be higher by 1.5%, which is well inside the uncertainty of the disk fraction derived. We also do not include in this analysis redder mid-infrared observations (i.e., 24 μm), so we are most likely missing some evolved disks that start dominating the emission at longer wavelengths (see Fig. B.2). However, most of the disk fractions available in the literature are also suffering from the same caveat, either because they do not include young objects with evolved disks in the derivation of disk fractions or because they do not include any 24 μm data in their analysis either.

If we compare the disk fraction of NGC 2244 with that of nearby SFRs (see e.g., Fig. 5 of Michel et al. 2021), we observe that most nearby SFRs of ages ∼2 Myr have disk fractions ≥50%, although our result in NGC 2244 is still in agreement within 1σ with regions such as Taurus or IC 348. Pfalzner et al. (2022) found that protoplanetary disks have median disk lifetimes of 5–10 Myr, and they proposed that the main reason behind the low disk mean lifetimes found until recently (1–3 Myr, Cieza et al. 2007; Fedele et al. 2010), is that they were found targeting clusters located further away having therefore an over-representation of higher-mass stars that are known to disperse their disks earlier (Williams & Cieza 2011; Ribas et al. 2015). However, in this work we restricted the analysis to objects with SpT later than K0 and even extend into the BD regime.

One possible explanation for the low prevalence of disks found in NGC 2244 is the presence of OB stars which would help disperse the protoplanetary disks faster with their strong UV radiation (see review by Winter & Haworth 2022). The population of young stellar objects with disks in NGC 2244 has been previously studied using Spitzer data (Balog et al. 2007; Poulton et al. 2008). Balog et al. (2007) found NGC 2244 to have a disk fraction of 44.5%, while they found a lower disk fraction of 27% when considering only the stars closer to O-type stars (< 0.5 pc). This result hints that OB stars may be affecting the evolution of protoplanetary disks closest to them. Our result agrees within errors with both results found in Balog et al. (2007), but we did not find a lower disk fraction for the members closest to the OB stars.

5. Summary and conclusions

In this work we have presented an analysis of the low-mass population of the core of NGC 2244 (central 2.4 pc²) using new, deep NIR KMOS/VLT spectroscopy of faint candidate members. We have also presented a new reduction of Spitzer data in the region that goes into the BD regime. In the stellar regime, the members come from a catalog built using optical to NIR photometry and Gaia proper motions (Mužić et al. 2022), VIMOS/VLT optical spectroscopy, and a catalog built using infrared excess and X-rays (Meng et al. 2017). We performed the analysis of the KMOS spectra, including the derivation of the stellar parameters (SpT, A_V) and the assessment of their membership to NGC 2244 using gravity-sensitive spectral indices. We have used the Spitzer photometry to evaluate the infrared excess of all the candidate members of the cluster. From the KMOS spectra, we find 26 new low-mass members and nine possible members (< 0.4 M_⊙) of NGC 2244, including the first spectroscopically confirmed free-floating BDs beyond 1 kpc, to the best of our knowledge.

We derived the IMF of the core region of NGC 2244, and find it to be well represented by a single power law between 0.045–0.2 M_⊙ and 0.045–0.4 M_⊙. The slope of the IMF is heavily dependent on the fitted mass range, with typical values of α = 0.7 − 1.1. These α slopes are at the high end, but still in agreement with those found in nearby SFRs (α = 0.5 − 1). We find the IMF to also be in agreement with the results found in Mužić et al. (2019) for the same region using only photometry. While our census is complete down to ∼0.045 M_⊙, most studies in the literature calculate the star-to-BD number ratio down to 0.03 M_⊙. In order to make a proper comparison with other works, we corrected our census using the IMF power-law fits and, in so doing, we obtain a star-to-BD number ratio between 0.03 and 1 M_⊙ of 2.8 ± 0.5 and 2.2 ± 0.3 by correcting for completeness using the IMF fits below < 0.2 M_⊙ and < 0.4 M_⊙, respectively. These results are in agreement with previous results in nearby SFRs (2-5), although they are on the low end of typical values, confirming the rich population of BDs in the core of the cluster. As discussed in Mužić et al. (2019), the BD population of the massive young cluster NGC 2244 seems to agree very well with the results coming from nearby SFRs. If there are any variations in the low-mass mass function, they must be subtler than the error bars allow us to discern.

We found that BDs in NGC 2244 are preferentially closer to the OB stars than stars. Whitworth & Zinnecker (2004) proposed that BDs can form due to the photoevaporation of the outer layers of the forming protostar by a nearby OB star. This mechanism is considered inefficient (a massive core is necessary to form a single BD) and cannot be a dominant process in BD formation (many BDs are observed in nearby SFRs where there are no OB stars), but it could still influence the production of BDs in centers of massive clusters rich in OB stars such as NGC 2244. BDs that formed via this mechanism are initially in the vicinity of OB stars, so our finding seems to suggest that OB stars have stimulated the formation of BDs in the central region of NGC 2244. However, in this scenario, at least a fraction of the BDs we have found must have formed following the typical BD formation process in nearby SFRs. These BDs should have been distributed randomly along the cluster, and not preferentially closer to the OB stars. Given the current velocity dispersion of the cluster (∼1.4 km s⁻¹, Mužić et al. 2022) a typical NGC 2244 member takes ∼1.1 Myr to cross the core of the cluster (i.e., the region studied in this work, 2.4 pc²). This implies that the initial spatial configuration of the members of NGC 2244 can be significantly different than the one we have observed. If there are BDs in NGC 2244 that have formed thanks to the effect of nearby OB stars, we would expect the BD population of NGC 2244 to be richer than that found in nearby SFRs where there are no OB stars. However, the range of sizes of the known BD populations is very wide, and the increase in BD formation due to OB stars could still keep NGC 2244 in agreement with previous results. A possible explanation for the lack of BDs in some regions of the cluster is that these regions present higher extinction, which prevented us from finding BDs. However, the spatial distribution of the faint objects (J > 17 mag) with KMOS spectra that we find to be contaminants in Sect. 3.1.1 is similar to that of stars. We also do not find evidence of systematic higher extinctions in the regions where we do not find BDs. Other important caveats of this result are that we are missing the line-of-sight positions of the cluster members and that the total number of BDs in this study is relatively low.

Another scenario proposed where the formation of BDs is facilitated by the presence of OB stars is by the gravitational collapse of globulettes. Globulettes are very small clumps of gas found in HII regions that would give birth to objects below the hydrogen-burning limit if they undergo gravitational collapse. These globulettes can be seen as a scaled-down version of the well-known proplyds (O’dell et al. 1993; Johnstone et al. 1998), although they do not always have an object inside (Mäkelä et al. 2014; Grenman & Gahm 2014). It has been proposed that these objects may be the detached or fragmented remnants of pillars located on the edges of HII regions (Gahm et al. 2007; Haworth et al. 2015) and a rich population of globulettes has been found in the Rosette nebula (Gahm et al. 2007, 2013; Mäkelä et al. 2014). However, Haworth et al. (2015) found that globulettes alone are in principle stable against gravitational collapse, until they collide with the HII region shell. On some occasions, this impact may trigger the gravitational collapse of the globulette, therefore forming a BD. Following this reasoning, the collapse of these objects would occur far from the OB stars and would not be a possible scenario to explain our results. Another issue with this scenario is that the mass distribution of globulettes peaks below the deuterium-burning limit (∼13M_Jup, Gahm et al. 2007, 2013). If these objects undergo gravitational collapse, most of the objects formed would be on the planetary-mass domain.

Here we discuss other proposed formation scenarios of BDs and whether they can explain the observed spatial distribution of BDs in NGC 2244. Reipurth & Clarke (2001) proposed that the dynamic interaction between protostars or fragments within a massive core may result in the ejection of the smallest fragments. If the ejected bodies have masses below the hydrogen-burning limit, they may become a BD. The ejection velocity of these objects would be on the order of 3 km s⁻¹, traveling a distance of 6 pc in 2 Myr (age of NGC 2244), which is inconsistent with our results. Bate (2012) found that the velocity dispersion of the ejected clumps should also be larger than 1 km s⁻¹. BDs could also form in protoplanetary disks of stars (e.g., Stamatellos & Whitworth 2009; Kratter et al. 2010; Vorobyov 2013). Basu & Vorobyov (2012) showed that clumps in protoplanetary disks could potentially be ejected by internal dynamical interactions between the clumps, and they would be ejected with lower velocities (∼0.8 km s⁻¹); however, these velocities would still kick these BD clumps far from their original positions by the age of NGC 2244. Ejection velocities below 0.3 km s⁻¹ are needed in order to keep BDs that were formed by any of these mechanisms close to their birth position. Overall, these formation scenarios of BDs do not provide a successful explanation for the observed spatial distribution of BDs either.

In order to further study the nature of this finding, we need to study a larger region and confirm the membership of a larger sample of BDs. If this behavior extends to other OB stars of the cluster, we would be able to confirm whether we are indeed seeing the effect of OB stars on the formation of BDs.

We found the disk fraction of NGC 2244 to be 39 ± 9%. This value is lower than the disk fraction found in nearby SFRs of similar ages, although it is still within errors with the disk fraction of Taurus or IC 348. The low prevalence of disks in NGC 2244 was already reported in Balog et al. (2007), where the disk fraction was also found to decrease when considering only the members closest to the O-type stars. The disk depletion timescale sets the time that is available for the formation of planets, and it is therefore a fundamental parameter in order to understand the formation of planets. The low fraction of NGC 2244 members that still have a protoplanetary disk can be a consequence of faster disk dispersal due to the presence of OB stars. Overall, we found that the OB stars in NGC 2244 may have had an important role in the early evolution of the cluster, by potentially increasing the efficiency of BD formation and decreasing the lifetime of protoplanetary disks around both stars and BDs members of the cluster.

Acknowledgments

V.A-A., K.M. and K.K. acknowledge funding by the Science and Technology Foundation of Portugal (FCT), grants No. IF/00194/2015, PTDC/FIS-AST/28731/2017, UIDB/00099/2020, SFRH/BD/143433/2019 and PTDC/FIS-AST/7002/2020. R.S. acknowledges financial support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709) and financial support from national project PGC2018-095049-B-C21 (MCIU/AEI/FEDER, UE).

References

Appendix A: Gravity-sensitive indices decision boundary

We perform a derivation of decision boundaries for the three gravity-sensitive indices used for the membership classification in Sect. 3.1.4. These boundaries give for each SpT bin, the spectral index value that separates young low-gravity members from nonmembers. We use the data set presented in Almendros-Abad et al. (2022). We used the decision function given by the Support Vector Machine (SVM) method applied to the spectral index and SpT using the radial basis function kernel in order to define the decision boundary between the “young” and “field” classes of the Almendros-Abad et al. (2022) data set. SVM maximizes the distance between the different classes and a separating hyperplane defined on a feature space that is transformed using a kernel. We use the same methodology to optimize the hyper-parameters and to obtain the decision function as in Almendros-Abad et al. (2022). The decision boundaries are given in Table A.1.

Appendix B: Photospheric infrared colors

We derive the K_S-IRAC photospheric colors for a compilation of known diskless members of the Cha-I (Luhman et al. 2008; Esplin et al. 2017), Upper Sco (Esplin et al. 2018), ρ-Ophiuchus (Esplin & Luhman 2020), Taurus (Esplin & Luhman 2019) and Corona Australis (Esplin & Luhman 2022) SFRs between K0 and L0. We restrict the analysis to low-extincted sources (Av< 2 mag) and performed a running average with 1 SpT window for each 0.5 SpT (see Fig. B.1). We provide in Table B.2 the derived colors. These are the colors used in Sect. 3.2 to obtain the infrared excess in each IRAC band for the NGC 2244 candidate members.

Fig. B.1. KS - IRAC color as a function of SpT for a compilation of diskless young sources from the literature (see text) for the four IRAC channels. The black dashed line represents the mean photospheric colors for each SpT bin we derived.

We also estimate what is the typical precision and completeness of the 0.25 mag excess boundary we have used to identify sources with infrared excess in Sect. 3.2. We add all the disk-bearing sources to the compilation of diskless members of the same regions. The disk-bearing sources all have classification of the type of disk they host from the literature. We divide the disk-bearing sources in “full” (full and transitional disks) and “evolved” (debris and evolved transitional) disks. In Fig. B.2 we show the color excess CDF for the three classes disk classifications (for diskless sources we show 1-CDF for visualization purposes). In Table B.1 we provide the percentage of sources above the 0.25 mag color excess threshold for the sources with full disks, with evolved disks and without disks. We observe that in the K_S−IRAC1 and K_S−IRAC2 colors there is a significant population of full disks that are not recovered. In the four colors, using a 0.25 mag color excess for classification we have a contamination of diskless sources < 3% but we are missing most of the evolved disk sources as we expected.

Fig. B.2. CDF of the KS - IRAC colors for the three disk classifications: Full (blue dotted line), evolved (gray dashed line) and diskless (black solid line, we show 1-CDF for visualization purposes).

Appendix C: SpT - T_eff scales

In Fig. C.1 we show the comparison between the T_eff we derived in Sect. 3.4 from comparing the KMOS spectra with the BT-Settl atmospheric models (Baraffe et al. 2015) and the T_eff converted from the SpT obtained in Sect. 3.1.2 using the SpT-T_eff relationships of Luhman et al. (2003), Pecaut & Mamajek (2013) and Herczeg & Hillenbrand (2014). We observe that while the three scales agree well with our T_eff there are certain important aspects to highlight. The scale from Luhman et al. (2003), left panel) systematically produces larger T_eff than what we find by fitting the spectra with the BT-Settl models. This effect is specially important above M5.5 (∼3000 K), which could be associated with the dramatic evolution of atmospheric models since this scale was defined. The Herczeg & Hillenbrand (2014), right panel) scale has good correlation with our T_eff, but from M6 (∼2800) it also provides larger T_eff. The scale from Pecaut & Mamajek (2013), middle panel) agrees best with our T_eff values despite having been defined for field objects. In order to be as consistent as possible, we have used the SpT-T_eff relationship from Pecaut & Mamajek (2013) to obtain T_eff for the sources with only SpT available, since we have found it to perform most similarly with the T_eff we have measured by comparing the KMOS spectra with the BT-Settl atmospheric models.

Fig. C.1. Comparison of the Teff derived in this work for the KMOS spectra (see Sect. 3.4) compared with the Teff derived using the scales from Luhman et al. (2003), L03), Pecaut & Mamajek (2013), PM13) and Herczeg & Hillenbrand (2014), HH14). The black line represents the linear relationship between the variables and the black dashed lines represent a ±100 K uncertainty.

Appendix D: NGC 2244 members

In Table D.1 we provide the ID, coordinates, SpT, A_V, T_eff, whether they have KMOS or VIMOS spectrum (Obs), membership assessment from Mužić et al. (2022) and Meng et al. (2017), the infrared excess key derived in Sect. 3.2 and the final membership assessment for all the members and candidate members used in this work (see Sect. 3.3).

All Tables

All Figures

Fig. 1. Color magnitude diagram (CMD) of the core of NGC 2244. We show the Flamingos-2 catalog (Mužić et al. 2019) with black dots, objects observed with KMOS are represented by orange empty diamonds, KMOS members as solid orange circles, KMOS possible members as blue solid circles, VIMOS observed objects as empty gray diamonds, and members from VIMOS, Mužić et al. (2022) and Meng et al. (2017) as gray circles. The black solid line represents the 2 Myr isochrone shifted to the distance of 1500 pc, and the dashed line represents the same isochrone reddened by the mean extinction toward NGC 2244 (i.e., AV = 1.5 mag, Mužić et al. 2019). The isochrone shown comes from the PARSEC stellar models (Bressan et al. 2012) down to 1 M⊙, and BT-Settl below 1 M⊙ (Baraffe et al. 2015), that is isochrone 1 in Mužić et al. (2019). The right panel represents the inset of the CMD to the region populated by the KMOS observations. The black arrow represents the effect of AV = 1.5 mag.

In the text

Fig. 2. Probability map of SpT and AV from comparison with spectral templates. Darker regions represent lower values of χ2 from the comparison with the young templates for object #40. We show the best fit SpT and AV in orange. The green line represents the SpT derived from a combination of various spectral indices (see Sect. 3.1.3) for the entire AV grid. At the right and bottom panels we also show the χ2 histogram of SpT and AV.

In the text

	Fig. 3. All the KMOS spectra of objects consistent with low-mass (sub)stellar nature, together with the best fit young templates (black) de-reddened by the best fit extinction. The selected SpT and AV (in magnitude units) are also shown. The color-coding represents the final membership of the targets: orange are members, blue are possible members and gray are nonmembers.
In the text

	Fig. 4. Histogram of SpT and AV of the objects observed with KMOS derived by direct comparison with spectral templates. We also show the histograms of the members (orange) and possible members (blue).
In the text

Fig. 5. Comparison of the spectral templates best fit SpT (Sect. 3.1.2) with the SpT coming from spectral indices (Sect. 3.3). The objects are color-coded based on the final membership discussed in Sect. 3.1.4: members (orange circles), possible members (blue squares) and nonmembers (black crosses). The solid line represents the 1:1 SpT relationship, and the dashed lines represent a ±1 SpT deviation from the 1:1 relationship.

In the text

Fig. 6. Gravity-sensitive indices used in Sect. 3.1.4 to evaluate the youth of the KMOS sources compared with the SpT. Dots represent the data set from Almendros-Abad et al. (2022) color-coded according to the age class of each object (young in orange, mid-gravity in blue and field in gray). The black solid line represents the decision boundary between young and field objects (see Appendix A). The objects with KMOS spectra are color- and symbol-coded as in Fig. 5. We also show with a black arrow the effect that an extinction of AV = 5 mag would have on the values of each index.

In the text

Fig. 7. KS – IRAC color excesses for the candidate members of NGC 2244 with SpT or Teff measurement. Orange circles represent KMOS bona fide members, blue squares KMOS possible members and gray diamonds represent VIMOS objects and Mužić et al. (2022) members. Big black circles represent the objects that we found to have infrared excess. The locus of sources from the literature with full disks is represented with orange contours, and the locus of sources without disks is represented with blue contours (see Appendix B). Black dashed line represents the 0.25 mag excess boundary.

In the text

Fig. 8. Teff derived in Sect. 3.4 compared with the SpT derived in Sect. 3.1.2 for the bona fide members (orange circles) and possible members (blue squares). The young Teff scale from Luhman et al. (2003, L03, red solid line), the field Teff scale from Pecaut & Mamajek (2013, PM13, gray dashed line), and the young Teff scale from Herczeg & Hillenbrand (2014, HH14, black dotted line) are also shown.

In the text

	Fig. 9. Hertzsprung-Russell diagram of all the members of NGC 2244. The color-coding is the same as in Fig. 1. The isochrones (black dashed lines, shown for 0.5, 2, 5, 10 and 100 Myr) and the lines of constant mass (green dashed lines) shown are from Baraffe et al. (2015).
In the text

Fig. 10. IMF of the central region of NGC 2244 considering both KMOS and Meng et al. (2017) possible members with P = 0.5 (see text), represented in equal-size logarithmic bins with size Δlog(m/M⊙) = 0.15. Two different IMF representations are shown: dN/dm (top panel) and dN/dlog(m) (bottom panel). The green dashed and the pink dotted lines in the top panel represent the power-law fits performed between 0.045–0.2 M⊙ and 0.045–0.4 M⊙ respectively. The orange solid line (both panels) represents the Kroupa segmented power-law mass function (Kroupa 2001), the blue solid line (bottom panel) shows the Chabrier mass function (Chabrier 2005) and the gray dash-dotted line (bottom panel) represents a log-normal function with a characteristic mass (mc) of 0.13 M⊙ and standard deviation (σ) of 0.45. The three representations of the IMF are normalized to match the total number of objects in the cluster. The black vertical dashed line at 0.03 M⊙, and the red line at 0.045 M⊙ represent the KMOS survey limit and the true mass completeness of our sample respectively.

In the text

Fig. 11. Dependence of the star-to-BD number ratio on cluster surface density. 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 number 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 result shown for NGC 2244 is the one calculated correcting for completeness with the IMF power-law fit performed between 0.045 and 0.4 M⊙.

In the text

	Fig. 12. Cumulative distribution function of the projected distance to the closest OB star of the stars (orange line) and BDs (cyan line) members of the core of NGC 2244. We also show the distribution of nonmembers with J > 17 mag (black line) as well as the KMOS sources with very low S/N that were not considered in any part of the analysis (gray line).
In the text

	Fig. 13. Sky distribution of all the members of the core region of NGC 2244. Orange circles represent stars (< 1.2 M⊙), blue circles represent BDs and yellow stars represent the OB stars in the region coming from Preibisch et al. (2008). On the background we show the KS-band image of the center of NGC 2244 taken with Flamingos-2/Gemini-S (Mužić et al. 2019).
In the text

	Fig. B.1. KS - IRAC color as a function of SpT for a compilation of diskless young sources from the literature (see text) for the four IRAC channels. The black dashed line represents the mean photospheric colors for each SpT bin we derived.
In the text

	Fig. B.2. CDF of the KS - IRAC colors for the three disk classifications: Full (blue dotted line), evolved (gray dashed line) and diskless (black solid line, we show 1-CDF for visualization purposes).
In the text

	Fig. C.1. Comparison of the Teff derived in this work for the KMOS spectra (see Sect. 3.4) compared with the Teff derived using the scales from Luhman et al. (2003), L03), Pecaut & Mamajek (2013), PM13) and Herczeg & Hillenbrand (2014), HH14). The black line represents the linear relationship between the variables and the black dashed lines represent a ±100 K uncertainty.
In the text