© The Authors 2026

1 Introduction

The study of star formation is of crucial importance in many aspects of astrophysics, for example cosmological evolution at large scales or the formation of planetary systems at smaller scales. The prevailing view is that stars form within star clusters as a result of the gravitational collapse of molecular clouds (Kroupa et al. 2024; Hennebelle & Grudić 2024). Within these clusters, stars and brown dwarfs (BDs) of varying masses emerge, spanning several orders of magnitude, making them an ideal environment for investigating the formation of low-mass stars and BDs.

The initial mass function (IMF) describes the distribution of stellar masses at birth and serves as a fundamental tool for comparing the outcome of star formation in various environments. Above 1 M_⊙, the IMF is often considered universal and follows the Salpeter power-law slope (dN/dM ∝ M^−α with α = 2.35; Salpeter 1955), although a slightly shallower slope has been observed in some regions (see Bastian et al. 2010 for a comprehensive review). To describe the sub-solar range, Kroupa (2001) proposed a broken power law with α = 0.3 in the 0.01–0.08 M_⊙ range and α = 1.3 in the 0.08–0.5 M_⊙ range. Additionally, Chabrier (2003, 2005) provides an alternative description, suggesting a log-normal function ($(dN/d\log m\propto e^{-\frac{{\left(\log m-\log m_c\right)}^2}{2{\sigma }^2}})$) with a peak at 0.25 M_⊙.

Following the power-law description, observational studies of the low-mass IMF report a range of values for α: α = 0.67 ± 0.23 and α = 0.29 ± 0.22 for Cha-I and Lupus 3, respectively, for the mass ranges 0.02–0.2 M_⊙ and 0.04–0.2 M_⊙ (Mužić et al. 2025), α = 0.7–1.1 for NGC 2244 for 0.045–0.2 M_⊙ or 0.045–0.4 M_⊙ (Almendros-Abad et al. 2023), α = 0.51 ± 0.15 for RCW 38 for masses 0.02–0.2 M_⊙ (Mužić et al. 2025), α = 0.18 ± 0.19 for masses 0.004–0.19 M_⊙ for σ Orionis (Damian et al. 2023), and α = 0.26 ± 0.04 for masses 0.011–0.4 M_⊙ for 25 Ori (Suárez et al. 2019). These differences in values could arise from uncertainties affecting the IMF derivation rather than intrinsic variations between clusters. Key factors influencing the determination of the IMF include assumptions about stellar age, cluster distance, theoretical models, membership selection methods, extinction laws, and varying mass ranges for the fits presented here (Scholz et al. 2013).

The low-mass IMF is closely linked to the BD production efficiency. The typical ratio of low-mass stars (0.075–1 M_⊙) to BDs (0.03–0.075 M_⊙) ranges from 2 to 6 (Mužić et al. 2015, 2017; Andersen et al. 2008; Almendros-Abad et al. 2023). Several mechanisms have been proposed for the formation of BDs, including turbulent fragmentation (Padoan & Nordlund 2002), ejection from multiple systems (Reipurth & Clarke 2001), formation in massive circumstellar discs (Stamatellos & Whitworth 2009), and the photoevaporation of pre-stellar cores by OB-type stars (Whitworth & Zinnecker 2004). However, studies of nearby regions suggest that most massive BDs likely form in a manner similar to stars (see Luhman 2012 for a review), implying the necessity of examining more massive clusters to probe potential environmental variations of the IMF. Several of these theories postulate that the BD formation efficiency depends on environmental factors such as stellar density, gas density, and the presence of massive stars (Luhman 2012). However, this remains an active area of research since it is not yet clear to which extent each of these factors affects the low-mass IMF.

One of the most massive star-forming regions in the Galaxy is the Carina Nebula molecular complex, which contains young massive clusters: Trumpler 14 (Tr 14), Trumpler 15, and Trum-pler 16 (Tapia et al. 2003). Tr 14 contains ∼20 OB-type stars (Shull et al. 2021), which could influence the IMF (Whitworth & Zinnecker 2004). Moreover, we estimate the stellar density of Tr 14 to be ∼10³ stars/pc² (see Sect. 4.4), which is several times denser than the Orion Nebula cluster as found in Almendros-Abad et al. (2023). Estimates of Tr 14’s age vary between <1 Myr and 5 Myr (Tapia et al. 2003; Ascenso et al. 2007; Hur et al. 2012; Itrich et al. 2024), with evidence of two distinct populations: a younger core (<1 Myr) and an older halo (∼5 Myr; Ascenso et al. 2007). The latest spectroscopic analysis, based on MUSE integral field unit observations, places the age of Tr 14 at ∼ 1 Myr (Itrich et al. 2024), which we adopted in our study. The cluster distance was taken as 2.37 ± 0.15 kpc, based on Gaia Early Data Release 3 (Shull et al. 2021), which is in agreement with the results of Göppl & Preibisch (2022), which are based on Gaia Data Release 3 (Gaia Collaboration 2023, Gaia Collaboration 2016). While Ascenso et al. (2007) found a normal extinction law (R_V = 3.1) with spatial variations, Hur et al. (2012) suggested an anomalous extinction law towards Carina. The extinction law that we use throughout the paper is from Wang & Chen (2019): R_V = 3.1.

