Dynamical masses of brightest cluster galaxies – II. Constraints on the stellar IMF Free

ABSTRACT

We use stellar and dynamical mass profiles, combined with a stellar population analysis, of 32 brightest cluster galaxies (BCGs) at redshifts of 0.05 ≤|$z$| ≤ 0.30, to place constraints on their stellar initial mass function (IMF). We measure the spatially resolved stellar population properties of the BCGs, and use it to derive their stellar mass-to-light ratios (⁠|$\Upsilon _{\star \rm POP}$|⁠). We find young stellar populations (<200 Myr) in the centres of 22 per cent of the sample, and constant |$\Upsilon _{\star \rm POP}$| within 15 kpc for 60 per cent of the sample. We further use the stellar mass-to-light ratio from the dynamical mass profiles of the BCGs (⁠|$\Upsilon _{\star \rm DYN}$|⁠), modelled using a multi-Gaussian expansion and Jeans Anisotropic Method, with the dark matter contribution explicitly constrained from weak gravitational lensing measurements. We directly compare the stellar mass-to-light ratios derived from the two independent methods, |$\Upsilon _{\star \rm POP}$| (assuming some IMF) to |$\Upsilon _{\star \rm DYN}$| for the subsample of BCGs with no young stellar populations and constant |$\Upsilon _{\star \rm POP}$|⁠. We find that for the majority of these BCGs, a Salpeter (or even more bottom-heavy) IMF is needed to reconcile the stellar population and dynamical modelling results although for a small number of BCGs, a Kroupa (or even lighter) IMF is preferred. For those BCGs better fit with a Salpeter IMF, we find that the mass-excess factor against velocity dispersion falls on an extrapolation (towards higher masses) of known literature correlations. We conclude that there is substantial scatter in the IMF amongst the highest mass galaxies.

1 INTRODUCTION

One of the biggest unknowns of stellar population and mass studies of galaxies is the stellar initial mass function (IMF), the distribution of stellar masses at birth. Limited information regarding the shape and properties of the IMF can be constructed from first principles, e.g. using analytical models (Hopkins 2013) or hydrodynamical simulations (Krumholz et al. 2016), and must mainly be empirically determined. First thought to be universal (Bastian, Covey & Meyer 2010), the last decade has seen increasing evidence suggesting that the IMF varies between galaxies, as well as within individual galaxies (e.g. Conroy & van Dokkum 2012; Cappellari et al. 2013; Ferreras et al. 2013; Martín-Navarro et al. 2015a; Lyubenova et al. 2016; Davis & McDermid 2017).

Since some stellar absorption line features are stronger in dwarf than giant stars (Spinrad & Taylor 1971), one method to constrain the IMF (from integrated stellar light) is stellar population synthesis and direct measurements of IMF-sensitive spectral absorption features (Conroy & van Dokkum 2012; Spiniello et al. 2012; La Barbera, Ferreras & Vazdekis 2015; Rosani et al. 2018; La Barbera et al. 2019). Another method is to dynamically measure the stellar mass-to-light ratio (Υ_⋆), by assuming the system is dynamically relaxed and by making assumptions about the contribution of the dark matter and supermassive black hole to the gravitational potential of the galaxy (Thomas et al. 2011; Cappellari et al. 2013; Tortora, Romanowsky & Napolitano 2013). A third independent method to constrain the IMF is through gravitational lensing (Auger et al. 2010; Treu et al. 2010; Dutton et al. 2013; Smith, Lucey & Conroy 2015; Collier, Smith & Lucey 2018).

Studies utilizing a combination of these independent techniques have found evidence for an IMF that is ‘heavier’, i.e. with larger stellar mass-to-light ratios than those predicted by Chabrier (Chabrier 2003) or Kroupa (Kroupa 2001) IMFs, in the most massive early-type galaxies (Treu et al. 2010; Cappellari et al. 2012; Conroy & van Dokkum 2012; La Barbera et al. 2013; Tortora et al. 2013). The excess mass appears to increase with the stellar velocity dispersion (a proxy for the mass, or similarly with correlated properties e.g. metallicity or [α/Fe]-enhancements) of the galaxy (Conroy & van Dokkum 2012; La Barbera et al. 2015; Rosani et al. 2018). This qualitative agreement between different techniques has lent confidence to claims of a non-universal, heavy IMF for massive early-type galaxies.

However, Smith et al. (2015) investigated early-type galaxy strong lenses from the SINFONI Nearby Elliptical Lens Locator Survey (SNELLS), and based on their lensing masses, concluded that these galaxies have a stellar IMF consistent with that of the Milky Way, and not the bottom-heavy IMF typically reported for massive early-type galaxies. Leier et al. (2016) also exclude a very bottom-heavy IMF in their study of 18 massive early-type galaxies from the Sloan Lenses ACS Survey (SLACS). Smith, Lucey & Edge (2017) further present dynamical modelling of the BCG in Abell 1201. By using a combination of lensing and stellar dynamics, and by imposing a standard NFW dark matter density profile (Navarro, Frenk & White 1996), they recover a stellar mass-to-light ratio that is consistent with a Milky Way-like IMF.

Smith (2014) pointed out that the IMF constraints obtained for individual galaxies using different techniques do not always agree, and that different studies do not agree on the underlying principal galaxy property tied to the IMF trends (see also Newman et al. 2017). However, the study by Lyubenova et al. (2016), which constrained the IMF of 27 early-type galaxies using stellar populations and dynamics from CALIFA data, suggested the opposite, i.e. that there are no disagreements on a case-by-case basis.

In this paper, we focus on the most massive early-type galaxies. We use a large sample of 32 nearby (⁠|$z$| < 0.3) brightest cluster galaxies (BCGs), from the Multi-Epoch Nearby Cluster Survey (MENeaCS) and Canadian Cluster Comparison Project (CCCP) cluster samples. The BCG sample spans M_K = –25.7 to –27.8 mag, with host cluster halo masses M₅₀₀ between 2.0 × 10¹⁴ and 1.5 × 10¹⁵ M_⊙ (Herbonnet et al. 2020). In an accompanying paper (Loubser et al. 2020), we have modelled the stellar and dynamical mass profiles of 25 out the 32 BCGs (excluding the seven for which we found multiple nuclei or significant substructure in their nuclei) using an adapted multi-Gaussian expansion (MGE; Monnet, Bacon & Emsellem 1992; Emsellem, Monnet & Bacon 1994; Cappellari 2002) technique and Jeans Anisotropic Method (JAM; Cappellari et al. 2006; Cappellari 2008) for an axisymmetric case, deriving the stellar mass-to-light ratio (⁠|$\Upsilon _{\star \rm DYN}$|⁠), and stellar velocity anisotropy (β_|$z$|).

One major uncertainty in stellar dynamical modelling is the contribution of the dark matter halo to the potential of the galaxies. We were able to estimate the dark matter mass (M_DM) from weak lensing observations (Herbonnet et al. 2020) for 23 of the 25 MENeaCS and CCCP clusters, and thereby limit the number of free parameters in the dynamical models. This, together with our large sample size, is a considerable improvement on previous studies modelling the stellar dynamics of BCGs.

