The Initial Mass Function in the Coma Berenices Dwarf Galaxy from Deep Near-infrared HST Observations

We utilized the WFC3/IR instrument in the F110W and F160W filters. We imaged the galaxy in a 2 × 2 tile mosaic, and reduced and analyzed the four tiles independently of each other. The program consisted of 44 HST orbits, equally divided among the 4 tiles. The 11 orbits per tile were split with 4 orbits in F110W and 7 in F160W. Each orbit was split into 2 exposures, thus resulting in 8 dithered exposures in F110W and 14 dithered exposures in F160W per tile. The resulting exposure times are 10994 s in F110W and 18989 s in F160W per tile.

WFC3/IR images in filters that are sensitive to the Earth's atmospheric Hei 1.083 μm line afterglow can suffer from a time-variable background. This is due to the change in the depth of the atmosphere along the line of sight, and is especially severe at the ingress/egress of HST from occultation during its orbit around the Earth. This problem may affect F110W observations, and thus impact the up-the-ramp fitting performed by the calwfc3 pipeline, which assumes a constant background, thus producing images with artificially increased noise. To reduce the problem, the individual up-the-ramp non-destructive reads of a WFC3/IR Multiaccum exposure can be equalized by subtracting the average sky value measured from the difference between two consecutive reads. We used software developed by a former WFC3 team member (B. Hilbert 2018, private communication) to treat the images before running the standard calwfc3 reduction process (see also Brammer 2016, for a summary of the He I related problem and a very similar mitigation strategy, as well as publicly available software performing similar pre-processing of WFC3/IR images).

The images were aligned and drizzled using the standard Drizzlepac tools (Gonzaga et al. 2012), using the F110W band as the astrometric reference. Images were drizzled to rotation = 0 (North up), a scale of 60 mas/pix, and pixfrac of 0.8. The latter controls the fraction by which the input pixels are shrunk before being drizzled onto the output image grid, i.e., the drizzle "drop" size. The final drizzled images were converted to count images in preparation for use with the standalone DAOPHOT-II package (Stetson 1987).

PSF-fitting photometry was conducted using the standalone DAOPHOT-II package to obtain F110W and F160W magnitudes calibrated in the STMAG system. Thresholds in the χ and sharpness values are used to distinguish point sources from resolved background objects. The combined color–magnitude diagram (CMD) for the four tiles is shown in Figure 1.

Zoom In Zoom Out Reset image size

To test our ability to resolve background sources, we compare the luminosity function of resolved objects, i.e., all the photometric detections that do not make the χ and sharpness cut, to that of the Hubble Ultra deep Field (HUDF), using the catalog made available by Rafelski et al. (2015). Because there is no F110W observation for the HUDF, we simulate F110W data by linearly interpolating the fluxes in the F105W and F125W.

Figure 1 shows the HUFD sources, scaled to the total area of our observations. The luminosity function of our extended, background sources matches that in the HUDF very well down to about F160W = 28.6 mag, below which our extended background source luminosity function falls below the HUFD one. This means that we fail to characterize background objects as extended sources beyond this point. Even though our 50% completeness limit for point sources is about 0.8 mag fainter than 28.6 mag in F160W (see Section 2.4), we choose to preserve the purity of our catalog by limiting our analysis to stars brighter than F160W = 28.6 mag, where we can safely rely on our χ-sharpness selection. This magnitude value corresponds to about 0.23 M_⊙ at the distance and extinction of Com Ber. The mass at the 50% completeness limit, F160W ∼29.4 mag, is instead about 0.17 M_⊙, thus the loss in mass range corresponding to the adopted conservative cut in magnitude is small but not negligible. The half-mass–radius for a typical dwarf galaxy, 300 pc, is unresolved in WFC3/IR beyond ∼1 Gpc, corresponding to a distance modulus of 40 mag and to a redshift, z ∼ 0.15. Dwarf galaxies of intrinsic brightness M₁₆₀ ∼ −10 mag at such a distance contribute to the unresolved background sources seen in Figure 1, with brighter dwarf galaxies contributing up to even larger distances. This problem will be alleviated (by a factor of three at fixed wavelength) when the James Webb Space Telescope becomes operational.