Previous studies (Ascenso et al. 2007; Rochau et al. 2011; Hur et al. 2012) of the IMF in Tr 14 have primarily focused on the stellar content, without robust coverage of the low-mass and substellar populations. For example, Ascenso et al. (2007) detect sources down to K_S ∼18.5 mag, but they are only 90% complete to K_S ∼14.5 mag, corresponding to ∼1–2 M_⊙ at 1 Myr depending on extinction and evolutionary models. Hur et al. (2012) report an observational limit of ∼1.6 M_⊙. In contrast, Rochau et al. (2011) reach significantly deeper, down to ∼0.25 M_⊙. With its rich massive-star population and young age, Tr 14 provides a valuable setting to investigate BD formation and the effect of environmental factors in shaping the low-mass end of the IMF.

In this work we present the most in-depth substellar survey of Tr 14 to date. We used high-resolution near-infrared imaging from the Gemini-South Adaptive Optics Imager. Our analysis focuses on the central ∼1.4′ × 1.4′ region of the cluster, corresponding to ∼0.98 × 0.98 pc² for a distance of 2.4 kpc. We aimed to characterise the low-mass IMF and evaluate the efficiency of BD formation in a high-radiation environment. We present our observations, the data reduction, and the photometry in Sect. 2. The statistical determination of cluster membership is described in Sect. 3. The derivation of individual object masses and the IMF is presented in Sect. 4. A summary of our results and conclusions are given in Sect. 5.

2 Observations, data reduction, and photometry

2.1 Observations

We observed the Tr 14 cluster in the JHK_s filters using the Gemini-South Adaptive Optics Imager (GSAOI; McGregor et al. 2004; Carrasco et al. 2012) assisted by the Gemini Multi-Conjugate Adaptive Optics System (GeMS; Rigaut et al. 2014; Neichel et al. 2014), as a part of programme ID GS-2019A-DD-107 (PI: M. Andersen). Each exposure was 60s, and a total of 25, 15, and 15 exposures were obtained, in J, H, and K_S, respectively. Additionally, short exposures of 10s each were taken, with 10 exposures per filter. Three natural guide stars were utilised for the adaptive optics correction; their coordinates are given in Appendix A. All observations were obtained under clear sky conditions. The image quality in the final images for J, H, and K_S filters, is 0.14″, 0.13″, and 0.11″, respectively.

2.2 Data processing

The images were processed using standard procedures with the THELI pipeline (Schirmer 2013), and the GSAOI-specific adaptations in Schirmer et al. (2015). Pre-processing included dark subtraction, flat-fielding with domeflats, and the construction of background models using median-stacking of the dithered science images after source masking.

The astrometry and distortion correction were based on the Gaia Early Data Release 3 reference catalogue. The final co-added image was resampled to the nominal GSAOI plate scale of 0.0197″ per pixel. We show the J-band image in Fig. 1.

2.3 PSF fitting photometry

We used SOURCEEXTRACTOR (Bertin & Arnouts 1996) and PSFEx (Bertin 2013) to obtain the point-spread-function (PSF) model and then extracted sources and measured the photometry using PSF-fitting. We tested various parameters to identify the optimal PSF fitting, including BASIS_NUMBER, PSF_SIZE, PSFVAR_KEYS, PSFVAR_GROUPS, PSFVAR_DEGREES, SAMPLE_MINSN. PSFVAR_KEYS and PSFVAR_GROUPS represent the image parameters with which the PSF varies (e.g. X and Y pixel position in the image) and the polynomial group, respectively. PSFVAR_DEGREES stands for the degree of the polynomial of the PSF. We used the same parameters consistently on all the images; the values are given in Table B.1.

Since our field is crowded, for each band, both in long-and short-exposure images, we made two separate extractions, following the concept of ‘hot’ and ‘cold’ detection as outlined in Rix et al. (2004). With just one extraction we were not able to extract both the bright and the faint sources at the same time. Thus, the first extraction (called extraction 1 hereafter) parameters were optimised to gain more sources that are brighter and easier to resolve, and the second extraction (or extraction 2) was done to obtain fainter sources in crowded areas. The most significant parameters for extraction (BACK_SIZE, BACK_FILTERSIZE, DETECT_THRESH, and DEBLEND_MINCONT) affected which sources were detected in the images. The BACK_SIZE parameter defines a cell size in which the local background is estimated using the mode value, while the BACK_FILTERSIZE parameter smooths the background map by applying a median filter. The parameter DETECT_THRESH dictates the threshold above the sky background in units of standard deviation of the background noise. The DEBLEND_MINCONT parameter dictates the contrast for de-blending multiple sources.

The obtained catalogues were cleaned for PSF_MAG=99, FLAGS>=4, and SPREAD_MODEL values deviating by more than 3σ from the mean (see Table B.1 for the exact values). PSF_MAG=99 represents a bad PSF magnitude extraction, for example a saturated source, while SPREAD_MODEL excludes all the deviated and ellipsoidal sources. FLAGS>=4 rejects sources on the edges of the images and with saturation, or corrupted photometric aperture (see SOURCE-EXTRACTOR documentation for more details; Bertin & Arnouts 1996).

2.4 Photometric calibration