Here, we present the spatially resolved stellar population properties of (all 32 of) the BCGs, and use it to estimate their stellar mass-to-light ratios (⁠|$\Upsilon _{\star \rm POP}$|⁠), under the assumption of both a Salpeter and a Kroupa IMF. We then compare the stellar mass-to-light ratios derived from the two independent methods (⁠|$\Upsilon _{\star \rm POP}$| with |$\Upsilon _{\star \rm DYN}$|⁠) for a subsample of BCGs (i.e. those with no young stellar populations and a non-variable |$\Upsilon _{\star \rm POP}$| over our radial range) and use it to constrain their stellar IMF.

We use H₀ = 73 km s⁻¹ Mpc⁻¹, Ω_matter = 0.27, and Ω_vacuum = 0.73 throughout. We refer to velocity dispersion as σ, to rotational velocity as V, and to |$\sqrt{V^{2} + \sigma ^{2}}$| as the second moment of velocity (ν_rms). We also derive the stellar mass-to-light ratios, |$\Upsilon _{\star \rm DYN}$| and |$\Upsilon _{\star \rm POP}$|⁠, in the rest-frame r band. In Section 2, we briefly summarize our data as well as the dynamical modelling procedure used to derive the stellar mass-to-light ratios (⁠|$\Upsilon _{\star \rm DYN}$|⁠) in Loubser et al. (2020). In Section 3, we present the spatially resolved stellar population properties, and the predicted stellar mass-to-light ratios (⁠|$\Upsilon _{\star \rm POP}$|⁠), assuming a Salpeter IMF (Salpeter 1955). We then also fit the stellar populations using a Kroupa IMF (Kroupa 2001), and discuss the constraints on the IMF in Section 4. We correlate the mass-excess factor (⁠|$\alpha = \log _{10} (\Upsilon _{\star \rm DYN}) - \log _{10} (\Upsilon _{\star \rm POP})$|⁠) with other galaxy properties in Section 5, and summarize our conclusions in Section 6.

2 DATA AND STELLAR MASS-TO-LIGHT RATIOS FROM DYNAMICS (⁠|$\Upsilon _{\star \rm DYN}$|⁠)

We use spatially resolved, long-slit spectroscopy for 14 MENeaCS and 18 CCCP BCGs, observed on the Gemini North and South telescopes, and r-band imaging observed on the Canada–France–Hawaii telescope (CFHT). We also use host cluster properties derived from Chandra/XMM–Newton X-ray data, and cluster masses measured through weak lensing (Mahdavi et al. 2013; Hoekstra et al. 2015; Herbonnet et al. 2020).

The stellar population analysis and star formation histories of the CCCP BCGs are presented in Loubser et al. (2016), but since we present the stellar population analysis for the MENeaCS BCGs here for the first time, we briefly summarize the relevant properties of the long-slit spectroscopic data used for the MENeaCS analysis. The 14 MENeaCS BCGs were observed using the Gemini Multi-Object Spectrograph (GMOS) detector during two semesters (GS2009A, GN2009A, GN2009B). The instrumental configuration consisted of the B600 grating at a central wavelength of 4600 Å, and a slit width of 1 arcsec. We used 2 × 2 binning, corresponding to an instrumental resolution of 71 km s⁻¹. The spatial apertures (i.e. 0–5 and 5–15 kpc) were chosen to ensure sufficient S/N for reasonable errors in the stellar population analysis for all BCGs, with the S/N of the MENeaCS spectra generally higher than that of the CCCP spectra presented in Loubser et al. (2016).

In Loubser et al. (2020), we present the Gauss–Hermite higher order velocity moments h₃ and h₄ for the BCGs, and find that the central measurements of h₄ are positive for all our BCGs. We then model the stellar and dynamical mass profiles of 25 out of the 32 BCGs (excluding the seven for which we found multiple nuclei or significant substructure in their nuclei), using an adapted MGE and JAM for an axisymmetric case (for both cylindrically and spherically aligned models), deriving the stellar mass-to-light ratio (⁠|$\Upsilon _{\star \rm DYN}$|⁠), and anisotropy (β_|$z$|),¹ where the dark matter mass was constrained from weak lensing results. Our fits to the observed kinematics are restricted to the galaxy centre, where the stellar component is the dominant contributor to the mass (<20 kpc).

Our dynamical modelling revealed that the stellar anisotropy and velocity dispersion profile slope (η, from Loubser et al. 2018) are strongly correlated. The BCGs with rising velocity dispersion profiles showed tangential anisotropy parameters, whereas the BCGs with decreasing velocity dispersion profiles showed radial anisotropy parameters. The positive velocity dispersion gradients (i.e. profiles rising from the centre of the BCG), can also arise from a significant contribution from the intracluster light (ICL). For a small subset of BCGs with positive velocity dispersion gradients, a variable |$\Upsilon _{\star \rm DYN}$| can also explain the profile shape, instead of tangential anisotropy or an ICL contribution.

3 STELLAR POPULATION ANALYSIS

In Loubser et al. (2016), we identified plausible star formation histories for the 18 CCCP BCGs for which we have long-slit spectra by fitting simple stellar populations (SSPs) and composite populations consisting of a young stellar component and an intermediate or old stellar component. In this section, we repeat the analysis for the 14 MENeaCS BCGs, located at lower redshifts. We use the University of Lyon Spectroscopic analysis Software (ulyss; see Koleva et al. 2009; Groenewald & Loubser 2014; Loubser et al. 2016) to fit stellar population models to the BCG spectra, taking into account the internal kinematics of the galaxies. We use the Vazdekis et al. (2010) stellar population models, based on the MILES library (Sánchez-Blázquez et al. 2006), first for a Salpeter IMF (Salpeter 1955), and then for a Kroupa IMF (Kroupa 2001). Any emission lines are masked, and the rest of the spectrum is then used to determine the best-fitting stellar population model. We identify both the best-fitting spectrum resulting from a single SSP, as well as the best-fitting composite stellar model comprising of a young stellar population superposed on an intermediate or old stellar population, and determine which best describe the observed spectrum.² We further describe the BCGs with young stellar components, and the stellar mass-to-light ratios (⁠|$\Upsilon _{\star \rm POP}$|⁠) below.

To assess possible systematic errors in our stellar population analysis, and the possible effect on |$\Upsilon _{\star \rm POP}$|⁠, we also use the full spectral fitting code firefly (Fitting IteRativEly For Likelihood analYsis) as described in Wilkinson et al. (2017), with the MaStar stellar population models (Maraston et al. 2020), in Section 4.2, and show that the modelling approach does not affect our main conclusions.

3.1 BCGs containing young stellar population components