We further select our data by excluding the sources on the red giant branch (RGB). In practice, we adopt a bright magnitude cut at F160W = 23 mag. In our simulations versus the data comparison, this cut allows us to ignore the RGB region, where the stellar models tend to not be able to reproduce the observed colors. Given the faster evolution on the RGB relative to the main sequence, this bright limit has no significant impact on the mass range probed. The final numbers of selected stars per tile are 223, 163, 149, and 182, respectively, for a total of 717.

The selected data are highlighted in a CMD in Figure 2, where the IR data presented in this paper are also compared to the optical Com Ber data first presented by Brown et al. (2014), and used for deriving the IMF in Gennaro et al. (2018). The figure shows how the new data reach a lower stellar mass, thus giving more mass leverage for the IMF determination with respect to the existing optical data. The optical observations cover, however, a larger area on the sky, and thus the number of Com Ber members is larger in the optical than in the IR CMD.

Zoom In Zoom Out Reset image size

We perform artificial star tests to obtain an estimate of the noise distribution of our sources. We use these experiments in our synthetic CMD simulation when fitting for the IMF parameters. The artificial star tests were produced by randomly populating isochrones at [Fe/H] = −1.6, −2.0, and −2.4 dex and age = 10, 11, 12, and 13 Gyr. Each pass of artificial star tests inserted 1000 stars into the images and blindly recovered them, comparing output to input. Specifically, each pass drew 800 stars from the isochrones, with significant scattering (to fill in the CMD parameter space), plus another 200 stars at random locations over the entire CMD (to avoid having gaps in the artificial star CMD coverage). The result places the most weight where we need it (i.e., near real isochrones, where most stars should fall) but characterizes the uncertainties and incompleteness over a broad swath of CMD space. The same χ and sharpness criteria applied to the real stars are applied to the artificial stars, and if the latter do not pass the χ-sharpness test, they are considered undetected.

Figure 3 shows the luminosity functions of the inserted stars (input magnitudes in each band) and the luminosity function of the detected ones (again input magnitudes), as well as the ratio of the two. The vertical dashed line in the m₁₆₀ bottom panel corresponds to 28.6 mag, i.e., the limit at which we lose the ability to detect extended background sources in our data. This is the faint cut we have to impose to our catalog in order to avoid contamination from unresolved background galaxies. It is clear that our ability to detect point sources extends well past this limit, with the 50% completeness level being about 29.4 mag in F160W. However, we choose to keep the highest possible purity in our catalog, at the cost of losing some depth.

Zoom In Zoom Out Reset image size

During the analysis stage we observed that the photometric uncertainty encoded in the artificial stars was slightly underestimated. In particular, the scatter of simulated CMDs, built using the artificial star tests and realistic metallicity and age distributions for Com Ber, showed a narrower main sequence than what is observed. This extra error may be due to the limitations of the PSF model constructed from the undersampled data. Another possibility is a metallicity effect. While we are using the most up-to-date metallicity distribution function (MDF) available for Com Ber (Brown et al. 2014), if the true intrinsic distribution is slightly broader, it may produce extra CMD width, with respect to what we can reproduce using our adopted MDF.

We determined that an extra error term of 0.0175 mag (standard deviation) in each photometric band was enough to force the simulations to agree with the data, in terms of simulated main sequence width. This extra error is implemented in our subsequent analysis by further scattering each artificial star's output magnitude using a draw from a normal distribution with a width equal to 0.0175 mag. It is worth noting that this error term is quite small, on the order of the photometric repeatability quoted by the WFC3/IR instrument team (∼0.5%–1%, see the WFC3 data handbook⁸ ), thus our heuristic approach of adding an extra random Gaussian term is appropriate and sufficient for our purposes.

Brown et al. (2014) showed a very good agreement of the stellar models with the M92 cluster ridge line in the HST/ACS optical (F606W–F814W, F814W) CMD. M92 is the most metal-poor cluster for which we have good WFC3/IR data from the HST Galactic Bulge Treasury Program (GO: 11664, PI: T. M. Brown). Such data can be used for calibration of the stellar models in the F110W and F160W bands. We do not observe the same agreement of the models with the M92 IR data. However, rather than fine-tuning the isochrone parameters to try and obtain a better fit of the M92 ridge line, we assume that the cluster parameters for which the models give good agreement in the optical are valid: (m − M)₀ = 14.62 mag, E(B − V) = 0.023 mag, age = 13.2 Gyr, [Fe/H] = −2.3 dex. We derive an empirical calibration of the IR models with respect to the M92 ridge line in the following steps:

• 1.  
  We generate a ridge line for the M92 (F110W, F160W) data, using the technique described in Correnti et al. (2016); this step is shown in the top right panel of Figure 4.
• 2.  
  We subtract the extinction in each band and distance modulus from the M92 ridge line using (m − M)₀ = 14.62 mag, E(B − V) = 0.023 mag, R_V = 3.1, and assuming that the extinction law can be described by a Cardelli et al. (1989) model. This results in (A₁₁₀, A₁₆₀) = (0.0228, 0.0144) mag at the respective effective wavelengths of each of the F110W and F160W filters. The distance and extinction-subtracted ridge line is shown as "Absolute M92 RL" in the top right panel of Figure 4.
• 3.  
  We compute the isochrone for the nominal M92 values of age (13.2 Gyr), [Fe/H] (−2.3 dex) and [α/Fe] (0.4 dex), with [O/Fe] = +0.64 dex. This is shown as "original" in the top right panel of Figure 4.
• 4.  
  We identify "characteristic points" on both the M92 ridge line and the original isochrone. These are visible features in the two curves, and we assume that, for a calibrated model, their color and magnitudes need to match exactly. These characteristic points are marked by dots in the top right panel of Figure 4 (the dots for the corrected model lie, by definition, on the ridge line). Some of these points are identifiable with significant evolutionary stages (i.e., the main sequence turnoff, the low-mass main sequence kink, the base of the RGB). Others have been chosen based on morphological characteristics of the isochrones (e.g the mid-point in the sub-giant branch, slight bends along the main sequence between the kink and turnoff).
• 5.  
  We compute the differences in F160W magnitude and F110W–F160W color between the points corresponding to the same stages along two curves.
• 6.  
  We interpolate such corrections as a function of $\mathrm{log}g$ in the model.
• 7.  
  We assume that the same correction in $\mathrm{log}g$ applies for all ages and metallicities.

Zoom In Zoom Out Reset image size

The above sequence of steps does not yet provide a complete agreement when the models for the peak of the Com Ber MDF are computed and overplotted on the data.

We thus define a second calibration step, this time in apparent magnitude space. The steps are as follows:

• 1.  
  We build a ridge line for the Com Ber population. This is shown in the bottom left panel of Figure 4.
• 2.  
  We compute the isochrone for the bulk of the Com Ber population, corresponding to the peak of the metallicity distribution of the galaxy, [Fe/H] = −2.6 dex and at an age of 14.0 Gyr, and at its distance, (m − M)₀ = 17.956 mag, and extinction, A_V = 0.124 (all values from Brown et al. 2014). From A_V we derive (A₁₁₀, A₁₆₀) = (0.0398, 0.0251)mag. This is the isochrone labeled as "original" in the lower right panel of Figure 4.
• 3.  
  We apply the empirical calibration described in Section 3.1 to the original isochrone, thus obtaining the isochrone labeled as "1st corr." in the lower right panel of Figure 4.
• 4.  
  We compute residual color offsets between the "1st corr." isochrone and the Com Ber Ridge Line.
• 5.  
  We interpolate the residual offset as a function of m₁₆₀.

The combination of the two calibration steps gives us models that agree very well with the Com Ber ridge line. The residual offsets visible in the lower right panel of Figure 4, between the final model (orange) and the ridge line(gray), are of the order of ≲0.01 mag.

It must be noted that the extent of the discrepancy between the uncorrected models and the M92 or Com Ber ridge lines is of the order of 0.05 mag maximum, and only in the coolest regions of the CMD. This is probably related to uncertainties in the model atmospheres of cool stars in the near-IR. Around 0.3 M_⊙, where the most serious discrepancies arise, the mass–luminosity relation in our passbands is of about 0.1 M_⊙/mag, thus an error of 0.05 mag in color would correspond to only .005 M_⊙. The impact of such errors on the IMF determination is negligible.