Once we obtained the photometry (Sect. 2.3), the extraction 1 and 2 catalogues were combined and then calibrated using the catalogue from Ascenso et al. (2007). The calibration catalogue contains data from two different instruments, with different spatial resolutions. We used only the low-spatial resolution instrument, to avoid possible errors (e.g. matching two stars from our catalogue to one star in the calibration catalogue) in the cross-matching process, and matched it with our catalogue with a tolerance of 0.3″. Moreover, for calibration purposes, we excluded the sources in our catalogue that have another source within 0.5″, to avoid mismatches between the two catalogues. The calibration catalogue from Ascenso et al. (2007) has been calibrated against the Two Micron All Sky Survey (2MASS) photometric system (Skrutskie et al. 2006). We calibrated the photometry using the zero-points (ZPs) only and neglected the colour terms due to them being consistent with zero within the errors. To derive the ZP values for photometric calibration, we took the difference of the instrumental magnitudes from our data and reference magnitudes from the calibrated catalogue, taking measurement errors into account and using 3-sigma clipping. The final ZP was calculated as a weighted average of the surviving data points, with weights determined by the inverse squared errors of both datasets, providing the ZP and its associated uncertainty. As a result, for the long exposure, we get ZP_J = 25.97 ± 0.11 mag, ZP_H = 26.18 ± 0.15 mag, ${\text{ZP}}_{{\text{K}}_s}$ = 25.64 ± 0.16 mag, and for the short exposure ZP_J = 25.99 ± 0.11 mag, ZP_H = 26.33 ± 0.11 mag, ${\text{ZP}}_{{\text{K}}_s}$ = 25.75 ± 0.10 mag. The photometric uncertainties in individual J, H, and K_S catalogues were calculated by combining the uncertainties of the ZPs and the measurement uncertainties supplied by S

ource

-E

xtractor

, and are shown in Fig. 2.

The final catalogue for each band was obtained by matching the long and short-exposure catalogues. For the J-band, we used 3-pixel separation, while the rest had 2.5-pixel separation due to the smaller full-width-at-half-maximum of the sources. For the sources that appear in both catalogues, we kept the values with the lower magnitude error. Finally, we cross-matched the catalogues in the three photometric bands using 2.5-pixel separation. We cut 400 pixels off each edge to exclude sources with lower signal-to-noise ratio. The final catalogue contains 1207 sources with photometry in all three photometric bands, 3 additional sources in only J and H bands, 12 sources in only J and K_S bands, and 356 sources in only H and K_S bands. The first three lines of the catalogue are given in the Table 1 and the full catalogue is provided electronically.

2.5 Completeness

We evaluated the photometric detection completeness through artificial star experiments. Artificial stars were randomly positioned in our images using the source function of Artificial Stellar Populations (A

RT

P

OP

) Python package (Greco & Danieli 2022). We injected the stars using the PSF model obtained with PSFEx and the parameters listed in Sect. 2.3. We inserted 50 stars at a time and repeated it 10 times, to improve statistics, in steps of 0.2 mag. The extraction used the same parameter values, PSF model, and procedure as in the photometric calibration. Figure 3 shows the fraction of recovered stars in comparison to the inserted artificial ones, as a function of magnitude. The 20% completeness limit is at J = 21.3 mag, H = 20.9 mag and K_S = 20.3 mag. This corresponds to the mass of ∼0.01 M_⊙, assuming the distance of 2.37 kpc, 1 Myr for the age and combined PARSEC and BT-Settl model (see Sect. 4.1), and with the derived median extinction of the cluster, A_V = 2.2 mag (see Sect. 4.1). The 90% completeness limit is at J = 18.4 mag, H = 17.9 mag and K_S = 17.7 mag, which corresponds to masses of 0.06–0.08 M_⊙. Table 2 summarises our results for the 20%, 50% and 90% limits. The previously deepest data in Tr 14 go down to ∼0.065M_⊙ (Itrich et al. 2024) at 30% completeness, not including the core of the cluster.

3 Cluster membership

To examine the cluster population, we show two colour-magnitude diagrams (CMDs) of the observations in the first two panels of Fig. 4. The cluster sequence can be fairly well distinguished, with colours J−K_S ≈ 1–1.7 mag, but there is a sparsely populated foreground population (to the left of the cluster sequence), and background contamination (to the right) not following the shape of the isochrone and located at the fainter end of the CMD, both of which probably partially overlap with the cluster population. The background population seen in the CMDs in the cluster direction appears to be significantly redder than the cluster sequence. The cluster seems to be located in front of a molecular cloud as seen from the Earth (Ascenso et al. 2007) and it is therefore expected that the background sources suffer larger extinction.

Therefore, in our field of view, we see both the cluster members and field stars contaminating our dataset. We had to be able to differentiate between the two to do a proper analysis of the cluster. Thus, we needed to perform a membership determination. In this section we describe the decontamination of the observations performed using two independent control field (CF) datasets: the Besançon Galaxy model¹ and a CF obtained using VIRCAM data from the Visible and Infrared Survey Telescope for Astronomy (VISTA) telescope (Preibisch et al. 2014). The setups for the CFs are explained in Sects. 3.1 and 3.2, respectively, followed by the common decontamination procedure in Sect. 3.3.

3.1 Decontamination setup with a Galaxy model