We detected prominent young (∼200 Myr) stellar populations in 4 of the 18 CCCP BCGs (see Loubser et al. 2016). We repeat the analysis on the MENeaCS BCGs, with the same spatial bins as used in Loubser et al. (2016), i.e. 0–5 and 5–15 kpc, and find 3 of the 14 MENeaCS BCGs (Abell 780, 1795, and 2055) show very young stellar populations (<200 Myr) in their inner apertures. Therefore, in total 7/32 (22 per cent) of the BCG sample show a young stellar population component. Two of these BCGs, Abell 780 and Abell 2055, have intermediate/old SSP-equivalent ages in their outer bins (and thus relatively large SSP-equivalent age gradients), whereas Abell 1795 has a fairly young (∼2 Gyr) SSP-equivalent age in its outer aperture. Only two MENeaCS BCG, Abell 1991 and 2319, have very old (>10 Gyr) stellar populations in both the inner and outer apertures. The nine remaining MENeaCS BCGs all present intermediate stellar populations in their inner and outer bins. See Loubser et al. (2016) for the complete discussion on intermediate SSP-equivalent ages in BCGs. The stellar population results for the MENeaCS subsample are presented in Table 1 (here), and the CCCP subsample in Loubser et al. (2016, their table 2).

3.2 Stellar mass-to-light ratios from stellar populations (⁠|$\Upsilon _{\star \rm POP}$|⁠)

In addition to identifying and constraining young stellar populations, we use the stellar population results to determine the stellar mass-to-light ratios (⁠|$\Upsilon _{\star \rm POP}$|⁠), which depend sensitively on the stellar IMF. The BCG stellar mass-to-light ratios from dynamics (⁠|$\Upsilon _{\star \rm DYN}$|⁠), which are treated as a (constant) free parameter in the dynamical mass models, can be directly compared to |$\Upsilon _{\star \rm POP}$| (for different IMFs) to place constraints on the IMF.

We use the ages and metallicities derived (Table 1), and photometric predictions from the Vazdekis/MILES (Sánchez-Blázquez et al. 2006; Vazdekis et al. 2015) models with the Girardi et al. (2000) isochrones with solar [α/Fe]-enhancement, to determine |$\Upsilon _{\star \rm POP}$| for the r filter in the inner (0–5 kpc) and outer (5–15 kpc) apertures. For the BCGs with younger components, we derive a light-weighted average |$\Upsilon _{\star \rm POP}$| of the composite stellar populations. If the |$\Upsilon _{\star \rm POP}$| derived in the inner and outer bins agree within their 1σ errors (as propagated from the errors on the ages and metallicities), then we consider the |$\Upsilon _{\star \rm POP}$| to remain constant with radius within this spatial range in the BCG.

We present these results (for a Salpeter IMF) in Table 2 for all 32 BCGs. We find that 19/32 (60 per cent) of the BCGs have constant |$\Upsilon _{\star \rm POP}$| over this radial range (0–15 kpc). Since the difference between the stellar population properties derived for the inner and outer apertures are relative, whether |$\Upsilon _{\star \rm POP}$| is relatively insensitive to the systematic uncertainties connected to the stellar population models or the assumed IMF.

3.3 Discussion of stellar population modelling results

We plot the stellar mass-to-light ratio (⁠|$\Upsilon _{\star \rm POP}$|⁠) against central cooling time (t_c,0) in Fig. 1. Blue symbols indicate BCGs with young components, green intermediate aged, and red old stellar populations. The solid symbols are those galaxies for which the |$\Upsilon _{\star \rm POP}$| is constant between the inner and outer apertures (within the errors), and the empty symbols non-constant (primarily driven by age gradients between the inner and outer apertures). From Fig. 1 it follows, similar to our conclusions in Loubser et al. (2016), that the BCGs with young stellar populations are located in host clusters with short cooling times, with the exception of the BCG in Abell 2055. The BCG in Abell 2055 (top left corner in Fig. 1) shows optical emission lines in its spectrum, a blue core (Bildfell et al. 2008), and a young stellar component, however, the host cluster has the longest cooling time of all the clusters considered here. This BCG hosts a BL Lac point source that is the dominant contributor to the observed emission (Green et al. 2017).

Unlike the BCGs in Abell 780, 1795, and 2055, which show young components as well as emission lines in their spectra, the BCG in Abell 646 shows some emission lines in its spectrum, but a stellar component in the inner aperture of age ∼3.7 Gyr. It also has a blue core as determined by Bildfell et al. (2008). There is no significant difference between the χ²-values for the stellar population fits in the inner aperture for an SSP and a two-component composite stellar population fit, so we default to an SSP fit (as described in Loubser et al. 2016). This is thus an intermediate component galaxy, but with a stellar population component (in the inner bin) on the limit of where a small fraction of young stars cannot be confidently detected (as described in detail in Loubser et al. 2016).

We do not find any strong correlations between stellar mass-to-light ratio (⁠|$\Upsilon _{\star \rm POP}$|⁠) and K-band luminosity, central velocity dispersion (σ₀), or M₅₀₀ (not shown here). There are also no clear correlations between the stellar population properties and the stellar velocity anisotropy or velocity dispersion slope (Loubser et al. 2018, 2020).

Lastly, it is worth noting that most (but not all) BCGs with young stellar components have non-constant |$\Upsilon _{\star \rm POP}$|⁠. One of the BCGs with a young stellar population component (Abell 383) has constant |$\Upsilon _{\star \rm POP}$| (i.e. a younger population distributed through the inner and outer apertures).

4 COMPARING |$\Upsilon _{\star \rm POP}$| TO |$\Upsilon _{\star \rm DYN}$| AND CONSTRAINTS ON THE IMF

For a given stellar population age and metallicity, the value of |$\Upsilon _{\star \rm POP}$| depends strongly on the IMF, since an excess of low-mass stars (or stellar remnants) contributes to the mass without significantly contributing to the luminosity. To compare to |$\Upsilon _{\star \rm DYN}$|⁠, we use the derived |$\Upsilon _{\star \rm POP}$| (averaged over both bins, 0–5 and 5–15 kpc), first for a Salpeter IMF (Salpeter 1955), and then for a Kroupa IMF (Kroupa 2001) in the r band (as described in Section 3.2).³ We therefore use two different functional forms of the IMF, namely a single (Salpeter) and double (Kroupa) power law. The Salpeter IMF is equivalent to a single power-law IMF slope of Γ = 1.3, whereas a Kroupa IMF is indistinguishable from a double power-law IMF (two power laws joined by a spline) with a slope of Γ_b = 1.3 above 0.6 M_⊙ and tapered for masses below 0.5 M_⊙. Both IMFs have a lower and upper mass cutoff of 0.1 and 100 M_⊙, respectively.

To determine the |$\Upsilon _{\star \rm POP}$| values for the Kroupa IMF, we refit the stellar population parameters, similar to the analysis in Section 3, assuming a Kroupa IMF. We propagate the 1σ errors from the stellar population parameters to estimate errors on the resulting |$\Upsilon _{\star \rm POP}$|⁠. We discuss possible systematic errors on |$\Upsilon _{\star \rm POP}$| in Section 4.2.