From Figure 6 we can see that the three IMF models give very similar results. The m₁₆₀ luminosity functions for the best-fit parameters for each of the three models (left columns in the figure) do trace the observed luminosity function, even with some differences in the details. From the right panels of Figure 6, the best-fit parameters for the three IMF models produce very similar PDMFs (green), regardless of the different underlying IMFs (orange). We note some discrepancy between the models, especially at the high-mass end, where there is the fewest stars. Another fact emerges quite clearly: no obvious turnover of the IMF is observed. The SPL model gives as good a fit to the data down to our limiting mass as the BPL and LN ones; if anything, the best-fit SPL model seems to be the one producing the luminosity function giving the best qualitative agreement at faint magnitudes.

To provide a more quantitative assessment of the goodness of fit (gof) and see whether any of the three models is inconsistent, in a gof sense, with the data, we devised an ad hoc strategy. Our fitting method uses an ABC approach, and thus bypasses an explicit calculation of the likelihood. Methods like Bayes Factors, Akaike Information or Bayesian Information criteria (AIC, BIC) cannot therefore be used here. Our acceptance of a proposed parameter set, ${{\boldsymbol{\theta }}}_{i}$, relies on generating a CMD using ${{\boldsymbol{\theta }}}_{i}$, which we call y_i and determine whether y_i is close to the observed CMD, y_obs. The distance between the simulated and observed CMD is given by Equation (10) of Gennaro et al. (2018); we will call this distance ρ(y_i, y_obs).

Our gof criterion is a modification of the gof statistic based on the posterior predictive distribution proposed by Lemaire et al. (2016), and it consists of the following:

• 1.  
  Take two subsets of draws from the posterior, a training set $\{{{\boldsymbol{\theta }}}_{j}\}{}_{j\in J}$, and a test set $\{{{\boldsymbol{\theta }}}_{k}\}{}_{k\in K}$, possibly from two disjointed parts of the MCMC chain. In practice, we use the first 250 draws for J and the last 1000 for K.
• 2.  
  Compute, for each k ∈ K, the average distance of the corresponding CMD, ${{\boldsymbol{\theta }}}_{k}$ to the $\{{{\boldsymbol{\theta }}}_{j}\}{}_{j\in J}$:

  $$\begin{eqnarray}&&{D}_{\mathrm{post},k}=\displaystyle \frac{1}{{n}_{J}}\displaystyle \sum _{j\in J}\rho ({y}_{k},{y}_{j}),\end{eqnarray} \tag{ 4 }$$

  where n_J is the size of the training set.
• 3.  
  Similarly, compute the average distance to the training set for the observed CMD:

  $$\begin{eqnarray}&&{D}_{\mathrm{post},\mathrm{obs}}=\displaystyle \frac{1}{{n}_{J}}\displaystyle \sum _{j\in J}\rho ({y}_{\mathrm{obs}},{y}_{j}).\end{eqnarray} \tag{ 5 }$$
• 4.  
  See where D_post,obs falls within the distribution of D_post computed from the test set.
• 5.  
  If the average distance for the observations falls outside some pre-defined quantile (e.g., if the observations are on average more distant than 95% of the $\{{{\boldsymbol{\theta }}}_{k}\}{}_{k\in K}$), rule out that model.

Figure 7 shows the results for each of the three IMF models. In each case the D_post,obs falls well within the bulk of the distribution of D_post. For none of the three models is the fit outstandingly better or worse than that for the others. All three models are not inconsistent with the data according to our gof test. This confirms the qualitative conclusion that could be drawn from the luminosity functions of Figure 6: all three IMF models can produce CMDs that are consistent with the observed CMD.

Zoom In Zoom Out Reset image size

When compared to the results of Gennaro et al. (2018), our current results for the SPL and LN model IMF are in good agreement, namely the SPL slope and the LN m_c and σ are within the 68% C.I. (see Figure 8). This is an important sanity check for both works, demonstrating that any bias in the results due to the different mass ranges probed, as well as to uncertainties in the adopted color–temperature and/or the mass–luminosity relations, if present, must be negligible with respect to the uncertainties in the parameter estimate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We cannot perform a comparison with the broken power-law model, since Gennaro et al. (2018) did not attempt a BPL modeling of the IMF to fit their optical data. This was because their data only reached down to ∼0.4–0.5 M_⊙, on average, for the six galaxies of their sample.