To address the contamination with the Besançon Galaxy model presented in Czekaj et al. (2014), as an input we used the central position (l = 287.4107^◦, b = −00.5794^◦) of our cluster with an area 24.5 times larger than our cluster area on the sky, which was scaled later in the process of decontamination (Sect. 3.3) with the photometric errors given as in Appendix C. We took this size of the area to get a significant number of the expected field stars. The output of the obtained Besançon Galaxy model contains the J-band magnitude, J–H, J–K, and H–K colours, age, mass, coordinates in right ascension and declination, distance and extinction. We transformed the model magnitudes from the Koornneef system to 2MASS using the relations given in Carpenter (2001).

Since the Besançon Galaxy model provides the distances of the simulated stars, we divided the model output into foreground and background populations by comparing those values with the cluster’s location at 2.37 kpc. The foreground population was used for decontamination as is, i.e. without applying any extinction. For the background population, we estimated the values of extinction by reddening the CF photometry until it coincided with the background contaminants seen in the cluster CMD (see Fig. 5). We thus decontaminated the cluster from the background after reddening it with the extinction values of A_V = 5, 6, and 7 mag. The decontamination procedure, following the setup described here, is detailed in Sect. 3.3.

3.2 Decontamination setup with VISTA

With the 4-m VISTA of ESO, Preibisch et al. (2014) mapped an area of 6.7 deg² around the Carina Nebula in the near-infrared J-, H-, and K_S -bands, a portion of which we used as a CF. To address the contamination with the CF, we took the shape of a rectangle centred at 10h40m48.00s RA, –59^◦31′12.00″ Dec, and 24.5 times bigger than the size of the cluster area. The CF is ∼ 0.4^◦ away from our cluster and they have a similar galactic latitude (b_CF = −1.118^◦, b_Tr14 = −0.575^◦). We took this position and size to get a significant number of field stars. To separate the foreground and background populations of the CF, we used Gaia parallaxes.

We cross-matched our catalogue with Gaia Data Release 3 (Gaia Collaboration 2016, 2023; Babusiaux et al. 2023) to obtain the proper motions and parallaxes, and we performed two cuts. Firstly, we made cuts in proper motion in right ascension and declination, keeping only the most probable members of the cluster, which should behave similarly in the proper motion diagram (PM_RA = −6.49±0.60 mas yr⁻¹, PM_Dec = 2.14±0.58 mas yr⁻¹). We did this by applying sigma clipping using a threshold of 3σ with a maximum of five iterations to iteratively exclude outliers that deviate significantly from the central distribution. This ensures that the derived statistical parameters reflect the intrinsic kinematic properties of the population, minimising the impact of contaminants or measurement errors. Secondly, we applied an additional brightness cut on the subset obtained from the first selection, keeping only the brightest (<16 mag) objects to calculate the weighted mean parallax and parallax error. This value was then used to identify and reject foreground objects in the Tr 14 catalogue based on a parallax zeta score >2. We calculated the zeta score with the expression $\zeta =\left|\omega -\overline{\omega }\right|/\sqrt{{\sigma }_{\omega }^2+{\sigma }^2}$, where ω and σ_ω are the parallax and the uncertainty for each star, while the $\overline{\omega }$ and σ are the weighted mean and uncertainty of the whole sample. We were left with a dataset of our cluster that contains probable cluster members as well as the background contaminants.

We then made the same cut based on the parallax zeta score in the VISTA CF catalogue, which we previously cross-matched with Gaia Data Release 3, to keep only the subset that corresponds to a probable background. We show the results of these cuts in Fig. 6. To this VISTA CF catalogue with removed foreground objects, we applied different amounts of extinction (A_V = 6.5, 7, and 7.5 mag; see Fig. 5), which we then used in the procedure explained in Sect. 3.3.

3.3 Statistical determination of membership

To statistically determine the cluster membership, after the independent setups for the CFs in Sects. 3.1 and 3.2, we applied the procedure from Mužić et al. (2019), which is summarised as follows:

1. The CMD containing the cluster and the CF, as in Fig. 5 (upper panels if we use the Besançon model as the CF and lower panels if we use VISTA as the CF), is divided in a grid with cell sizes ∆colour and ∆magnitude.

2. For each cell, the expected number densities of stars in the cluster field of view and in the CF (either from the Besançon model or from VISTA) are calculated. The densities are then scaled for different on-sky areas of the cluster and the CF.

3. he expected number of field objects in each cell is then randomly removed from the cluster population.