From our dynamical mass models, we use the best-fitting value for the |$\Upsilon _{\star \rm DYN}$| parameter (using the stellar, central, and dark matter mass components ‘⋆ + CEN + DM’, from Loubser et al. 2020; see summary in Section 2), where |$\Upsilon _{\star \rm DYN}$| was kept constant over our fitting range (<20 kpc).⁴ The central mass component in our dynamical mass models represents a supermassive black hole (with mass estimated from the M_nBH–σ relation) as described in Loubser et al. (2020). To estimate the dark matter mass component, we assume an NFW profile (Navarro et al. 1996) and take M_DM and r₂₀₀ from weak lensing observations (Herbonnet et al. 2020). We use M_DM = αM₂₀₀, where α = Ω_M/(Ω_M − Ω_b), and the total mass M₂₀₀ is obtained from weak lensing. We assume the baryon fraction within r₂₀₀ is equal to the cosmological baryon fraction, and therefore M₂₀₀ ∼ 1.2M_DM, as described in detail in Loubser et al. (2020). The dynamical models are free from any assumptions about the stellar populations. We show various robustness tests of our dynamical mass models (e.g. influence of the mass and radius of black hole, influence of the point spread function, and sensitivity of our derived parameters to the dark matter mass distribution and concentration value used for the dark matter halo) in the appendices in Loubser et al. (2020).

We aim to eliminate possible uncertainties in the stellar population and dynamical modelling results by comparing |$\Upsilon _{\star \rm POP}$| to |$\Upsilon _{\star \rm DYN}$| for the BCGs for which we have dynamical models, but eliminate: (i) the two BCGs (Abell 963 and 2055) for which we find the best-fitting β_|$z$| in the ‘⋆ + CEN + DM’ case to indicate extreme tangential anisotropy β_|$z$| < −1 (see Loubser et al. 2020); (ii) the BCGs where we detect young stellar components; and (iii) the BCGs where we detect significant age gradients between the inner and the outer stellar population bin (i.e. |$\Upsilon _{\star \rm POP}$| is non-constant within 15 kpc). Composite stellar populations would increase |$\Upsilon _{\star \rm POP}$| compared to SSPs (Cappellari et al. 2006; Trager, Faber & Dressler 2008), although this is much more relevant for low-mass early-type galaxies. 14 BCGs, indicated with a ‘^⋆’ in Table 2, satisfy these criteria. We plot the comparison between |$\Upsilon _{\star \rm POP}$| for a Salpeter (orange symbols) and Kroupa (blue symbols) IMF, and |$\Upsilon _{\star \rm DYN}$| in Fig. 2. We have included the dark matter mass component, estimated from weak lensing measurements for these clusters, in our best-fitting values for |$\Upsilon _{\star \rm DYN}$|⁠. This decreases |$\Upsilon _{\star \rm DYN}$| on average by 8.3 ± 2.9 per cent over our kinematic range. As expected, a ‘heavier’ IMF shifts the data points towards larger values for |$\Upsilon _{\star \rm POP}$| (x-axis). If we would have added even more dark matter in the mass models, it would push |$\Upsilon _{\star \rm DYN}$| (y-axis) to lower values (see analysis in Loubser et al. 2020).

4.1 Discussion on constraining the IMF

From Fig. 2, we can see that, for a Salpeter IMF, the majority (10) of the BCGs fall above the 1-to-1 line, but there is a subset of four BCGs below the 1-to-1 line (Abell 754, 1689, 1763, and 1942). The BCGs above the line (with |$\Upsilon _{\star \rm DYN} \gt \Upsilon _{\star \rm POP}$|⁠) are better described by the bottom-heavy Salpeter, or an even ‘heavier’, IMF due to a higher fraction of low-mass stars. For relatively old stellar populations (as expected for passively evolving, massive early-type galaxies), only stars with masses below ∼1 M_⊙ are present (see Martín-Navarro et al. 2015b; Lyubenova et al. 2016). This limits stellar populations based IMF analysis of massive early-type galaxies to its low-mass end. As seen from Fig. 2, a Kroupa IMF (which has fewer stars with masses below 0.5M_⊙ than a Salpeter IMF) describes the four BCGs below the 1-to-1 line better, and is closer to reconciling the stellar mass-to-light ratios measured from stellar population and kinematic properties.

Most studies find a bottom-heavy IMF (e.g. Salpeter) for massive early-type galaxies (Cappellari et al. 2013; La Barbera et al. 2015; Martín-Navarro et al. 2015a). However, as discussed in Section 1, the SNELLS galaxies yielded lensing masses in strong disagreement with a bottom-heavy IMF for massive early-type galaxies (Newman et al. 2017; Smith et al. 2017), instead measuring mass-to-light ratios consistent with a Milky Way-like IMF (e.g. Kroupa), both in low-resolution ground-based, and in high-resolution space-based observations (Collier et al. 2018). However, the SNELLS systems also show other peculiar properties as we discuss further in Section 4.3.

We compare the fraction of dark matter mass (to the total mass), as derived from our dynamical models where the dark matter mass was constrained from weak lensing measurements (Herbonnet et al. 2020), with previous studies. We find an average dark matter fraction of ∼8 per cent within 0.38R_e as described in Loubser et al. (2020; see Loubser et al. 2020 for the values for individual clusters). This is consistent with the median dark matter fraction found for massive early-type galaxies in comparable studies (e.g. Cappellari et al. 2013; Lyubenova et al. 2016). Including more mass attributed to dark matter will bring some of the BCGs above the 1-to-1 line closer to the line, but it will move the four BCGs better described by a Kroupa IMF further below the 1-to-1 comparison line. A universal IMF is therefore not only inconsistent with the weak lensing mass measurements, but would also imply high dark matter fractions for some BCGs (in the very centres of the galaxies) and none for others.

4.2 Systematic errors on |$\Upsilon _{\star \rm DYN}$| and |$\Upsilon _{\star \rm POP}$|

4.2.1 |$\Upsilon _{\star \rm DYN}$|

BCGs can be classified as oblate, triaxial, or prolate (Krajnović et al. 2018). We therefore also use the axisymmetric Jeans equations under the assumption of an anisotropic (three-integral) velocity ellipsoid that is aligned with a spherical polar coordinate system (Cappellari 2020). A comparison between the solutions obtained from JAM with spherical polar coordinates and JAM with cylindrical polar coordinates allows for a robust assessment of the best-fitting parameters from the dynamical modelling. For BCGs with decreasing velocity dispersion profiles, |$\Upsilon _{\star \rm DYN}$| (spherical coordinates) is up to ∼15 per cent lower, and for BCGs with rising velocity dispersion profiles, |$\Upsilon _{\star \rm DYN}$| (spherical coordinates) is up to ∼15 per cent higher (Loubser et al. 2020). We indicate this possible uncertainty in Fig. A1 in Appendix  A, and find that it cannot account for the scatter we observe above or below the 1-to-1 line.

The IMF may also vary radially within high-mass early-type galaxies, becoming bottom heavier towards the central regions (van Dokkum et al. 2017; Oldham & Auger 2018; Parikh et al. 2018; Sarzi et al. 2018; La Barbera et al. 2019). For M87, Sarzi et al. (2018) found that the IMF drops from a low-mass excess in the core to a Milky Way IMF at 0.4R_e (for our BCG sample this corresponds to ∼15 kpc, on average). Vaughan et al. (2018) find that the IMF for the BCG in the Fornax cluster NGC 1399 is heavier than the Milky Way and remains constant out to 0.7R_e, before it decreases to become marginally consistent with a Milky Way IMF. Sonnenfeld et al. (2019), for their massive galaxies from the Baryon Oscillation Spectroscopic Survey (BOSS) constant mass (CMASS) sample, find that the region where the IMF is significantly heavier than that of the Milky Way is smaller than the scales probed by the Einstein radius of the lenses in their sample (5–10 kpc).