The binary fractions are generally in worse agreement between the current work and Gennaro et al. (2018). The derived binary fraction is quite sensitive to the observed width in the CMD. As illustrated in Section 2.5, we were forced to add an extra error in order for our simulated CMDs to reproduce the observed CMD width. We used the same extra error in each band and at all magnitudes, because a more detailed model for it would have implied a better "a priori" knowledge of what the width of the simulated CMD should have looked like, i.e., knowledge of the binary fraction, among other things. We speculate that the necessarily qualitative-only assessment of the extra errors may have caused an extra spread of the simulated main sequence, which in turn translates into a smaller binary fraction, to compensate for the extra width. On the other end, in Gennaro et al. (2018), Com Ber was the only UFD with a best-fit binary fraction higher than 50%, while for the other UFDs the best-fit binary fraction was always smaller than 30%. Thus, it could be the case that the binary fraction was overestimated in that work as well. We do not try to draw further conclusions from discrepancies between the derived binary fractions.

If they were limited to the SPL parameterization, our results would suggest that the mean IMF slope of Com Ber is somewhat shallower than the Milky Way classical Salpeter slope of −2.35. At the same time, given the limited mass range and the small number of stars (717), the significance of such a result is limited. The Salpeter value falls outside the 95% credible interval, but inside the 99% one.

The BPL and LN parameterizations give results that are more similar to the Milky Way for the same type of functional forms, i.e., Kroupa (2001) for the BPL and Chabrier (2003) or Bochanski et al. (2010) for the LN form. While the best-fit parameters for both the BPL and LN model would result in a slightly bottom-light IMF for Com Ber, with respect to the Galaxy, the Kroupa (2001) and Chabrier (2003) values are within the 68% confidence interval.

Gennaro et al. (2018), based on slope values from SPL model fits, claimed that, on average, a group of six UFDs showed a more bottom-light IMF with respect to the Milky Way, i.e., fewer dwarf stars. It is possible that the differences between the IMF for the UFDs with respect to the Galaxy claimed by Gennaro et al. (2018) are at least partly an artifact of using a SPL model. Gennaro et al. (2018) studied the mass range above 0.4–0.5 M_⊙, at the edge where Kroupa (2001) placed the power-law break, and claimed to not expect a significantly biased result for an SPL model. On the other hand, El-Badry et al. (2017) showed how adopting an SPL model could artificially lead to shallow slopes if the true underlying IMF has an LN shape and data do not probe down to the characteristic mass. El-Badry et al. (2017) obtained a slope of ∼−1.55 by fitting a Chabrier (2003) IMF with an SPL model, but limiting the fit to the (0.4–0.77) M_⊙ range, similar to what was obtained by Gennaro et al. (2018) for Com Ber. The fact that the results of the current work, for the SPL model, using an extended mass range that reaches very close to the Chabrier (2003) characteristic mass, are very similar to Gennaro et al. (2018), shows, however, that any bias, if present, must be smaller than the uncertainty introduced by the limited number of stars used.

It is worth noting that in Gennaro et al. (2018) Com Ber was, together with Boo I, the UFD for which the IMF best-fit parameters for both an SPL and LN model were closer to those of the Milky Way. Other galaxies in Gennaro et al. (2018), in particular Hercules, Leo IV, and CVn II, showed a much more significant discrepancy, and mean slope values close to −1. Also, the characteristic mass for the LN model for these three galaxies was much higher than the Chabrier (2003) value, thus for those galaxies the discrepancy with the Milky Way was present even when using the LN parameterization (although, again, with lower significance than when using an SPL model).

Gennaro et al. (2018) noted that Com Ber and Boo I, being the closest UFDs in their sample, subtended a larger angle in the sky, and thus their CMDs could have been more contaminated by unresolved background galaxies, especially at the faint end, making their IMFs artificially steeper. In the current work we have tried to minimize the impact of background interlopers by limiting ourselves to m₁₆₀ < 28.6 mag, where we can still discriminate between extended background galaxies and point sources. This does not ensure 100% purity in our sample, and thus may be part of the reason for having an IMF that does not differ from the Milky Way as significantly as for other UFDs.