Each object in the CMD is represented with a Gaussian probability distribution in the space of colour and magnitude, with the width being determined by the photometric uncertainties. The part of the Gaussian that lies within the cell is then integrated and summed to get the expected number of objects for each cell. Cell sizes are determined in a way to have good statistics of cluster members and CF objects in each cell. We then removed from the cluster field as many sources as computed in the CF in each cell. The specific objects to be removed are chosen randomly. We used different cell sizes with ∆magnitude = (0.4, 0.5, 0.6) mag, which results in N_mag = 3 different magnitude bin sizes and ∆colour = (0.3, 0.4) mag, yielding N_col = 2 colour bin sizes. Additionally, to improve accuracy for the objects laying at the edges of our cells, for each magnitude-colour combination, we shifted the cell centre for 0.0, 1/3, and –1/3 of the cell size in each dimension. This created N_s = 3 · 3 = 9 cell centre configurations. We repeated the decontamination procedure for each combination of various parameters, namely cell sizes and cell centres, as well as the applied reddening. We looked at N_cmd = 2 CMD types (J, J–K,J−H) for $N_{A_V}$ = 3 extinction values as shown in Fig. 5. This results in N_mag · N_col · N_s · N_cmd · $N_{A_V}$ = N_L = 324 lists of potential cluster members, which contain IDs of catalogue sources that we later employed in the procedure for the IMF determination of the cluster. We calculated the probability of each star in the catalogue being a cluster member by counting the number of times it appears in all lists and then dividing by the total number of lists. This methodology follows the similar approach previously applied in other young clusters such as NGC 2244 (e.g. Mužić et al. 2019), and has been validated through spectroscopic analysis by Almendros-Abad et al. (2023), demonstrating its robustness in separating cluster and field populations. The resulting probabilities are shown in Fig. 7, where we use the VISTA (bottom panel) and the Besançon model (upper panel) as CFs in the decontamination method. In the bottom panel, the foreground population appears as high-probability members since for the decontamination with VISTA we cross-matched our catalogue with Gaia, which is much shallower than our photometry, and thus the sample is limited. We then removed the foreground objects, which are known to contribute to contamination based on their Gaia parallaxes, as explained in Sect. 3.2.

4 Results and discussion

4.1 Extinction and masses

To derive the IMF of our cluster, we had to estimate the masses of the cluster members. We used evolutionary models to get these values. To obtain extinction and masses, we de-reddened our sources to the J, J–H isochrone in the CMD, as in Mužić et al. (2017). The full CMD and colour-colour diagram (CCD) are shown in Fig. 4. We combined the PARSEC and BT-Settl models at 1 M_⊙ to a 1 Myr isochrone (Bressan et al. 2012; Chen et al. 2014; Baraffe et al. 2015) as in Mužić et al. (2019).

First, we generated 1000 magnitude values for each source assuming a normal distribution, where the standard deviation is the photometric uncertainty. We checked whether the source with the uncertainties was inside or crossing the model limits (solid red and dotted red lines in Fig. 4 in the CCD); otherwise, the mass was not derived. The extinction was then derived by de-reddening in the CMD to the 1 Myr isochrone in these 1000 runs. The mass was obtained by interpolating the mass-magnitude relation corresponding to the model isochrone.

An example of the mass and extinction derivation procedure of one of the sources is shown in Fig. 8. The right extinction line delimits we would expect to see T Tauri (Meyer et al. 1997) and Herbig AeBe (Hernández et al. 2005) stars (we would expect to see them to the left of this line). However, as argued in Mužić et al. (2019), the excess coming from the circumstellar disc of the classical T Tauri or Herbig AeBe stars appears to have a negligible impact on the IMF; thus, we did not attempt to correct for it. Additionally, mass and extinction histograms are plotted. In these runs, it is possible to obtain negative A_V; however, this is not physical. To correct for it, without skewing the distributions, only data points with A_V >= 0 are used to estimate the IMF.

As a result of this procedure, we get the mass distributions, median mass, median and mean extinction with its standard deviation for each source in the catalogue. The mass distributions are later used to obtain the IMF.

4.2 Initial mass function

We aimed to derive the IMF from a set of mass distributions, obtained as described in Sect. 4.1, using the lists containing IDs of catalogue sources that remain after the decontamination procedure described in Sect. 3.

Each star’s mass is represented by a mass distribution that is smoothed using a Gaussian kernel density estimation. From the mass distribution of each star, we drew a sample of N₁ = 100 values per object; this resulted in 100 different mass distributions for the entire sample. Next, for each of the N₁ resampled sets, bootstrap (random samplings with replacement) is performed N₂ = 100 times to account for statistical fluctuations. This results in N₁ · N₂ individual mass distributions. During each iteration, we evaluated the IMF as a histogram with variable bin sizes, chosen to contain roughly the same number of objects and save it. The same bins were used in each iteration. We corrected this mass function for completeness by dividing each star in the bin individually with its corresponding completeness value. We excluded all the sources with completeness <= 20% from the calculation of the IMFs because we did not want to be biased by the faint objects that have low completeness. To ensure the robustness of the derived IMF, we then averaged the N₁ · N₂ · N_L IMFs obtained and calculated the standard deviation to estimate the uncertainty of each bin.

The results for the IMFs are shown in Fig. 9 and summarised in Table 3. The obtained slope α (dN/dM ∝ M^−α) when fitting between 0.01 and 0.2 M_⊙ differs depending on whether we do the decontamination with the Besançon model or the VISTA CF: using the Besançon model, we get α = 0.03 ± 0.24, whereas using the VISTA CF, we obtain α = 0.14 ± 0.19. When fitting the IMF between 0.2 and 4.5 M_⊙, we find α = 1.72 ± 0.04 for both decontamination methods.