As a first test, we also investigate the comparison between |$\Upsilon _{\star \rm DYN}$| and |$\Upsilon _{\star \rm POP}$| (for a Salpeter and Kroupa IMF), using |$\Upsilon _{\star \rm POP}$| for just the central 5 kpc of the BCGs (i.e. just the inner aperture, for the same 14 BCGs). We also expect the dark matter mass component to contribute very little to the total mass in the central 5 kpc, so we use |$\Upsilon _{\star \rm DYN}$| values for the (‘⋆ + CEN’) mass models from Loubser et al. (2020), i.e. where a dark matter mass component is not included in the dynamical modelling. We show this in Appendix  B (Fig. B1) and find that it does not change our conclusions.

As an additional test, we estimate a parametrized |$\Upsilon _{\star \rm DYN}$| (to vary as a function of radius) following the results for M87 from Sarzi et al. (2018; in Appendix  C). We estimate the r-band |$\Upsilon _{\star \rm DYN}$| ratio at 2.5 kpc (for the inner aperture 0–5 kpc) to be 50 per cent higher than at 10 kpc (outer aperture of 5–15 kpc). We rerun our dynamical modelling from Loubser et al. (2020), and show our findings in Fig. C1, and find that it does not change our conclusions. Our results strongly suggest that there is substantial scatter in the IMF among the most massive early-type galaxies. It is unlikely that any other systematic overestimation or underestimation of |$\Upsilon _{\star \rm DYN}$| from assumptions in our dynamical models (see e.g. the studies of Thomas et al. 2007a, b; Li et al. 2016) will change this result.

4.2.2 |$\Upsilon _{\star \rm POP}$|

To assess possible systematic errors in our stellar population analysis, and the possible effect on |$\Upsilon _{\star \rm POP}$|⁠, we also use the full spectral fitting code firefly as described in Wilkinson et al. (2017) and applied to 2 million SDSS DR14 and DEEP DR4 spectra in Comparat et al. (2017) and to MaNGA data in Goddard et al. (2017).

firefly is a minimization fitting code that fits combinations of single-burst stellar population models by using an iterative best-fitting process and Bayesian methods. No priors are applied, and all solutions (and their weight) are retained within a statistical cut. Moreover, no multiplicative or additive polynomials are used to adjust the spectral shape (as used in ulyss), and the continuum information is retained and used to determine the parameters. firefly is compared to starlight (Cid Fernandes et al. 2005), steckmap (Ocvirk et al. 2006), and vespa (Tojeiro et al. 2009) in Wilkinson et al. (2017).

We use the MaStar stellar population models (Maraston et al. 2020), built from the MaNGA stellar library (Yan et al. 2019), with the empirical (E-MaStar) stellar library (see Maraston et al. 2020) due to its coverage in age and metallicity parameter space and high resolution. We derive light-weighted SSP-equivalent ages and metallicities for a Salpeter and a Kroupa IMF, similar to our method in ulyss (but with no priors on age components), and we use the stellar population results to derive the |$\Upsilon _{\star \rm POP}$| in the rband.

In Appendix  D (Fig. D1), we illustrate and describe how Fig. 2 changes using a different stellar population model, stellar library, and full spectrum fitting method. The comparison in Appendix  D indicates that realistic errors on |$\Upsilon _{\star \rm POP}$| should be larger to include the systematic errors from using a different combination of stellar population model, library, and fitting method. Even though using a different stellar population analysis has a pronounced effect on the determination of |$\Upsilon _{\star \rm POP}$|⁠, it does not eliminate the variety of IMFs necessary to describe the BCGs. In Fig. D3, we show that the average (and standard deviation) of the two different determinations of |$\Upsilon _{\star \rm POP}$| still scatter above and below the 1-to-1 line. This standard deviation is used as error bars in Fig. 2 to indicate realistic systematic uncertainties on |$\Upsilon _{\star \rm POP}$|⁠. We further briefly compare the ages and metallicities derived using six different combinations of stellar population model, library, and fitting method in Appendix  D, and find that no single combination can consistently derive stellar population parameters that would reconcile the |$\Upsilon _{\star \rm POP}$| above as well as below the 1-to-1 line with |$\Upsilon _{\star \rm DYN}$|⁠.

4.3 Galaxies with low |$\Upsilon _{\star \rm DYN}$|

We now consider the four BCGs below the 1-to-1 comparison line in Fig. 2 individually. The BCG in Abell 1689 has |$\Upsilon _{\star \rm DYN}$| of 3.60 ± 0.09, which is not notably lower than what we expect for passively evolving BCGs, and Fig. D3 shows that if we take the average of two stellar population method/models (ulyss/Vazdekis/MILES and firefly/MaStar/E-MaStar) then the error bars on |$\Upsilon _{\star \rm POP}$| are larger and both the Salpeter and Kroupa IMF reach the 1-to-1 line reconciling the dynamical and stellar population estimates. For Abell 754 (⁠|$\Upsilon _{\star \rm DYN} = 2.56$|⁠), with one of the lowest dynamical mass estimates, and a corresponding low central velocity dispersion (295 ± 14 km s⁻¹), Fig. D3 also shows that the average of two stellar population methods/models reconcile the dynamical and stellar population estimates.

For two BCGs, the dynamical and stellar population estimates of Υ_⋆ could not be reconciled in any way. Abell 1942 (⁠|$\Upsilon _{\star \rm DYN} = 1.23$|⁠) has an SSP-equivalent age of ∼4 Gyr, and therefore had more recent star formation (but not enough or recent enough to constrain the young stellar component, see Loubser et al. 2016). Abell 1942 is also bright (M_K = –27.40 mag, compared to the average in our sample M_K = −26.52 mag), but with a below-average central velocity dispersion (296 km s⁻¹) compared to the other BCGs. For the BCG in Abell 1763 (⁠|$\Upsilon _{\star \rm DYN} = 1.89$|⁠), we also find that it has one of the lowest contributions of dark matter mass to total mass in the centre (3.4 per cent) on account of its brightness (M_K = –27.33 mag) and high stellar mass.

Therefore, the BCGs of Abell 1942 and 1763 are peculiar in that they are very bright (i.e. high stellar mass, and low |$\Upsilon _{\star \rm DYN}$|⁠) compared to other BCGs of similar velocity dispersion. This is similar to the findings of Newman et al. (2017, their fig. 2) for the SNELLS systems, which have a total Υ lower than that for galaxies with similar mass and fall below the expectation from Υ versus velocity dispersion scaling relation (i.e. a projection of the Fundamental Plane). Their results show that the SNELLS galaxies also have peculiar properties that are not related to possible issues with, e.g. the IMF and the stellar population analysis.

As mentioned in Section 1, Smith (2014) concluded that the IMF constraints for individual galaxies, determined using different techniques (⁠|$\Upsilon _{\star \rm POP}$| versus |$\Upsilon _{\star \rm DYN}$|⁠), do not always correlate on a galaxy-by-galaxy basis. However, Lyubenova et al. (2016) used the CALIFA sample to show that stellar populations and stellar dynamics give consistent results for a systematically varying IMF, and their results strongly suggest no case-to-case inconsistencies. They emphasize that inconsistencies found in other studies may be due to differences in stellar population models or aperture sizes used, or non-optimal dark matter halo corrections. In this paper, our direct comparison between |$\Upsilon _{\star \rm POP}$| and |$\Upsilon _{\star \rm DYN}$|⁠, and the search for possible correlations with other properties (Section 5), assumes a consistency between the dynamical and stellar population determined Υ_⋆, and no case-by-case variations.