In Section 6.2 we illustrated some possible concerns about the reliability of the estimated binary fraction, in relation to the extra error added to our artificial star tests in order to reproduce the observed width of the CMD. With those caveats in mind, it is still interesting to dissect the multi-dimensional posterior probability distribution functions in order to emphasize the correlation of the other parameters with the binary fraction itself.

For each of the SPL, LN, and BPL models, we have divided the posterior draws into 4 parts, corresponding to the 0% to 25%, 25% to 50%, 50% to 75%, and 75% to 100% quantiles of the marginal distribution of the binary fraction. The edges of the quantiles change slightly between the different models, according to the detailed difference in the marginal distribution of the binary fraction.

As visible in Figure 9, the two parameters that correlate the most with the binary fraction are the slope in the SPL model and the characteristic mass, m_c, in the LN model. The correlation goes in the same direction for both models: a smaller binary fraction corresponds to a more bottom-light IMF. At increasingly high binary fractions, the marginal distributions for the slope in the SPL model and m_c in the LN model become more consistent with what is observed in the Milky Way. A similar but less prominent correlation is observed for the slopes and the break mass in the BPL case. Even in this case, when considering a higher binary fraction, the IMF becomes more "bottom-heavy" or more Galaxy-like.

While this is an interesting exercise, we do not have any prior knowledge that allows us to independently constrain the binary fraction, and thus narrow down the uncertainty on the parameters that correlate with it. Our fitting procedure allows us to marginalize over the uncertain binary fraction and report more realistic uncertainties on the other parameters, than if we were to arbitrarily fix them to, e.g., the typical values observed in the Galaxy of about 30% (Raghavan et al. 2010).

A further caveat to consider is the way in which binary pairing is implemented in our simulations. As explained in Gennaro et al. (2018), our IMF represents the distribution of system masses. In simulating synthetic CMDs, we draw system masses, M_sys, from an IMF model, then decide whether a system should be a single star or a binary, according to a probability set by the binary fraction. If the system is labeled as a binary, a mass ratio, q = M₂/M₁, is drawn from a uniform distribution between 0 and 1. Individual masses are given by ${M}_{1}=\tfrac{{M}_{\mathrm{sys}}}{1+q}$ and ${M}_{2}=\tfrac{{{qM}}_{\mathrm{sys}}}{1+q}$ . This is different from drawing individual masses from the same IMF and pairing randomly. The qualitative argument for our choice is that it is unlikely that the two stars in a binary are completely independent draws from the same IMF. For progressively tighter binaries it is reasonable to assume that at formation the two stars affected each other's amount of accreted mass from the protostellar phase.

In determining the mass ratio distribution, we used the fact that for "solar-type" stars, the binary mass ratio is uniformly distributed between 0 and 1 (Duquennoy & Mayor 1991; Raghavan et al. 2010). Such a distribution may not be valid outside of the Milky Way or in different mass regimes. We thus stress that, while we use a reasonable assumption, the derived fraction of binaries depends on the adopted mass-ratio distribution.

Nevertheless, we performed a sanity check on the adopted rule of binary component pairing. We verified that our process for generating binaries, which we refer to as correlated draws, does not alter the inferred slope, with respect to a method in which the two stars in a binary are generated from the same single-star IMF, referred to as independent draws. We assume the same input slope for the single-star IMF and the system IMF in the independent and correlated cases, respectively. We then compare the single-star (primaries plus secondaries) system, as well as the primary star and secondary star distributions. Differences appear at the high and low ends of the mass spectrum. The differences at the low-mass edge are due to the fact that in the correlated case, secondaries can in principle have M₂ → 0, while a finite limit must be set for extracting individual masses from a power law in the independent case. At the high-mass edge there are differences due to the definitions of upper limits, i.e., we cannot enforce the same limit on both the maximum individual mass and the maximum system mass consistently for both the independent and correlated scenarios. When limiting the analysis to the mass region accessible to us, i.e., far from both the low-mass and high-mass end of the IMF, no significant differences appear for either the mass ratio distribution, the individual star distribution, or the system mass distribution between the two methods of generating binaries. As an example, a two-sample Kolomogorov–Smirnov test for 1000 draws from each of the independent and correlated cases, giving a p-value of 0.95; thus the null-hypothesis that the distributions of the two samples are the same cannot be rejected.