We note that the decontamination primarily affects the lowest-mass bin, whereas above this limit the results remain very similar. On the low-mass part of the IMF, the Besançon model may be overestimating the number of contaminants. This could originate from Tr 14 being a cluster located in a young star-forming region, implicating the presence of nebulosity and extinction. The approach using VISTA as the CF is more realistic because it directly samples the actual background and environmental conditions, whereas the Besançon model relies on assumptions that may not fully take into account local conditions. In particular, the model does not include explicit spiral-arm overdensities along this line of sight, such as the Carina arm, which could also cause it to underestimate the number of field stars near the cluster distance. However, even if such an over-density were included, it would likely be uniformly distributed across the CMD and thus would not affect our decontamination result, especially in the lowest-mass bin where the largest difference relative to the VISTA CF is observed. Meingast et al. (2016) find that the Besançon model overestimates number of faint (K_S > 19 mag) and bright (K_S < 10 mag) objects, which could explain the results that we find. However, our approach with the VISTA CF could, due to its position, sample different parts of the Milky Way and thus potentially not see general Carina contamination, in addition to underestimating the foreground contamination, since Gaia does not have any parallax measurement for faint objects as explained in Sect. 3.2. Overall, the true number of contaminants, and consequently the resulting IMF slope, likely lies between the values inferred from the VISTA and Besançon model decontamination approaches, with the two CFs yielding consistent results when the lowest-mass bin is excluded.

In the dN/dlogM representation of the IMF (Chabrier 2003), as shown in Fig. 10, we fitted a log-normal function. We derived a characteristic mass m_c = (0.20 ± 0.02) M_⊙ and σ = 0.45 ± 0.03 with the decontamination approach using the Besançon model, while for VISTA we get m_c = (0.20 ± 0.02) M_⊙ and σ = 0.50 ± 0.03. Rochau et al. (2011) studied the central ∼1 × 1 arcmin² of Tr 14. They report a higher characteristic mass of ∼0.52 M_⊙. For the dN/dM power law, they find α = 0.37±0.32 for masses 0.25–0.5 M_⊙ and α = 1.52 ± 0.3 for masses 0.5–3.2 M_⊙. In contrast to Rochau et al. (2011), we do not find statistically significant evidence of a break in the IMF slope around 0.5 M_⊙. This may be due to our extended mass coverage and improved completeness, especially below 0.5 M_⊙. Compared to Damian et al. (2021), who find m_c = (0.32 ± 0.02) M_⊙ and σ = 0.47 ± 0.02 from a sample of nine clusters, our characteristic mass is lower, while the width of the distribution (σ) is similar or slightly broader. This shift towards lower m_c may reflect differences in environmental conditions specific to Tr 14, or it may result from our deeper completeness and sensitivity to lower masses, particularly in the BD regime.

The results we obtain for the high mass part are similar to other studies in this mass range: for RCW 38, α = 1.48 ± 0.08 for masses 0.2–20 M_⊙ (Mužić et al. 2017); for σ Ori, α = 1.73±0.16 for masses 0.35–19 M_⊙ (Peña Ramírez et al. 2012); but shallower than: 25 Ori where α = 2.50 ± 0.11 for masses 0.4–13.1 M_⊙ (Suárez et al. 2019), NGC 2244 where α = 2.3 ± 0.3 for masses 1.5–20 M_⊙ (Mužić et al. 2022) and for the standard Salpeter slope (α = 2.35). Our derived substellar IMF slope is significantly shallower than the typical values reported in the literature. For example, α = 0.67 ± 0.23 in Cha-I for masses 0.02–0.2 M_⊙ and α = 0.51 ± 0.15 in RCW 38 for masses 0.02–0.2 M_⊙ (Mužić et al. 2025); α = 0.7 ± 0.12 for masses 0.045–0.2 M_⊙ in NGC 2244 (Almendros-Abad et al. 2023). Compared to these regions, our results suggest a flatter distribution of substellar objects, implying a lower relative BDs with respect to stars. However, we find α = 0.33 ± 0.22 and α = 0.26 ± 0.16 for the VISTA and Besançon model CF in the 0.03–0.2 M_⊙ mass range, similar to α = 0.29 ± 0.22 in Lupus 3 for masses 0.04–0.2 M_⊙ (Mužić et al. 2025), α = 0.26 ± 0.04 for masses 0.011–0.4 M_⊙ in 25 Ori (Suárez et al. 2019) and α = 0.18±0.19 for masses 0.004–0.19 M_⊙ in σ Ori (Damian et al. 2023). While our lowest-mass range 0.01–0.2 M_⊙ yields particularly shallow slopes, this result is based on the least massive point, where the incompleteness is highest and uncertainties are largest. If this point is excluded, the remaining slope values are consistent with those reported for other young clusters. The results are illustrated in Fig. 11 (see e.g. Hennebelle & Grudić 2024 for a more comprehensive summary).

4.3 Testing IMF assumptions

4.3.1 Dynamical evolution

To interpret the observed mass function in our field of view as IMF, it is necessary to assume that no significant fraction of cluster members resides outside our observed region, and that dynamical evolution has not substantially altered the mass distribution since formation. While such assumptions are generally justified for very young clusters, Tr 14 is both massive and compact, with a half-mass radius of 1.5–1.8 pc and a total stellar mass of several thousand solar masses (Shull et al. 2021). Its crossing time is estimated at ∼0.7–0.8 Myr, and its two-body relaxation time ranges from 40 to 80 Myr (Shull et al. 2021). Given the cluster’s age of ∼1 Myr, this places it just beyond one crossing time but still far from being dynamically relaxed. This suggests that dynamical effects such as the ejection of members or significant radial redistribution (especially for low-mass objects) are unlikely to have strongly influenced our sample. While Sana et al. (2010) report signs of mass segregation, most likely due to dynamical friction, among stars more massive than ∼10 M_⊙, Shull et al. (2021) argue that even the most massive stars (30–80 M_⊙) are unlikely to have undergone significant segregation within the cluster’s 1–3 Myr lifetime. While this supports the idea that our sample, limited to lower-mass stars and BDs in the central region, is likely to still reflect the IMF, some degree of primordial mass segregation cannot be ruled out. A wider field study would be needed to verify whether spatially dependent effects like mass segregation influence the mass function at larger radii.