5 CORRELATIONS OF THE MASS-EXCESS FACTOR WITH OTHER PROPERTIES

When comparing mass-to-light ratios from dynamical mass models to those from stellar populations, the constraint on the IMF is often expressed as the mass-excess factor |$\alpha = \log _{10} (\Upsilon _{\star \rm DYN}) - \log _{10} (\Upsilon _{\star \rm POP})$|⁠, where |$\Upsilon _{\star \rm POP}$| is most-commonly the predicted stellar mass-to-light ratio for some reference IMF, in our case for a Salpeter IMF (Treu et al. 2010). We present the mass-excess factor and the propagated uncertainty in Table 3 (columns 2 and 3) for the 14 BCGs from Fig. 2. If we assume no case-to-case inconsistencies as described above, a mass-excess factor of α = 0 implies the galaxy has a stellar mass-to-light ratio in agreement with the chosen reference IMF. An α > 0 indicates departures from this IMF which could be either a bottom- or top-heavy IMF due to a higher fraction of low-mass stars or stellar remnants, respectively. As mentioned in Section 4, for passively evolving, massive early-type galaxies consisting of old, low-mass stars, it points to a bottom-heavy IMF.

Although it is known that the mass-excess factor (α) increases with velocity dispersion (σ) in early-type galaxies (Cappellari et al. 2013; La Barbera et al. 2015), consensus has not yet been reached on if, and how, it varies with other galaxy properties (Conroy & van Dokkum 2012; McDermid et al. 2014; La Barbera et al. 2015; Martín-Navarro et al. 2015a, b). We plot the α of the 14 BCGs from Table 3 against velocity dispersion in Fig. 3. As the 14 BCGs are all very massive, our range in velocity dispersion is too narrow to measure a correlation with α.

However, we compare our data points to correlations from the literature for early-type galaxies, derived using different methods, following the compilation summarized in Barber, Crain & Schaye (2018, their fig. 1). In Fig. 3, we include, in blue, the observed relation (for r band) between mass excess (also for a Salpeter reference IMF) with velocity dispersion measured in Cappellari et al. (2013) for massive elliptical galaxies in the ATLAS3D survey. Also shown are observed trends from Conroy & van Dokkum (2012), La Barbera et al. (2013), and Spiniello et al. (2014). We also show fits from Li et al. (2017) for elliptical and lenticular galaxies using two different stellar population models and a Salpeter IMF. The dashed blue line shows the extension to the ATLAS3D (Cappellari et al. 2013) sample by including SLACS galaxies (Posacki et al. 2015). This mass-excess–velocity dispersion relation from Posacki et al. (2015) is slightly steeper than from the ATLAS3D sample alone, suggesting a non-linear relation that depends on the range of velocity dispersion probed.

From Fig. 3, we see that the variations in α are larger than any measurement errors and strongly suggests inconsistency with a single universal IMF. We note that our BCGs are generally more massive than the velocity dispersion range covered by previous studies of massive ellipticals. We emphasize that we excluded the BCGs where we see evidence for a radially variant Υ_⋆ in their stellar population properties. We also note that our quantities (σ₀, α) were measured in the central part of the BCGs (which is less affected by the dark matter contribution), whereas the quantities from the reference correlations were measured within the effective radius (R_e). Nevertheless, dark matter mass, as constrained from weak lensing observations, are included in our dynamical models (see Loubser et al. 2020). The velocity dispersions measured for the reference correlations (σ_e within R_e) trace the overall galaxy potential (mass) instead of the detailed kinematics. For the BCGs, we found that the velocity dispersion gradients are very diverse (Loubser et al. 2018), and can be steep within their large effective radii (Loubser et al. 2020). However, we use this information and do an aperture correction to derive the velocity dispersion within the half-light radius, σ_e, for the BCGs. We show this plot in Appendix  E, and find that it does not influence our conclusions. For the BCGs better described by a Salpeter (or heavier) IMF, our data points fall on an extrapolation of the correlations, also suggesting a systematic variation of the IMF for these galaxies in that an increasingly bottom-heavy IMF is needed for more massive galaxies (higher σ), as opposed to case-by-case inconsistencies between the dynamical and stellar population methods. This figure also emphasizes the substantial scatter in the IMF among the most massive galaxies.

For the 14 BCGs from Section 4, we further show (Fig. 4) the mass-excess factor (α) against redshift, K-band luminosity, and M₅₀₀, and find no correlations. Some previous studies e.g. McDermid et al. (2014) and Davis & McDermid (2017) do not find any correlations between their IMF mass-excess parameter and galaxy dynamical or stellar population properties. Other studies find a correlation with metallicity, e.g. Martín-Navarro et al. (2015a, b) and Zhou et al. (2019), or find a correlation with [α/Fe] abundances, e.g. Conroy & van Dokkum (2012), both properties that are strongly correlated with the velocity dispersion. These correlations, or lack thereof, are better studied using large samples of galaxies over a wide mass range (Posacki et al. 2015), or cosmological simulations (Sonnenfeld, Nipoti & Treu 2017; Barber et al. 2018).

In future, it will be interesting to use a physical parametrization instead of the mass-excess factor. An example is F_0.5, defined as the fraction of stars with masses below 0.5 M_⊙, often used for massive early-type galaxies (La Barbera et al. 2013, 2015; Martín-Navarro et al. 2015b, 2019; Lyubenova et al. 2016). However, the interpretation of this is vastly complicated by the fact that spectroscopic IMF studies are sensitive to the present-day stellar populations, and the large variety of parametrizations of IMF variations makes comparison between different methods difficult. Cosmological, hydrodynamical simulations (e.g. Clauwens, Schaye & Franx 2016; Barber et al. 2018) are needed to aid our interpretation of physical parametrizations derived from observed properties.

6 CONCLUSIONS

We investigate the stellar and dynamical mass profiles as well as stellar populations of BCGs, and use the results to place constraints on the stellar IMF. In an accompanying paper (Loubser et al. 2020), we have modelled the stellar and dynamical masses of 25 BCGs using the MGE (Emsellem et al. 1994; Cappellari 2002) and JAM (Cappellari et al. 2006; Cappellari 2008) methods, deriving the stellar mass-to-light ratio (⁠|$\Upsilon _{\star \rm DYN}$|⁠), and stellar velocity anisotropy (β_|$z$|), where the dark matter mass was constrained from weak lensing results. Here, we study the spatially resolved stellar population properties of the BCGs, and use it to calculate their stellar mass-to-light ratios (⁠|$\Upsilon _{\star \rm POP}$|⁠) assuming a (single power-law) Salpeter and (double power-law) Kroupa IMF. We compare the stellar mass-to-light ratios derived from the two independent methods (⁠|$\Upsilon _{\star \rm POP}$| with |$\Upsilon _{\star \rm DYN}$|⁠) and use it to constrain the IMF. We summarize our main conclusions as follows:

• Of the 14 MENeaCS BCGs, we find three BCGs (Abell 780, Abell 1795, and Abell 2055) with very young stellar populations (<200 Myr) in their inner (0–5 kpc) apertures. Together with the 4 of the 18 CCCP BCGs for which we have detected young stellar populations in Loubser et al. (2016), it constitutes 22 per cent of the full sample of 32 BCGs, equally distributed in redshift. From Fig. 1 it follows, similar to our conclusions in Loubser et al. (2016), that the BCGs with young stellar populations are located in host clusters with short cooling times.

• We use the ages and metallicities derived from our spectra to determine |$\Upsilon _{\star \rm POP}$| in the r filter in an inner (0–5 kpc) and outer (5–15 kpc) aperture, and find that 19/32 (60 per cent) of the BCGs have constant |$\Upsilon _{\star \rm POP}$| over this radial range (0–15 kpc). The non-constant |$\Upsilon _{\star \rm POP}$| in these BCGs is primarily driven by age gradients between the inner and outer apertures (i.e. inner aperture significantly younger than the outer aperture).

• To place constraints on the IMF, we eliminate: (i) the two BCGs for which we find extreme tangential anisotropy (β_|$z$| < −1, see Loubser et al. 2020); (ii) the BCGs where we detect young stellar components; and (iii) the BCGs where we detect significant age gradients between the inner and the outer stellar population bin (i.e. |$\Upsilon _{\star \rm POP}$| is non-constant within 15 kpc). We compare |$\Upsilon _{\star \rm POP}$| versus |$\Upsilon _{\star \rm DYN}$| (with the inclusion of a dark matter mass component) for 14 BCGs in Fig. 2.

  From this comparison (Fig. 2), it follows that the majority of the BCGs are better described by a bottom-heavy Salpeter or heavier IMF, but there is a small number of data points below the 1-to-1 line, for which a Kroupa IMF describes the data much better. This agrees with the studies for massive early-type galaxies that find a ‘bottom-heavy’ (Salpeter-like) IMF (Cappellari et al. 2013; La Barbera et al. 2015; Martín-Navarro et al. 2015a) as well as with the SNELLS galaxies which instead measured Υ_⋆ consistent with a Milky Way (Kroupa-like) IMF (Newman et al. 2017; Smith et al. 2017; Collier et al. 2018), and confirms substantial scatter in the IMF among the most massive galaxies.

  We test various possible systematic effects on this direct comparison (e.g. only considering the inner aperture, or using a variable |$\Upsilon _{\star \rm DYN}$| in our dynamical models) in the Appendices. Although using different stellar population analysis models and a different fitting method can bring some data points in better agreement with the 1-to-1 line, none of these systematic effects can consistently reconcile all the |$\Upsilon _{\star \rm POP}$| and |$\Upsilon _{\star \rm DYN}$| measurements above and below the 1-to-1 line.

• We plot the mass-excess factor against velocity dispersion in Fig. 3, and compare it to correlations from the literature for massive elliptical galaxies, derived using different methods (also see Barber et al. 2018). For the BCGs better described by a Salpeter (or heavier) IMF, our data points fall on an extrapolation of the correlations, also suggesting a systematic variation of the IMF for these galaxies (as opposed to a case-by-case inconsistencies).

In summary, we find substantial scatter in the IMF among the most massive galaxies. For most BCGs, a Salpeter (or even more bottom-heavy) IMF is required. For one BCG, a Kroupa IMF is preferred, and for another two BCGs an IMF even lighter than a Kroupa IMF is preferred. Our dark matter fractions are consistent with previous studies (on average), and even though including more mass attributed to dark matter will bring some of the BCGs above the 1-to-1 line closer to the line, it will move the four BCGs better described by a Kroupa IMF further below the 1-to-1 comparison line. A universal IMF will therefore not only be inconsistent with the weak lensing mass measurements, but will also imply very high dark matter fractions for some BCGs (within our limited radial range of <15 kpc) and none for others.

ACKNOWLEDGEMENTS

We gratefully acknowledge the constructive report from the anonymous reviewer. This research was enabled, in part, by support provided by the bilateral funding agreement between the National Research Foundation (NRF) of South Africa, and the Netherlands Organisation for Scientific Research (NWO) to SIL and HH. AB acknowledges support from the Natural Sciences and Engineering Research Council (NSERC), Canada, through the Discovery Grant program. HH acknowledges support from NWO through VICI grant 639.043.512. YMB acknowledges funding from the EU Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement 747645 (ClusterGal) and the NWO through VENI grant 639.041.751.

Based, in part, on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), and Ministério da Ciência, Tecnologia e Inovação (Brazil).

Based, in part, on observations obtained at the CFHT which is operated by the National Research Council of Canada, the Institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with support from the Canadian Space Agency.

Any opinion, finding, and conclusion or recommendation expressed in this material is that of the author(s) and the NRF does not accept any liability in this regard.

DATA AVAILABILITY

Based, in part, on observations obtained at the Gemini Observatory (GN-2008A-Q103, GN-2008B-Q5, GN-2009A-Q107, GN-2009B-Q118, GS-2008A-Q21, GS-2008B-Q4, GS-2009A-Q82). We use versions of the publicly available ulyss (Koleva et al. 2009), and firefly (Wilkinson et al. 2017) software packages, modified for our data and analysis. This research made use of astropy,⁵ a community-developed core python package for Astronomy (Astropy Collaboration 2013, 2018). The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: CYLINDRICALLY OR SPHERICALLY ALIGNED JEANS AXISYMMETRIC MODELS

In addition to the axisymmetric Jeans equations for cylindrically aligned coordinates, we also use the axisymmetric Jeans equations for spherically aligned coordinates (Cappellari 2020), since a comparison between the two solutions allow for a robust assessment of the modelling results and dynamical parameters. We adapt the spherically aligned JAM models (abbreviated as JAM_sph) for our purpose by modifying the models to fit our data, and to include a dark matter mass component.

In Loubser et al. (2020), we found that neither model is significantly, or consistently, better or worse than the other. However, we did find that JAM_sph is relatively insensitive to anisotropy (i.e. a bigger change from β = 0 is required to best fit the observed kinematics). Corresponding to this systematic change in velocity anisotropy in JAM_sph, there is a systematic change in best-fitting |$\Upsilon _{\star \rm DYN}$|⁠, where |$\Upsilon _{\star \rm DYN}$| is lower for decreasing ν_rms profile BCGs (i.e. radial anisotropy, or β positive), and higher for increasing ν_rms profile BCGs (i.e. tangential anisotropy, or β negative). These changes are larger than the statistical error on the parameters, and can correspond to a change of up to 15 per cent in best-fitting |$\Upsilon _{\star \rm DYN}$|⁠. We indicate a 15 per cent uncertainty, and the direction of the change (whether it is higher or lower |$\Upsilon _{\star \rm DYN}$|⁠) in Fig. A1, and find that it will move data points both above and below the 1-to-1 line higher and lower, and that this assumption we made in the dynamical modelling does not affect our conclusions.

APPENDIX B: THE INNER 5 KPC