4.3.2 Crowding

To assess the potential influence of crowding and the halos of bright stars on the derived IMF, we compared the IMFs extracted from two distinct regions within our field of view. We defined the most crowded region by drawing a rectangle over the central area of the image, while the rest of the field encompassing less crowded surroundings was treated as a separate region for comparison. We computed the photometric completeness separately for each region and subsequently derived their respective IMFs. We found only small differences between the two IMFs. For the mass range 0.01–0.2 M_⊙, we measured slopes of α_VISTA,inside = 0.07 ± 0.18 and similarly α_{VISTA,outside} = 0.16 ± 0.22, while for the Besançon model CF α_model,inside = −0.07 ± 0.22 and α_{model,outside} = 0.05 ± 0.26. In the mass range 0.2–4.5 M_⊙, we find α_VISTA,inside = 1.68 ± 0.06 and α_{VISTA,outside} = 1.76 ± 0.04, while for the Besançon model CF α_model,inside = 1.65 ± 0.05 and α_{model,outside} = 1.77 ± 0.05. These results indicate that crowding and the presence of bright stellar halos do not significantly bias our IMF determination.

4.3.3 Evolutionary models

To be able to cover the full range from the substellar to high stellar masses, in our analysis we used a BT-Settl and PARSEC isochrones combined into a single isochrone at 1 M_⊙. To estimate the impact of this combination, we evaluated the IMF with the two isochrones separately, within the common mass range. The results are given in Table 4. We see that with the PARSEC (Bressan et al. 2012) evolutionary model we have an overabundance of low-mass sources in comparison to the results obtained using BT-Settl evolutionary models (Baraffe et al. 2015). Using the PARSEC isochrone, in the process of deriving masses, a large portion of objects will be crossing the lowest point on the isochrone, creating the artificial abundance of low-mass objects. We thus conclude that using a combined isochrone does not produce any anomalous effects at the intersection of the isochrones, and that the IMF is consistent with both models.

4.4 Star-to-BD ratio

To calculate the star-to-BD ratio, we used the same method as in the IMF derivation, up to the point where we binned the data. At this point, we instead counted the number of sources in the mass ranges 0.03–0.075 M_⊙ for the BDs and 0.075–1 M_⊙ for the stars and correct for completeness, restricting only to sources with completeness ≥20%. The median star-to-BD that we obtain using the VISTA CF is 4.0, with a 95% confidence interval of 2.8–5.8, while for the Besançon model it is 4.3 with the 95% confidence interval of 3.1–6.0. As stated previously, we consider the result with the VISTA CF to be more realistic. This result is in agreement with the star-to-BD ratios in RCW 38 (2.0±0.6), Cha-I (3.2–4.8), and Lupus 3 (2.1–4.5) from Mužić et al. (2025). Our result is in agreement with Andersen et al. (2008), who report star-to-BD ratios in different regions in the range 3.3–8.5.

To place Tr 14 in context of different environments, we estimated its surface density and far-UV (FUV) flux from massive stars. We calculated a mean surface density of ∼1000 stars/pc² by dividing the number of sources in our field of view by the area. To estimate the FUV radiation, we followed the process outlined in Mužić et al. (2025). We derive a flux of approximately 10^5.54 G₀.

The Fig. 11 in Almendros-Abad et al. (2023) shows that the star-to-BD ratio does not significantly vary with stellar surface density, nor with the presence of OB stars. Their analysis includes clusters with surface densities in ∼10–1000 stars pc⁻² range, including the Orion Nebula cluster, RCW 38, NGC 2244, Lupus 3, Cha-I, NGC 1333 and IC 348. Tr 14 is on the high end of this range, placing it among the densest clusters examined and only surpassed by RCW 38.

In terms of FUV radiation, Fig. 8 of Mužić et al. (2025) presents the median FUV flux for Lupus 3, Cha-I, NGC 2244, CrA, NGC 1333, and RCW 38. Tr 14, with its estimated flux of 10^5.54 G₀, exceeds all the clusters in their sample, including RCW 38 (∼10^5.3 G₀).

Despite this extreme environment, high density and high FUV flux, our measured star-to-BD ratio evaluated for the mass range 0.03–0.075 M_⊙, lies between ∼3 and 6. This range is comparable to those observed in both low- and high-density environments in Almendros-Abad et al. (2023) and Mužić et al. (2025), suggesting that the star-to-BD ratio in Tr 14 does not significantly differ from other regions in this mass range. However, we note that the IMF we derived is significantly flatter than those in the comparison regions. This apparent contradiction can be explained when considering the slope uncertainties and excluding the lowest-mass bin, which drives much of the flattening; with that bin removed, the slope becomes more consistent with other clusters. However, if we were to include lower-mass BDs (<0.03 M_⊙), where our IMF shows a clear deficit, we would likely start to see deviations from the trend, indicating a possible environmental influence on the formation efficiency of the lowest-mass substellar objects. Interestingly, while Almendros-Abad et al. (2023) propose that OB stars may enhance BD formation in NGC 2244, our results from the denser and more UV-intense environment of Tr 14 do not support this scenario, pointing instead to a potential suppression. It may be that a combination of high density and intense radiation fields inhibits the fragmentation processes needed for efficient formation of very low-mass objects.

5 Summary and conclusions

Using GSAOI/GeMS observations, we have presented the most complete mass function in Tr 14 to date, covering objects from the substellar regime up to 4.5 M_⊙. Our final catalogue contains 1207 sources in the J, H, and K_S bands. The data are 50% complete down to 0.02 M_⊙ and 20% complete down to 0.01 M_⊙. Because there are both field stars and cluster members in our dataset, we statistically decontaminated the catalogue using either the Besançon model or observational data from the VISTA telescope as a CF. Masses and extinctions were obtained from comparisons to the evolutionary models from Bressan et al. (2012), Chen et al. (2014), and Baraffe et al. (2015). The resulting IMF below 0.2 M_⊙ differs depending on whether we use the Besançon model or VISTA for the decontamination process. We believe that the Besançon model overestimates the contamination and therefore prefer the results obtained using the observational CF.

Using the VISTA CF, we find the slope of the high-mass part of the IMF (∼0.2–4.5 M_⊙) to be α = 1.72 ± 0.04, which is similar to other young star-forming regions but shallower than the standard Salpeter slope. For the low-mass end of the IMF (∼0.01–0.2 M_⊙), we find a relatively shallow slope of α = 0.14 ± 0.19. While this appears significantly flatter than values reported for other regions, the comparison should be treated with caution: not only do we probe a slightly different mass range than most studies, but the lowest-mass bin in this mass range is also the most affected by incompleteness, leading to larger uncertainties. If the least massive point is excluded, the remaining slope values are consistent with those reported for other young clusters.

Nonetheless, the observed flattening could suggest that the number of BDs declines steeply with decreasing mass, especially below 0.03 M_⊙, in turn suggesting a lower formation efficiency for substellar objects in this environment. The cause of this suppression is unclear but may relate to the specific environmental conditions of Tr 14. Its mean surface density of ∼1000 stars/pc⁻² lies between that of NGC 2244 (lower) and RCW 38 (higher), yet the substellar IMF we find is much flatter than in either of those regions. In terms of FUV flux, Tr 14 is exposed to stronger radiation fields than both NGC 2244 and RCW 38, and this may inhibit the formation or survival of its lowest-mass BDs. Though we cannot yet conclude this definitively, our findings highlight the importance of considering both stellar density and radiative feedback in shaping the low-mass end of the IMF.

Finally, we find that the star-to-BD ratio is 4.0 (with a 95% confidence interval of 2.8–5.8) in the mass range 0.03–1 M_⊙, which is within the expected range of ∼3–6 found in the literature. This places Tr 14 near the middle of the observed distribution of star-to-BD ratios, suggesting that in this specific mass range, BDs form at least as efficiently as low-mass stars, consistent with other environments. However, our IMF shows a significantly flatter slope and a deficit of objects at even lower masses (below 0.03 M_⊙), which are not included in this ratio. This may indicate that environmental effects primarily influence the formation efficiency of the lowest-mass BDs.

Data availability

Table 1 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/706/A58.

Acknowledgements

This work was supported by/supported in part by the Croat-ian Science Foundation under the project number HRZZ-MOBDOK-2023-9579. K.M. acknowledges support from the Fundação para a Ciência e a Tecnologia (FCT) through the CEEC-individual contract 2022.03809.CEECIND, and the Scientific Visitor Programme of the European Southern Observatory (ESO) in Chile. V.A.-A. acknowledges support from the INAF grant 1.05.12.05.03. Based on observations (programme ID GS-2019A-DD-107, PI: M. Andersen) obtained at the international Gemini Observatory, a programme of NSF NOIR-Lab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the U.S. National Science Foundation on behalf of the Gemini Observatory partnership: the U.S. National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministe-rio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

---

Appendix A Natural guide stars

In Table A.1 we list the positions of the natural guide stars.

Appendix B Photometry parameters

In Table B.1 we list the parameters used in the extraction of photometry using SOURCEEXTRACTOR and PSFEX.

Appendix C Besançon Galactic model photometric errors

The expressions for photometric errors, in corresponding bands, given to the Besançon Galactic model for the expected foreground population:

σ_J = 0.116 + exp(1.15 · J − 28.5)

σ_H = 0.13 + exp(1.3 · H − 31)

${\sigma }_{K_s}$ = 0.13 + exp(1.3 · K_S − 30).

The photometric errors, for corresponding bands, for the background population are

σ_J = 0.116 + exp(c_J · J_r − b_J)

σ_H = 0.13 + exp(c_H · H_r − b_H)

${\sigma }_{K_s}=0.13+\exp (c_{K_s}\cdot K_{s_r}-b_{K_s}).$

The coefficients are given in Table C.1 for the specific extinction value that we use and the index r notes that the values of the magnitudes used are reddened for the A_V.

All Tables

All Figures