We also investigate the comparison between |$\Upsilon _{\star \rm DYN}$| and |$\Upsilon _{\star \rm POP}$| (for a Salpeter and Kroupa IMF), using |$\Upsilon _{\star \rm POP}$| for just the central 5 kpc of the BCGs. We expect the dark matter mass component to contribute very little to the total mass, so we use |$\Upsilon _{\star \rm DYN}$| values for the (‘⋆ + CEN’) mass models from Loubser et al. (2020), i.e. where a dark matter mass component is not included in the dynamical modelling. We still find that the majority of the BCGs are better described by a ‘heavier’ IMF, but that four BCGs below the 1-to-1 line are better described by a Kroupa IMF.

APPENDIX C: DYNAMICAL MODELS WITH VARIABLE |$\Upsilon _{\star \rm DYN}$|

As described in Section 4.2, the IMF may also vary radially within high-mass early-type galaxies, becoming bottom heavier towards the central regions (van Dokkum et al. 2017; Oldham & Auger 2018; Parikh et al. 2018; Sarzi et al. 2018; La Barbera et al. 2019). In addition to the test we do in Appendix  B, we estimate a parametrized |$\Upsilon _{\star \rm DYN}$| (to vary as a function of radius) following the results for M87 from Sarzi et al. (2018). Following their results (their fig. 11), we estimate the r-band |$\Upsilon _{\star \rm DYN}$| ratio at 2.5 kpc (for the inner aperture 0–5 kpc) to be 50 per cent higher than at 10 kpc (outer aperture of 5–15 kpc). We note that Vaughan et al. (2018) find a constant IMF up to 0.7R_e for the BCG NGC 1399, which will imply a constant |$\Upsilon _{\star \rm DYN}$| over our radial range. Nevertheless, we rerun our dynamical modelling from Loubser et al. (2020), and instead of adding additional free parameters, we find the best-fitting |$\Upsilon _{\star \rm DYN}$| (along with the best-fitting β_|$z$|) if we assume |$\Upsilon _{\star \rm DYN}$| is 50 per cent larger in the inner aperture than in the outer aperture. We show our findings in Fig. C1.

We find that the new best-fitting inner |$\Upsilon _{\star \rm DYN}$| is higher than the constant overall (constant) best-fitting |$\Upsilon _{\star \rm DYN}$| previously, and the outer |$\Upsilon _{\star \rm DYN}$| is lower than the constant overall best-fitting |$\Upsilon _{\star \rm DYN}$| previously, and that in general the new best-fitting β_|$z$| is slightly lower than previously but it does not influence any conclusions from Loubser et al. (2020). A parametrized, variable |$\Upsilon _{\star \rm DYN}$| does not systematically reconcile |$\Upsilon _{\star \rm DYN}$| with |$\Upsilon _{\star \rm POP}$| for all the BCGs, in the inner or the outer aperture, and for the Salpeter or the Kroupa IMF.

APPENDIX D: |$\Upsilon _{\star \rm POP}$| DERIVED FROM MASTAR STELLAR POPULATION MODELS USING firefly

As described in Section 4.2, we illustrate how Fig. 2 changes using a different stellar population model, stellar library, and full spectrum fitting method. We use firefly) as described in Wilkinson et al. (2017), and the MaStar stellar population models with the empirical (E-MaStar) stellar library (Maraston et al. 2020). We derive light-weighted SSP-equivalent ages and metallicities for a Salpeter and a Kroupa IMF, similar to our method in ulyss (but with no priors on age components), and we use the stellar population results to derive the |$\Upsilon _{\star \rm POP}$| in the r band (Fig. D1).

Of the four BCGs below the 1-to-1 line (Abell 754, 1689, 1763, and 1942), one BCG (Abell 1689) is now on the 1-to-1 line for a Salpeter IMF, and another BCG (Abell 754) now has error bars that also encompass the 1-to-1 line. For the two other BCGs (Abell 1763 and 1942), the Kroupa IMF is still more comparable to the 1-to-1 line, and in addition one other BCG (Abell 1650) is now also below the 1-to-1 line (for a Salpeter IMF). Of the 14 BCGs (using a Salpeter IMF), seven |$\Upsilon _{\star \rm POP}$| determined using ulyss/Vazdekis/MILES and firefly/MaStar/E-MaStar agree within the errors. Of the other seven, four |$\Upsilon _{\star \rm POP}$| are smaller using firefly/MaStar/E-MaStar, and three larger. We also illustrate how Fig. 3 changes using a different stellar population model, stellar library, and full spectrum fitting method in Fig. D2. Two BCGs now fall below the known correlations, but the scatter at higher mass-excess factor (α) is now more pronounced.

This comparison (ulyss/Vazdekis/MILES and firefly/MaStar/E-MaStar) indicates that realistic errors on |$\Upsilon _{\star \rm POP}$| should be larger to include the systematic errors from using a different combination of stellar population model, library, and fitting method. Even though using a different stellar population analysis has a pronounced effect on the determination of |$\Upsilon _{\star \rm POP}$|⁠, it does not eliminate the variety of IMFs necessary to describe the BCGs. In Fig. D3, we show that the average (and standard deviation) of the two different determinations of |$\Upsilon _{\star \rm POP}$| still scatter above and below the 1-to-1 line.

Our goal is not to do a comprehensive comparison between stellar population models, but to understand the effect different models have on determining |$\Upsilon _{\star \rm POP}$|⁠. We further briefly directly compare ages and metallicities derived with six different combinations of fitting methods, stellar population models, and stellar libraries. We use: fitting methods ulyss (Koleva et al. 2009), and firefly (Wilkinson et al. 2017); the stellar population models Vazdekis (Vazdekis et al. 2010), PEGASE-HR (Le Borgne et al. 2004), MaStar (Maraston et al. 2020), and M11 (Maraston & Strömbäck 2011); the stellar libraries MILES (Sánchez-Blázquez et al. 2006), ELODIE 3.1 (Prugniel & Soubiran 2001), and E-MaStar and Th-MaStar (Yan et al. 2019; Maraston et al. 2020).

• ulyss + Vazdekis + MILES

• ulyss + Pegase + ELODIE

• firefly + MaStar + E-MaStar

• firefly + MaStar + Th-MaStar

• firefly + M11 + MILES

• firefly + M11 + ELODIE

We find that the average SSP-equivalent age and metallicity determined for six different combination fall within the average age and metallicity determined using ulyss/Vazdekis/MILES and firefly/MaStar/E-MaStar, and no single combination can consistently derive an age and metallicity that would reconcile the |$\Upsilon _{\star \rm POP}$| above as well as below the 1-to-1 line with |$\Upsilon _{\star \rm DYN}$|⁠.

APPENDIX E: THE VELOCITY DISPERSION WITHIN σ_e

We do an aperture correction using the velocity dispersion profiles measured in Loubser et al. (2018), and R_e measured in Loubser et al. (2020), to derive the velocity dispersion within the half-light radius, σ_e and plot the aperture corrected plot in Appendix  E. However, our α measurement is limited to a 15 kpc aperture. Fig. E1 shows that an aperture correction does not change our overall conclusion of substantial scatter in the IMF for the most massive galaxies.