Recovery of the low- and high-mass end slopes of the IMF in massive early-type galaxies using detailed elemental abundances

Star formation in the early Universe has left its imprint on the chemistry of observable stars in galaxies. We derive elemental abundances and the slope of the low-mass end of the initial mass function (IMF) for a sample of 25 very massive galaxies, separated into brightest cluster galaxies (BCGs) and their massive satellites. The elemental abundances of BGCs and their satellites are similar, but for some elements, satellite galaxies show a correlation with the global velocity dispersion. Using a subset of derived elemental abundances, we model the star formation histories of these galaxies with chemical evolution models, and predict the high-mass end slope of the IMF and star formation timescales. The high-mass end IMF slope of the satellite galaxies correlates with the global velocity dispersion. The low- and the high-mass end IMF slopes are weakly correlated in a general sense that top heavy IMFs are paired with bottom heavy IMFs. Our results do not necessarily imply that the IMF was simultaneously bottom and top heavy. Instead, our findings can be considered consistent with a temporal variation in the IMF, where, for massive galaxies, the high-mass end IMF slope is representative of the very early age and the low-mass end slope of the later star formation. The small but noticeable differences between the BCGs and the satellites in terms of their elemental abundances and IMF slopes, together with their stellar kinematical properties, suggest somewhat different formation pathways, where BCGs experience more major, gas-free mergers.


INTRODUCTION
The elemental abundances in the atmospheres of stars are indicative for the chemical composition of the gas from which these stars formed.Different chemical elements can be related to different types of stellar end-products, which generally correspond to different time scales.Understanding the origin of these elements therefore helps interpreting the evolution of galaxies (e.g.Tinsley 1979).Type II, or core collapse (CC) supernovae (SNe) produce primarily  elements (e.g.Na, Mg, Si, Ca).Type Ia supernovae mainly contribute iron peak elements (e.g.Fe, Cr, Co, Ni, Cu, Mn).The third important channel for the nucleosynthesis of light elements is the formation of these elements in asymptotic giant branch (AGB) stars.Depending on the initial mass of these stars, they can be important contributors ★ E-mail: dkrajnovic@aip.de of C or N and many of the s-process elements e.g.Sr, Y, Nb, Ba (e.g.Nomoto et al. 2013).
The time scales at which the ISM is enriched by these processes vary: CC SNe, such as type Ib, type Ic and type II SNe, happen almost instantly (∼few Myr) after star formation starts.Type Ia SNe happen on much longer time scales (Maoz & Badenes 2010), probably with a delay time distribution that has a slope declining somewhat steeper than  −1 .Observationally the delay time distribution has however been difficult to constrain.Theoretically, delay time distributions can be predicted from binary population synthesis calculations.Depending on whether the progenitor of the SN Ia consists of a single degenerate (SD) or a double degenerate (DD) object, different time scales are predicted.For example, Mennekens et al. (2010) quote 50 Myr as earliest supernova for the DD mechanism, whereas their fig.4 shows a ∼ 150 − 200 Myr delay for the onset of supernovae through the SD channel.
AGB enrichment sets in after the heaviest progenitors evolve away from the main sequence (∼ 8M ⊙ , corresponding to ∼ 40 Myr).
The high equivalent widths of CO and CN indices have been used to argue extended star formation histories (e.g.Sánchez-Blázquez et al. 2003;Mármol-Queraltó et al. 2009), as higher equivalent widths can result from enrichment by AGB stars.
The observed amount of  elements has been widely used as a measure of star formation duration.This is usually done using Mg as a proxy and normalising it by the solar abundance ([Mg/H]), with respect to the amount of iron peak elements ([Fe/H]).The  abundance is then written as the solar abundance-normalised logarithm [Mg/Fe].It is now well established that there exists a relation between  abundance and the mass (or velocity dispersion) of early-type galaxies (e.g.Worthey et al. 1992;Davies et al. 1993;Carollo & Danziger 1994;Maraston et al. 2003;Thomas et al. 2003a).Over the past 15 years it has however become clear that abundances and ages are not the only stellar population parameters that vary with galaxy mass; more and more evidence is pointing in the direction of a non-universal initial mass function (IMF; e.g.Hopkins 2018, for a review).
The IMF in the Milky Way was inferred for the first time by Salpeter (1955) using bright stars (  < 3.5).The Salpeter (1955) IMF is described by a power law with a negative slope of  = 2.35 1 .Later work, in particular the papers by Miller & Scalo (1979) and Kroupa (2001), changed the idea of a single power law describing both the high-mass end and the low-mass end of the IMF in the Milky Way.The high-mass end of the IMF has long been suspected to be top heavy (i.e.having a power law with index shallower than Salpeter's  = 2.35) at early times, based on several lines of evidence (see e.g.Larson 1998;Calura & Menci 2009;De Masi et al. 2019).Observations of nearby star forming galaxies (Gunawardhana et al. 2011) as well as observations of galaxies at high redshift (Nanayakkara et al. 2017(Nanayakkara et al. , 2020) ) confirm that the high-mass end depends on star formation rate.
The analysis of old stellar populations suggests that in the centres of early type galaxies, the slope of the IMF with respect to mass is steeper than the IMF of the Milky Way, i.e., they show an excess of low-mass stars (bottom heavy), with the most bottom heavy IMFs occurring in the more massive galaxies.The evidence for this is based independently on mass-to-light ratios obtained from stellar population modelling and Jeans models (e.g.Cappellari et al. 2012Cappellari et al. , 2013)), from lensing studies (e.g.Treu et al. 2010;Spiniello et al. 2011;Oldham & Auger 2018), from stellar population analysis of dwarf sensitive features in spectra (e.g.van Dokkum & Conroy 2010;Spiniello et al. 2012;La Barbera et al. 2013;Gu et al. 2021;Feldmeier-Krause et al. 2021), but also from alternative methods that count the amount of dwarf stars (e.g.van Dokkum & Conroy 2021).It is, however, not fully understood what is the main driver of the variable low-mass end of the IMF.The IMF slope is known to correlate with velocity dispersion and [/Fe] abundance (Cappellari et al. 2013;McDermid et al. 2014).It is also known that IMF appears to be most bottom-heavy in the centre of galaxies and bottomlighter at larger radii (Martín-Navarro et al. 2015;van Dokkum et al. 2017;Sarzi et al. 2018;Parikh et al. 2018;Barbosa et al. 2020;Feldmeier-Krause et al. 2021), and it is therefore not clear how to reconcile this with the rather flat [Mg/Fe] gradients observed in early-type galaxies (Mehlert et al. 2003;Sánchez-Blázquez et al. 2007;Kuntschner et al. 2010;Greene et al. 2015;Krajnović et al. 2020;Parikh et al. 2021;Feldmeier-Krause et al. 2021).We note 1 Throughout the paper we will use the convention where the IMF is given by the formula  () =    ∼  −  , as in Bastian et al. (2010).
that although there exists evidence from multiple methods for a variable IMF, each of these methods carries a considerable amount of uncertainty with it.Both lensing and stellar dynamics studies have to deal with the problem that they have to separate any non-stellar mass contribution from a stellar mass contribution.In addition, non-homology in elliptical galaxies may cause variation in mass-tolight ratio with increasing galaxy mass that is not always captured by the dynamical models at the required accuracy.For example, axisymmetric models can severely underestimate the mass-to-light ratio of triaxial galaxies (e.g.Thomas et al. 2007;den Brok et al. 2021) and also the assumed shape of the velocity ellipsoid may bias the stellar mass-to-light ratio (Thater et al. 2022).Clauwens et al. (2015) argued that the IMF trends seen in the ATLAS 3D sample may in fact arise from a universal IMF combined with underestimated modelling errors.The simultaneous analysis of mass-to-light ratio excesses from SSP modelling and from kinematics by Smith (2014) showed no evidence for any correlation between the two methods; a similar result was found for the lensing/SSP analysis by Newman et al. (2017).The lensing analysis of 23 strong lenses by Sonnenfeld et al. (2019) shows that on scales of 5-10 kpc an IMF heavier than in the Milky Way is ruled out.On somewhat smaller scale (< 5 kpc), Collier et al. (2018a), using strong lensing, find no evidence for a large mass excess in two lensed galaxies; Collier et al. (2018b) show that this is even the case for a lensed galaxy for which the mass is constrained within  < 1.5kpc.From a theoretical perspective, there are several reasons to expect an environmental dependence for the IMF.Padoan & Nordlund (2002), Hennebelle & Chabrier (2008) and Hopkins (2012) derive an IMF based on the turbulent Jeans mass of molecular clouds.Nam et al. (2021) show that flatter turbulent spectra produce top-heavier IMFs.Weidner et al. (2013a) proposed a temporal variation of the IMF, which changes during the formation of early-type galaxies from a short duration top-heavy IMF to a long duration bottom heavy IMF, motivated by the increased pressure in the interstellar medium from the high amount of high-mass stars formed because of the top-heavy IMF.Lacey et al. (2016) show that semi-analytic galaxy formation models do indeed fit better to data when they allow for a somewhat top-heavy IMF during starburst phases.Barber et al. (2018) implement a pressure dependent IMF in the EAGLE simulations, which leads to weak trends of IMF slope with global metallicity (Barber et al. 2019a) and negative radial IMF gradients (Barber et al. 2019b).
Contrary to the direct measurement of the low-mass slope of the IMF in early-type galaxies (ETGs), which can be determined with advanced stellar population models directly from the observed spectra, the determination of the slope of the high-mass end of the IMF can almost never be done directly, as high-mass stars have lifetimes that are much shorter than the typical age of the stellar populations in elliptical galaxies.Through their latest evolutionary stages, the high-mass stars have, however, left an imprint on the chemistry of the ISM, which is still observable in the low-mass stars formed out of this gas.By analyzing the abundances of the low-mass stars and using chemical evolution models, it is, therefore, still possible to get an inference on the high-mass IMF slope, even after the high-mass stars have died.Moreover, this inference can be directly compared to the low-mass slope of the IMF.
In this paper we derive abundances and stellar population parameters for a sample of 25 very massive ETGs (stellar mass  ⪆ 10 12 M ⊙ ), which we observed with the Multi Unit Spectroscopic Explorer (MUSE, Bacon et al. 2010) at the Very Large Telescope, and derive slopes for the high and low-mass end of the IMF.The sample selection and data reduction are presented in Section 2.
The methods used for the analysis of the data and the chemical modelling are presented in Section 3. The resulting measurement of the elemental abundances, predictions of the chemical evolution modelling, and the derivation of the low-and high-mass end IMF parameters are presented in Section 4. We discuss the uncertainties on the measurements and the chemical evolution models, as well as the implication of the results in Section 5, presenting conclusions in Section 6.

SAMPLE AND DATA
The data used in this work come from the MUSE observations of massive galaxies, collected in the M3G sample (Krajnović et al. 2018).Both the sample selection as well as the data reduction of the M3G galaxies has been presented in Krajnović et al. (2018) and will be elaborated in more details in Krajnović et al (in prep.).We briefly repeat that the sample consists of 25 very massive early-type galaxies found in some of the densest environments.Galaxies are located in central region of rich clusters, with 14 galaxies spread over the three main clusters of the Shapley Super Cluster (SSC), including the three brightest cluster galaxies (BCGs).The remaining 11 galaxies were chosen to be BCGs from a sample of galaxy clusters with richness above 40 (Abell et al. 1989).The sample, therefore consists of 14 BCGs (3 in SSC and 11 in various clusters) and 11 non-BCGs, which we will refer to as satellites (all in SSC).Among satellites there are also second brightest galaxies in the three SSC clusters covered by the survey.Fig. 1 presents mass -size relation for the M3G galaxies, divided into BCGs and satellites, together with galaxies from the MASSIVE Survey (Ma et al. 2014) and the more massive galaxies from the ATLAS 3D Survey (Cappellari et al. 2011) sample.Noteworthy for this work is that, while there is some overlap, M3G BCGs are typically more massive and larger than the M3G satellites.
Galaxies were observed with MUSE, with exposure times between 2 to 6 hr.All data were reduced using v1.2 to v1.6 of the MUSE pipeline (Weilbacher et al. 2020), except PGC 046832, for which the reduction has been presented in den Brok et al. (2021) and for which the data reduction included the use of the Zurich Atmospheric Purge (ZAP, Soto et al. 2016).Basic reduction included bias subtraction, flat-fielding, wavelength calibration, flux calibration, sky subtraction (using offset fields), alignment and drizzling on a common reference frame.The major difference in this work with respect to the steps presented in Krajnović et al. (2018) is that we switched of the telluric correction in the pipeline, but instead determined telluric absorption with the molecfit tool (v1.5.9;Smette et al. 2015), as in Pagotto et al. (2021).We do this by selecting foreground field stars, depending on availability, with a as featureless spectrum as possible, in each (non telluric subtracted) reduced and combined cube.For a few galaxies that have no available field stars, we use instead the centre of the galaxy as the input spectrum for molecfit to calculate a telluric absorption spectrum.We apply the telluric spectra produced by molecfit for correcting the final data cubes.This approach helped in reducing the residual telluric features, but we still track their location and mask them when necessary.
The M3G sample can be put in perspective by comparison with the ATLAS 3D and the MASSIVE samples (Fig. 1).The former is a sample of 260 nearby (<40Mpc) early-type galaxies, less massive than a few ×10 11  ⊙ , of which only one is a BCG (M49 in Virgo, but there is also M87 as the second brightest galaxy), while the sample galaxies are in typically less dense environments compared to the Mass -size relation for M3G galaxies (circles), divided into BCGs (orange) and satellites (light blue).Squares are galaxies from the MASSIVE survey (values taken from Ma et al. 2014), and pentagons show massive part of the ATLAS 3D Survey of nearby early-type galaxies (values taken from Cappellari et al. 2013).Sizes of galaxies from all surveys are based on the 2MASS observations.The masses of the MASSIVE galaxies are based on the 2MASS K-band absolute magnitude and calibration eq. ( 2) from Cappellari (2013).ATLAS 3D masses come from dynamical models.M3G masses are based on the virial estimate and taken from Krajnović et al. (2018).The dashed lines are the lines of the constant velocity dispersion, based on the Cappellari et al. (2006) parametrisation of the virial mass estimator.
clusters that host M3G galaxies.The latter is a sample of the most massive galaxies within 108 Mpc.The M3G galaxies are larger and more massive than galaxies in the ATLAS 3D sample, while there is an overlap with galaxy properties of the MASSIVE sample.M3G galaxies are also more distant, with a mean distance of about 190 Mpc.Therefore, M3G can be considered as an extension of these nearby galaxy samples, both in mass and size.

Abundances within an effective radius
We bin spectra for each galaxy within an elliptical aperture of a semi-major axis radius equal to the half-light (effective) radius of the galaxy.The shape of the elliptical aperture (position angle, ellipticity and size) is based on values from Table 1 of Krajnović et al. (2018).The same aperture was used by Krajnović et al. (2018) to extract the global stellar velocity dispersion,   , which we use throughout this work as a characteristic property of each galaxy.The combined central spectra have average S/N per pixel above 200 in all but one case where it is just under 100.The mean S/N is 235, while in three cases it is higher than 300.Although the data have a high S/N, the use of separate sky observations, and the data reduction methods for sky removal, nevertheless lead to imperfect sky emission-line subtraction.We solve this by using an outer-bin spectrum obtained by combining the spectra outside the central effective radius (elliptical) aperture.We first perform a stellar population fit for each galaxy on this outer-bin, and use it to estimate the residual sky-lines.These sky-lines are then removed from the central effective-radius bin before the final stellar population fit. Figure A1 shows a typical example of the residual sky spectrum and the effect of its removal from the galaxy spectra.See Appendix A for more details.
Two commonly used approaches for extracting information about the properties of stellar populations from galaxy spectra are the index fitting (e.g., Worthey 1994;Trager et al. 2000;Thomas et al. 2003b;Kuntschner et al. 2006;Schiavon 2007) and the full spectral fitting (e.g.Cappellari & Emsellem 2004;Cappellari 2017;Cid Fernandes et al. 2005;Ocvirk et al. 2006;Tojeiro et al. 2007;Koleva et al. 2009;Wilkinson et al. 2017) methods, while in the recent literature there are also hybrid methods combining spectra and photometric (broad or narrow band) magnitudes (Cappellari 2023), or defining specific indices but fitting each pixels within the index wavelength ranges (Martín-Navarro et al. 2019).All methods have their advantages, and in some cases limiting the range of the pixel fitting methods is beneficial (e.g.Vazdekis et al. 2015;Gonçalves et al. 2020), while in some cases index fitting might be a more appropriate method (e.g.Zibetti et al. 2024).Nevertheless, we select the full spectral fitting as it does not depend on the specific definition of the spectral features and their pseudo-continua, and benefits from the higher S/N of our spectra.Barbosa et al. (2021) present an educational overview of recent choices and an extensive references, advocating the full spectral fitting for estimating elemental abundances and the IMF of galaxies.A comparison of our results, and those from index fitting and full pixel fitting studies is provided in Appendix B and Fig B3.
We model the central effective-radius bin of each galaxy with alf (Conroy & van Dokkum 2012a,b;Conroy et al. 2018), which fits synthesized stellar population models to galaxy spectra.The underlying models are based on MIST isochrones (Choi et al. 2016), MILES stellar spectra (Sánchez-Blázquez et al. 2006a), IRTF stellar spectra (Villaume et al. 2017) and theoretical response curves to change the abundance of individual elements (see references in Conroy et al. 2018).We use the 'full' mode of alf.This means that besides metallicity, 19 elemental abundances are left as free parameters.In addition to a dominant old population, two parameters (age and light fraction) govern the presence of a second (subdominant) younger population with the same elemental abundances.We explore single power law IMFs with a slope fixed to the Salpeter slope ( = 2.35), as well as with a free power-law slope, both with a lower mass cut-off at 0.08 M ⊙ .Although it is possible to run alf with a more complex shape for the IMF, we chose not to do this for the reason of simplicity, but also because recent observations prefer uni-modal power laws over more complex shapes (e.g.Conroy et al. 2017;Zhou et al. 2019).Therefore, we fit for 23 free parameters (age, metallicity, light fraction of the younger component, 19 elements and the IMF slope).alf uses an implementation of the MCMC sampler emcee (Foreman-Mackey et al. 2013).We use nwalkers=1024 and nburn=20000, for the number of walkers used by the MCMC sampler and the duration of the burn-in phase.
For all galaxies we fit the spectral range between 4600 − 8800 Å.Despite the care we take in subtracting the sky lines and correcting the telluric absorption, several spectra still show not fully or over-subtracted sky emission, or over or under-corrected telluric absorption.We therefore mask features that are poorly subtracted (e.g. the [O I] emission at 6300Å) and additionally mask any remaining strong features by hand.We also mask H, the surrounding [N II] lines and the [S II] ∼ 150Å redward of H, all with a width of 800 km/s on each side of the systemic velocity of the galaxy.A few of the M3G galaxies show line emission from ionized gas (Pagotto et al. 2021).For these galaxies, we also masked the H, NaD, [O I] 6300Å and [O III] 4960,5008Å lines, with the same window width.
We show an example fit in Fig. 2, and for all other galaxies in Fig. D1.Clearly the model spectra reproduces the observed data very well.Some features are however visible: i.e. positive residuals are visible around ∼7630Å and ∼7680Å.We see this in almost all fits, and it may be due to the absorption feature around 7665Å being underestimated.The same holds for the feature at ∼7416Å.
Distribution of elemental abundances as a function of global velocity dispersion   (measured within the same central elliptical aperture) can be seen in Fig. B1.The abundance values for each galaxy are presented in Table B.

Chemo-dynamical modelling
As the second step in this work, we use the derived element abundances to constrain the chemical evolution of each galaxy using the code Chempy2 (Rybizki et al. 2017).In our default setup, which is very similar to the one in Rybizki et al. (2017), stars are formed according to a Γ(, ) distribution (eq.( 2) in Rybizki et al. 2017), with a shape parameter  and time scale .In this set up chempy is a single-zone open box chemical model.Gas with primordial abundances is added to a corona with the same rate as star formation.Corona is defined as a reservoir of a well mixed gas surrounding the zone that we are modelling chemically.It is slowly enriched by the stellar feedback and diluted by the inflow of primordial gas, and, therefore, it can be visualised as a circum-galactic medium surrounding our galaxies.Feedback from stars is distributed over the central star forming zone and the corona according to an outflow feedback fraction parameter  out .The gas mass in the corona at the beginning of the simulation starts out with a mass  cor times the Chempy divides a simulation up in discrete time steps.Since the stellar populations in our galaxies are mostly old, we limit the total simulation time length to 3.5 Gyr.The chosen simulation time length is tuned to reach the peak of the star formation epoch of today massive disks and elliptical galaxies at z ∼ 2 (Madau & Dickinson 2014).As visible from the top left panel of Fig. B1, this covers the span of stellar ages for most of our galaxies, except for a few youngest (see discussion at the end of Section 4.1).We also run a chemical model of 7 Gyr duration to test for possible effects of extended star formation (see Appendix C).We use 200 time steps, leading to a resolution of 17.5 Myr.At each time step, stars are formed according to the assumed star formation rate, with a powerlaw IMF with a slope  and (M min , M max ) = (0.08, 100) M ⊙ .At the end of its lifetime, every formed star enriches the ISM according to its mass, metallicity and explosion mechanism.Stellar lifetimes are taken from Argast et al. (2000).In our default setup, stars above 8 M ⊙ explode as core-collapse supernovae (CC-SNe).This ignores the fact that stars with masses between 8 and 10 M ⊙ are expected to explode as electron capture supernovae (EC-SN).These EC-SN will mainly contribute to iron peak elements.Yields for the CC-SN explosions, as well as for hypernova, are taken from Nomoto et al. (2013), original presented by Kobayashi et al. (2011).Above 25 M ⊙ , a fraction    of CC-SNe explodes as hypernovae.This fraction is a free parameter of the model.Stars with masses between 1 and 8M ⊙ evolve as AGB stars.We use the yields of Karakas & Lugaro (2016) for AGB enrichment.Yields for SN Ia are from Seitenzahl et al. (2013).
As we will see in Section 4.2, it is difficult to reproduce all elemental abundances for our galaxies, regardless of the input parameters for chempy (we have investigated several different set up, but present below only the most relevant).We, therefore, limit ourselves to reproducing (fitting) only selected abundances, [O/H], [Fe/H], [Mg/H] and [Na/H], in our chemo-dynamical run with chempy.

Stellar abundances
We use alf results to estimate the average abundances of elements in M3G sample galaxies.For each elemental abundance (Fig. B1), we median combine the values and estimate the 16th and 84th percentile limits of the distributions.We do this for all galaxies in the sample as well as separately for BCGs and satellites, and present the values in Table 1.We show in Fig. 3 the abundances of the BCGs and satellite galaxies with respect to the iron abundance, normalised to the abundance values of the Sun (Asplund et al. 2009).There are no significant difference between satellite galaxies and BCGs, which perhaps can be understood as a consequence of the small difference in mass between these two sets.The most notable (but still statistically insignificant) difference is observed for [Co/Fe], while for [Si/Fe], [K/Fe] and [Mn/Fe] there are BCGs with large elemental abundance values, which drive the upper limits of the distributions.In some cases, it is the satellite galaxies that show large spreads of elemental abundances, such as for [N/Fe], [Na/Fe] and [V/Fe].Nevertheless, the lack of significant differences points to a broadly similar formation scenarios for stars in BCGs and massive satellites at these mass scales, but we will examine in more detail the differences for individual cases below.
A more detailed look at Fig. 3 reveals high average values for some elemental abundances, especially [O/Fe] and [Na/Fe] (for more information about the trends for individual galaxies see also Figs. 4 and B1).In general, the -elements have relatively high values, as expected, as the contribution of CC-SNe are much more pronounced for galaxies with short formation time scales.We note that the -element abundances differ from the ones measured by Conroy et al. (2014), for their galaxies with largest velocity dispersions (their bins with  = 250 and 300km/s).Specifically, [Mg/Fe] is higher by up to ∼ 0.1 dex, while [Si/Fe] is significantly lower than in Conroy et al. (2014) (∼ 0.1 dex), but in the latter case there are galaxies in the M3G sample which reach values observed in the SDSS sample of galaxies by Conroy et al. (2014).Beverage et al. (2023) updated the Conroy et al. (2014) abundances for SDSS galaxies using the latest version of alf (we used the same version, but in a somewhat different set up, especially with respect to IMF parametrisation).The differences between our and their abundances are smaller, but with the same trend as for the Conroy et al. (2014) values.Given the higher value of [Mg/Fe], the values of [Si/Fe] between 0-0.16 found in this work seem somewhat low.It is unclear if this is related to observational issues, for example, because of the limited wavelength range of the data (i.e., the strongest response of Si is observed between 4100-4300Å) for this element, or if the Si abundance is indeed low.The Ca abundance is supposed to be close in value to the iron abundance, and this is indeed the case for our data.
We observe two peaks around the iron peak elements.The Co peak was also observed by Conroy et al. (2014), who also noted a large spread in values as a function of velocity dispersion, with the highest  bin galaxies having larger values than what we observe.As mentioned above, in our case this element shows the largest difference between the BCG and satellite galaxies.Contrary to Conroy et al. (2014), we also observe a peak in the vanadium abundance, with the values much larger than those reported by Conroy et al. (2014), specifically for galaxies with high velocity dispersions.This discrepancy seems to decrease with an updated values for [V/Fe] by Beverage et al. (2023), which are now positive for the most massive galaxies in the SDSS sample, and are within the distribution of values of our (satellite) galaxies.Another difference are the lower values for [Mn/Fe] compared to the Conroy et al. (2014), which in our case are mostly negative, except in a few BCGs, while the average of the Conroy et al. (2014) galaxies are positive at all velocity dispersions.
Abundance measurements with respect to iron as a function of the global velocity dispersion of each galaxy, again normalised by the solar values, can be used as a proxy to investigate trends with galaxy mass.We show them for M3G sample galaxies in Fig. 4 (left column), but only for elements selected because they show correlations between the parameters, while all elements are shown in Fig. B1.We estimate the Pearson correlation coefficient between the velocity dispersion, for which we ignore the small observational uncertainties, and the elemental abundance (taking into account the uncertainties returned by alf).Furthermore, we estimate a 95% confidence interval on the correlation coefficient by repeatedly drawing [X/Fe] values from the alf MCMC chains and bootstrapping different galaxies.The correlation coefficient uncertainty is determined by taking the lower 2.5% and upper 97.5% limit values of all bootstraps, therefore providing a 2 level of significance.
Most elements show little evolution with velocity dispersion.Exceptions to this are [Fe/H], [N/Fe], [Na/Fe] and [Ni/Fe], which we find to correlate, and [V/Fe] which instead anti-correlates with velocity dispersion.For reasons that we do not understand, some of the elements slightly less massive than Fe (Ti, V, Cr, Mn) show negative trends with velocity dispersion.The same was reported for V and Cr by Conroy et al. (2014).On the other hand, somewhat more massive elements (Co, Ni, Cu) show instead positive trends with velocity dispersion.These trends are however dominated by large spread and strong outliers.We find that the correlations between the abundance and the velocity dispersion are the most significant for Na and N in our data, but we note that we have a small range of the velocity dispersion.Nevertheless, the strong correlation trend in [Na/Fe] is also seen in the MASSIVE sample, but MASSIVE galaxies seem to have somewhat smaller [O/Fe] and [Mg/Fe] element abundances (Gu et al. 2021).A more detailed comparison is presented in Appendix B.
In some cases, there is a notable difference between the BCGs and satellites in terms of how strongly their elemental abundances correlate with the velocity dispersion.The largest differences are seen for [Fe/H], [Na/Fe], [Mg/Fe] and [V/Fe].In the case of [Fe/H] abundances, BCGs do not show any evidence for a correlation ( BCG = 0.09 +0.35  −0.48 ), but the satellite galaxies have a high Pearson correlation coefficient:  SAT = 0.82 +0.10  −0.11 .The situation is similar for the [Mg/Fe] abundances, where the observed weak correlation between [Mg/Fe] and   for all galaxies is due to the dilution of a relatively significant correlation for the satellite galaxies ( = 0.44 +0.24  −0.24 ) with non-correlating BCGs ( = 0.08 +0.28 −0.24 ).For .BCGs are denoted with orange points, satellite galaxies with light blue points.We plot here only elements for which a significant correlations, either for BCGs or satellite galaxies.The Pearson correlation coefficient  (for the full sample) is indicated on each panel, with uncertainties indicating the 2 level.The correlations for all other elements are in Figs.B1 and B2.For comparison we include data from the MASSIVE sample (Gu et al. 2022) with small green symbols.the [Na/Fe] abundance the situation is somewhat different as both BCGs and satellite galaxies show a correlation, but the correlation is not statistically significant for BCGs ( = 0.41 +0.26  −0.49 ), but it is for satellites ( = −0.81+0.10  −0.13 ).The situation is reversed for [V/Fe], where the satellites do not show a correlation ( = −0.18+0.33  −0.26 ), while BCGs have a significant one ( = −0.75+0.35  −0.14 ).These results point to a tentative difference in the chemical evolution of BCGs and satellite galaxies in our sample, as Na and Mg are -elements produced in CC-SN explosions, while V is produced both in CC-SN and SN Ia explosions.We summarise most significant correlations in Table B1.
In our nominal alf extraction we left the IMF slope as a free parameter, and in Fig. 5 we show the derived IMF slopes for our galaxies.As the currently alive stars in our massive galaxies are all low mass stars, we will refer to this slope as the low-mass IMF slope,  LM .Numerous studies, both based on the spectral analysis and dynamical arguments (e.g.Conroy  2012a; Cappellari et al. 2012;La Barbera et al. 2013;Ferreras et al. 2013;Spiniello et al. 2014), showed that the low-mass IMF slope correlates with velocity dispersion.Our velocity dispersion range is somewhat small, and the data in Fig. 5 do not show a clear correlation.The Pearson correlation coefficient for the full sample is  = 0.27 +0.30 −0.37 , and  = 0.24 +0.28 −0.36 and  = 0.32 +0.29 −0.40 for BGCs and satellites, respectively.While the correlation seem to be somewhat stronger than that reported for the MASSIVE sample (Gu et al. 2022) 3 , it is similarly not statistically significant, and the lack of correlation can be traced to the short   baseline.Furthermore, 8/25 (32%) galaxies have  LM larger than the Salpeter value, while 7/25 (28%) galaxies have  LM smaller than 2, significantly deviating from the Salpeter slope.This is in agreement with previous findings, e.g.where approximately one-third of the galaxies with  > 200 km/s in Conroy & van Dokkum (2012a) show evidence for an IMF with a steeper  LM (i.e. more bottom heavy) than Salpeter, while only about 1/4 of galaxies in the MASSIVE sample ( > 200 km/s) have the low-mass IMF slope below the Salpeter (Gu et al. 2022).We also note that the large spread of  LM values at the mean   for our sample ( ∼ 290 km/s) is also consistent with results for lensing studies (e.g.Newman et al. 2017;Collier et al. 2018a,b), which typically find lower (low-mass end) IMF slopes.
There is no difference in the  LM between BCGs and satellite galaxies, with average values of 2.1±0.3 and 2.1±0.4,respectively, even though 5 satellite galaxies and 3 BGCs have larger than Salpeter IMF slope.The low-mass IMF slope derived by alf, as well as other IMF slopes used in this work are presented in Table 2.
We show the individual correlations between elements and    for selected elements in Fig. 4 (right column) and for all elements in Fig. B2.As expected (Conroy & van Dokkum 2012b), we see a correlation between [Mg/Fe] and  LM , and we also find statistically significant correlations for [Sr/Fe] with a Pearson correlation coefficient of  = 0.71 +0.17 −0.33 .For these elements the correlations exist for the full sample and when separating into BCGs and satellite galaxies.In a few cases, the BCG galaxies (and not the satellites) have a noticeable trend between the element abundances and the low-mass end IMF slope ([N/Fe], [Ti/Fe], [Co/Fe], and [Ba/Fe]), but they are not statistically significant (i.e. in case of [Ti/Fe]  = 0.45 +0.14  −0.53 , or [Ba/Fe]  = 0.60 +0.01 −0.70 ).The most significant correlations are summarised in Table B1.
The population modelling with alf provided also the luminosity-weighted stellar ages for our galaxies.Most populations are very old, reaching the limits of the models and the age of the Universe (for a discussion why this is possible, see McDermid et al. 2015), but there are six galaxies with stellar ages lower than 10 Gyr.A particular case is PGC073000, the youngest galaxy in the sample with an age of 5.9 ± 0.3 Gyr.Stellar populations of this ages are just about old enough for the derivation of the low-mass end IMF slope based on spectral fitting (Conroy & van Dokkum 2012a).Nevertheless, this galaxy also has the highest value of  LM (2.77), while its overall metallicity is close to solar.Elemental abundances, in particular for elements created in CC SN, such as [Na/Fe] or [Mg/Fe], are moderate to high.The relatively low overall metallicity, low stellar age and the existence of emission-line gas in the nucleus, mass of which is the highest for M3G sample (Pagotto et al. 2021), suggest this galaxy experienced a rejuvenation event, and a significant secondary star formation.The high  LM is likely not related to the actual measured elemental abundance, except to possibly Na values, as we show in Section 5.1.Importantly, this galaxy is the only one in the sample with a clear evidence for multiple star-formation phases, and, in this respect, the simple set up of our chempy models is likely not applicable for this galaxy.This is not necessarily because the duration of the chemical evolution model is set to 3.5 Gyr (as we show in Appendix C), but because our model does not account for accretion beyond the corona, or any kind of mergers.While PGC073000 warrants a more detailed stellar population and chemo-dynamical modelling, we nevertheless keep the same chempy set up for all galaxies, and highlight PGC073000 in the text when necessary.

Chemical evolution models
In Fig. C1 we show the results of the parameters from the chemical evolution model using the Nomoto et al. (2013) yields for corecollapse SNe.The values are listed in Table C1.As expected, all galaxies show short star formation time scales, with profile shapes that are close to exponential, and decreasing with time.The models prefer no or little outflows to the corona.For most galaxies we find a fraction of CC-SN over all other supernovae (i.e.hypernovae plus other supernovae) close to 1, but we note that fixing this parameter to either 1.0 or 0.5 makes little difference for the results.The star formation efficiency is found to be close to one and similar for all galaxies.Most of other parameters are poorly constrained with large uncertainties, and there is little difference between BCGs and satellite galaxies.
In the same Appendix C and Fig. C1, we show the results of an extended duration chempy run, where the model set-up was kept the same as for the nominal run, except that the simulation time was extended to 7 Gyr.All model parameters were kept free, as before.There are no significant differences between the models, and we continue using the results of the 3.5 Gyr run chempy in the rest of this work.Furthermore we tested chemical models with respect to different outflow fraction, but there were no significant changes to the main results.We refer the discussion of the similarities between the chemical models to the appendix.
One parameter deserves some attention: the slope of the high- mass IMF ( HM ), which shows an anti-correlation with velocity dispersion.We show it in more details in Fig. 6.While the Pearson correlation coefficient for the full sample is  = −0.55+0.31 −0.19 and is statistically significant, BCGs do not show a significant correlation ( = −0.30+0.58  −0.39 ).It is the satellite galaxies that drive the anti-correlation of the high-mass IMF slope with velocity dispersion for the full sample ( = −0.84+0.12  −0.09 ).In the chempy model with an extended duration, as well with different levels of outflow fractions, the resulting  HM are consistent with the one presented here, conserving the same anti-correlations (see Appendix C).
The origin of this anti-correlation could be related to the noticeable correlations of the [Na/Fe] and [Fe/H] with velocity dispersion (Fig. B1).Specifically, as mentioned above, the satellite galaxies in our sample show a clear rise in the metallicity ([Fe/H]) with   , but also show a stronger correlation for [Na/Fe] with   compared to BCGs.Furthermore, Nomoto et al. (2013) show that [Na/Fe] abundance, weighted by the initial mass function and their yields, increases with the metallicity much more than any other element (their fig.10).Therefore, the chemical evolution models are trying to compensate for the increased metallicity and strong correlations in [Na/Fe] of the satellite galaxies, while trying to fit a slight increase in [Mg/Fe] (for both populations).
We note that the  2 values of the Chempy fits are well above one.One reason for this is the small value of the uncertainties on the abundances, which are of order ∼ 0.01 (Fig. B1), but also show the general inability of the chemical model to reproduce the "observed" element abundance obtained with alf.In the Section 5.2 we will comment more on these issues.
An inference of the slope of the high-mass IMF for individual galaxies obtained with chemical modelling allows us a simultaneous comparison of the low and high-mass end slopes, shown in Fig. 7.We note that the high-mass end of the IMF is always flatter than Salpeter regardless of the chemical evolution model set up; per definition the high-mass end is therefore always top heavy.BCGs and satellite galaxies seem to have somewhat different distributions  of the low and high-mass IMF slopes.This can be seen by calculating the Pearson correlation coefficient for the two samples.For BCGs,  = −0.09+0.36 −0.29 at the 95% confidence level clearly shows no correlation between the IMF slopes, while for the satellite galaxies  = −0.32+0.46  −0.28 suggest a statistically significant correlation, but not a particularly strong one.Combining the two subsamples results in a low and statistically insignificant Pearson correlation coefficient ( = −0.23 +0.41  −0.32 ) as one could judge from the figure by eye.Therefore, for the satellite galaxies there is a weak trend within their central effective radius that, in general, if a galaxy has a more bottom heavy IMF, it also had a more top heavy IMF; a galaxy that has a bottom light IMF, also had a top light IMF.Such a correlation does not seem to exist for BCGs, which for a relatively uniform high-mass IMF slope (centred on  HM = 2.1) exhibit a range of low-mass IMF slopes, from Kroupa-to super-Salpeter-like (1.4 <  LM < 2.7).Krajnović et al. (2018) showed that almost half of M3G galaxies have rotation around their major axis, what is sometimes described as prolate-like rotation and can be quantified with kinematic misalignment angle (Franx et al. 1991), the angle measured between the major axis and the orientation of the velocity field 4 .Prolate-like rotation in M3G galaxies was defined when the kinematic misalign-4 Alternatively, the angle between the minor axis and the projection of the angular mometnum vector.ment angle is greater than 75 • .Such rotation is mostly present in M3G BCGs, but some satellites also have it.In order to test if the kinematics will bring more information, we show in Fig. 7 galaxies with prolate-like rotation with pentagons.Galaxies with prolate-like rotation are distributed differently from both the BCGs and satellites, suggesting a positive correlation (top-heavy to bottom-heavy IMFs).However, the correlation is not significant ( = 0.17 +0.29  −0.35 ).On the other hand, the Pearson correlation coefficient for galaxies without prolate-rotation is  = −0.36+0.36  −0.26 , somewhat stronger and more significant compared to that of the satellites, but the differences are small.Conroy & van Dokkum (2012b) were the first to find a correlation between [Mg/Fe] and IMF slope, and they argued that the low-mass IMF slope could be a result of the short formation time scale of galaxies with high [Mg/Fe].Such as correlation is also found in our sample (Fig. B2).The Chempy models on M3G galaxies, however, suggest the star formation time scales are similar for galaxies with different velocity dispersions (Fig. C1, panel a), while the chemical evolution model compensates the differences in the elemental abundances by pairing different low-and high-mass end IMF slopes.The general abundance values, however, drive the general trend of pairing the bottom heavy (low-mass) with the top heavy (high-mass) IMFs.

Uncertainties in the values of the abundances and the low-mass end IMF slope
We found that both the average abundances of N and Na (Fig. 4), as well as, to a lesser extent, the low-mass end slope of the IMF (Fig. 5), change with increasing velocity dispersion.The possibilities exist therefore that an actual change in the IMF slope with velocity dispersion induces a non-existing trend in elemental abundance measured via absorption features, especially for those that are also sensitive to IMF (e.g.La Barbera et al. 2013;Spiniello et al. 2014).For example, are the trends for N or Na abundances due to shortcomings in the stellar population models, or, in the opposite way, that trends in abundances manifest itself as a change in IMF slope?To asses the first possibility, we have re-measured the abundances from the spectra with alf, but fixed the IMF slope to the Salpeter value ( = 2.35).The difference in the abundances between the free-IMF and fixed-IMF measurements are generally very small (the mean difference is typically less than 0.01 dex).We show the cases of [N/Fe] and [Na/Fe] in Fig. 8.It is thus unlikely that the free low-mass IMF slope changes the measured abundances in such a way that they imply a strong co-variance between the low-mass end and the high-mass end slopes.One interesting change in the measured parameters is, however, seen for the average age of the old component.In the first panel of Fig. 8, we show that this age is almost always lower for those measurements for which the IMF slope was fixed to the Salpeter value.Fixing the dwarf-to-giant ratio in stellar models also fixes the ratio of turn-off stars to giant stars.Therefore, imposing a Salpeter slope will give a higher number of main sequence stars and deeper hydrogen absorption lines, which the stellar models perhaps compensate for by decreasing the age of the population.As the low-mass IMF slope in our sample, when left free, is often shallower than Salpeter, the result is that for many galaxies the ages of old populations decrease somewhat.Nevertheless, it is surprising that the average abundances stay however roughly constant when varying the IMF slope.Additionally, galaxies for which we masked H, because of emission-line contamination, do not show evidence for their lowmass IMF slopes being outliers.
While relative precision is more relevant in this workl than the absolute accuracy of the derived values, we note that the statistical errors on element abundances are quite small (Table B).They range from 0.01 dex for [Fe/H], [Mg/Fe] or [Na/Fe], to 0.1-0.3 for abundances of elements with minor spectral features in our wavelength range (i.e.[Cu/Fe], [Sr/Fe] or [Eu/Fe]).Such small errors are nevertheless not uncommon in similar studies of massive galaxies.Gu et al. (2022) present uncertainties of ∼ 0.02 for both [Fe/H] and [Mg/Fe] for galaxies with  > 250 km/s also using alf.Working with the SDSS spectra, and limiting again to galaxies with  > 250 km/s, Walcher et al. (2015) derive uncertainties of 0.03 for both metallicity and -element abundances.While these two studies used full spectral fitting methods, Loubser et al. (2009) perform index fitting on a sample of BCGs (long-slit spectra) and obtain average uncertainties of 0.07 and 0.04 for the metallicity and element abundances, respectively.Kuntschner et al. (2010), working with IFS data and index fitting methods report errors of 0.02-0.04 and 0.02-0.05for metallicity and element abundances, respectively, of early-type galaxies.The S/N of our data are likely higher than in most of the mentioned studies, while our aperture of one effective radius is also larger, both of which contribute to somewhat smaller errors in our case.We present a more detailed comparison with literature values in Appendix B.
Such small statistical errors suggest that the actual error budget could be dominated by the systematic errors.Those are hard to measure, as they depend on a number of assumptions inherent to the modelling approach.The quality of the stellar population models, fitted wavelength range, assumptions on IMF and the general number of free parameters used in the fit.Regarding alf, some of these aspects were tested by Conroy et al. (2014, their figure 17), who conclude that one can expect a typical change in abundance patters of the order of 0.05 dex or less, while some cases can be as large as 0.1.We report similar tests in the lower panels of Fig. 8, where we fixed the slope of the IMF (as in the top panels), masked all Na lines within our spectral range and limited the fit to < 6900Å.Our conclusions are consistent with Conroy et al. (2014), where the systematic uncertainties range between 0.05 and 0.1 dex.This is highlighted by [Fe/H], [Mg/Fe], [O/Fe] and [Na/Fe], where the differences stay below or around 0.05 dex, except in the case of shorter fitting range, when they are closer to 0.1 dex.Note that these are the elemental abundances that were used to constrain the chempy fits.We further discuss the difference in the derived alf parameters in Appendix B.
The second case, in which abundances cause the low-mass part of the IMF to change, is not as easily invalidated.We do however note that one of the strongest indicators for the low-mass IMF slope is the Na line at 8190Å.The sensitivity of this line to the proportions of cool dwarfs and giants has been known since Spinrad & Taylor (1971), and has subsequently been used by several authors (e.g.Carter et al. 1986;Faber & French 1980;van Dokkum & Conroy 2010;Spiniello et al. 2012).To show that the low-mass IMF slope is not a reflection of the Na abundance, we re-fit the low-mass part of the IMF from the stellar spectra, but with the 8190Å Na line masked.We also re-run alf with all Na line features masked (including also the Na D doublet at 5892Å), even though various authors (e.g.La Barbera et

<690 nm
Figure 8. Top: The difference between the selected alf results allowing for the low-mass end IMF slope ( LM ) to be free or fixed to the Salpeter value ( = 2.3), respectively.From left to right: the difference in the log age, [N/Fe] and [Na/Fe] abundances (IMF free -IMF fixed).Orange symbols are for BCGs and light blue for satellite galaxies.Symbols with red edge on the left-most panel are galaxies with emission-line gas where we masked the H line in alf runs, as indicated on the legend.Note the small scatter in elemental abundance and relatively large scatter in age, where the age of the models with free IMF slope is systematically larger.Bottom: for each galaxy (indicated by their   values), we plot all free alf parameters, as the differences between the nominal fit (outlined in Section 3.1) values and three alf setups, which are, from left to right: IMF slope fixed to the Salpeter values (as in the top three panels), the full set of Na lines masked, and limiting the extraction to the "blue" region ( < 6900Å).Red, blue, green and brown lines highlight the differences in values of [Fe/H], [Mg/Fe], [O/Fe] and [Na/Fe], respectively.The single vertical bar in each panel highlights an error of ±0.05 dex.
this absorption features, although sensitive to IMF slope, is more sensitive to abundance.
The results are shown in Fig. 9, where only small differences are noticeable between the low-mass IMF slopes with partially or fully masked Na lines.Both Na masks result in a consistent lowmass IMF slope with respect to the one from the unmasked run.The largest discrepancy for both masks is observed for PGC073000, which has the youngest stellar populations in the sample.The lowmass end IMF slope measured without Na lines is close to the Salpeter value or lower.The reason for this change could be that the unusually high value is driven by the combination of Na abundance (which is not extreme) and young age, or due to a difficulties in fitting the line, due to emission.The panels in Fig. D1 pertaining to PGC073000 (bottom left), show a possible contribution of HeI 5875 emission-line.Similarly, the IMF slope is influenced in a few other galaxies with emission-lines (i.e.PGC047273, Pagotto et al. 2021), and we masked fully or partially the regions around their NaD doublets (for the nominal fit).Overall however, Fig. 9 results suggest that the observed trends are genuine and not driven by internal inconsistencies.

Uncertainties in the chemical evolution modelling
As we noted in Section 4.2, our Chempy fits to the elemental abundances produced by alf have large  2 values.We visualise how well the chemical evolution models can reproduce observed data for our galaxies in Fig. 10.We first fitted the Chempy resulting parameters (presented in Fig. C1) versus the velocity dispersion with straight lines and used these models to predict elemental abundances as functions of the velocity dispersion.In Fig. 10 we show the model predictions (as lines) and the alf abundances for the following elements: C, N, O, Na, Mg, Fe, Si, Ca and Ti.Note, however, that chempy models are only constrained by [O/Fe], [Na/Fe], [Mg/Fe] and [Fe/H] (see Section 3.2).Part of the reason of leaving out most elemental abundances for our chemical evolution modelling was to focus on elements for which the yields are more reliable.In general the uncertainties increase for heavier elements, which provides an explanation why we are not able to reproduce the abundances of Ti and Si.
The chemical models reasonably reproduce the elemental abundances used to directly constrain them, although [Na/Fe] and [Mg/Fe] are typically over or under predicted, respectively.When it comes to pure prediction of elemental abundances (those not used to constrain the models), the models typically do not fit them particularly well (with a possible exception of [Ca/Fe]).Using yields from Chieffi & Limongi (2004) or Limongi & Chieffi (2018), as a replacement for Nomoto et al. (2013), does not improve the predictions in general, with some possible exceptions (i.e.[C/Fe] for Chieffi & Limongi (2004) or [Ca/Fe] for Limongi & Chieffi (2018)).An interesting case is nitrogen, which seems to be bracketed between chemical models based on the Nomoto et al. (2013) and Chieffi & Limongi (2004) yields, and those of Limongi & Chieffi (2018).Nitrogen can be produced in multiple ways, originating both in AGB stars and in CC-SN (Kobayashi et al. 2020).It can be a by-product of the CNO cycle in stars, but could also be produced from initial hydrogen and helium present in a star.In the first case, the nitro-  The relation between the low-mass end IMF slopes derived using alf by masking Na lines (-axis) and not masking (-axis).We used two different masks, only applying it to the Na 8190Å line, or both on 8190Å line and NaD doublet.Only small differences are seen between two masks and there is an overall good agreement with the unmasked values.
gen is called secondary, as the production of nitrogen is dependent on the initial CNO abundance (Clayton 1983;Arnett 1996).This leads to a metallicity dependence on the nitrogen abundance.If on the other hand nitrogen is produced in a primary way, it would, to first order, follow the abundances of carbon and oxygen (Talbot & Arnett 1974).One way to produce sufficient nitrogen is therefore to include stellar rotation (Meynet & Maeder 2002a,b).Since the yields of Nomoto et al. (2013) do not include rotation this might be the reason why our models under-predict it.However, even without rotation (the  = 0 km/s models of set R) the Limongi & Chieffi (2018) yields already produce more than sufficient N, and are also able to reproduce the increase seen with velocity dispersion, as can be seen from Fig. 10.
The Limongi & Chieffi (2018) yields also over-predict the [C/Fe] abundance, possibly because a lot of carbon is ejected before entering the Wolf-Rayet phase for stars heavier than 10 M ⊙ .C is produced either in low mass (1-4M ⊙ ) AGB stars or in massive (> 10M ⊙ ) stars.As the high-mass C yields depend very sensitively on the treatment of convective overshooting, predictions for C yields vary strongly.Additionally, C yields may be particularly sensitive to binary evolution (Farmer et al. 2021).Fig. 10 also shows that it is difficult to obtain models with [Mg/Fe] high enough to explain the observations, especially so with Limongi & Chieffi (2018) yields.Difficulties with the underproduction of Mg in the Limongi & Chieffi (2018) models were previously also noted by Prantzos et al. (2018), who attribute them to presently poorly modelled physical phenomena.Keeping in mind the level at which the data can be reproduced, we continue the discussion about the only parameter of the chemical evolution model of our sample of very massive galaxies that shows a clear variability with the galaxy velocity dispersion, the high-mass IMF slope  HM (Fig. 6).
Our preferred yield sets of Nomoto et al. (2013) result in a slight trend of  HM with  LM for the satellite galaxies (Fig. 7).This trend relies largely on the abundances of Na and  elements.In Fig. 11, we show the two IMF slopes when we, instead of Nomoto et al. (2013) use the yields of Chieffi & Limongi (2004).The high-mass IMF slopes using the Chieffi & Limongi (2004) yields ( HM(CL04) ) are on average smaller (mean slope is 2.02) compared to the  HM based on the Nomoto et al. ( 2013) yields (mean slope 2.1).One galaxy, previously highlighted PGC073000, seems to be an outlier with  HM(CL04) ∼ 2.182.If we neglect this galaxy, and a systematic shift to lower values for  HM(CL04) , the distribution of low-mass and high-mass IMF slopes is very similar to that seen on Fig. 7.The Pearson correlation coefficient in Fig. 11 for the satellite galaxies is somewhat smaller, but still statistically significant with  = −0.28+0.36  −0.25 , while the coefficient for the BCGs is still not significant ( = −0.25 +0.32  −0.32 ).We conclude that using different yields does not change the results from our chemical evolution modelling, and that there might exist a tentative anti-correlation between lowand high-mass end IMF slopes, specifically among massive non-BCG galaxies.More objects are, however, needed to establish the real trends.
We also used the yields of Limongi & Chieffi (2018) to explore the high-mass end of the IMF.Leaving [Mg/Fe] out of the analysis (as it is poorly reproduced) and only using [O/Fe], [Fe/H] and [Na/Fe], we also find a trend between the high-mass end IMF slope and velocity dispersion, which is, however, even steeper and leads to more extreme values of  ∼ 1 for the galaxies with the highest velocity dispersions.We do not consider these models any further.
Recently, Cinquegrana & Karakas (2022) published yields for AGB stars at very high metallicities.Particularly at metallicities Z>0.04, these yields predict that AGB stars can produce very high amounts of 14 N and 23 Na, which could potentially explain the trends we see with velocity dispersion and might remove the need for topheavy IMFs at the highest dispersions.To test this scenario we have incorporated these yields in our distribution of Chempy, and supplemented them with the Karakas & Lugaro (2016) AGB yields at lower metallicities.In practice this does not seem to change the results of the chemical evolution models, as metallicities almost never reach values above Z = 0.04 and instead points at Na being mainly produced in massive stars, and returned via CC-SN (i.e.Kobayashi et al. 2020).The N abundance in the M3G galaxies correlates weakly with the Ba abundance.Although the latter is an s process element, suggesting a possible origin in AGB stars, we note for our galaxies N abundance does not correlate with Sr abundance (another element mostly produced in AGB stars), while (e.g.Limongi & Chieffi 2018) show that Ba can also be produced in (rotating) massive stars.Fig. C1 presents the set of free parameters using in the chemical modelling, for both the "short" (3.5 Gyr) and the "long" (7 Gyr) runs.It is telling that while the IMF slopes of the satellite galaxies in both cases show the consistent anti-correlation with the velocity dispersions, while this is not the case for any other parameter.We stress that all chempy parameters were kept free during modelling, and the resulting values are mostly flat across the range of the   covered by our sample, with a possible exception of the SN Ia delay times, which in the "long" chempy run show a weak trend with   .Furthermore, we run chempy models by setting the outflow fraction to zero and fixing it to x  = 0.2.This means that all enriched material is kept within the ISM, or a fixed fraction is forced to be diluted in the the corona, respectively.While some model parameters show minor changes, the high-mass end IMF slopes are essentially the same and keep the same correlations with the velocity dispersion as in our nominal run (Appendix C).These tests offer further support that  HM is not degenerate with respect to the chemical evolution parameters.
Our high-mass end of the IMF slope is thus mainly dependent on the almost constant [O/Fe] and [Mg/Fe] abundance and increasing [Na/Fe] with velocity dispersion.The decrease of the IMF slope with velocity dispersion for satellite galaxies is then not surprising, as [Na/Fe] strongly increase with stellar mass.We stress that this assumes that we can fully rely on the assumptions of the chemical evolution code (one zone, instant mixing, accretion of gas with primordial abundances) as well as all the uncertainties associated with the yields, and our specific parametrisation of the IMF in the first place.

Differences between BCGs and satellites
The differences in elemental abundances between BCGs and satellite galaxies seem to be minor and at a low statistical significance (Fig. 3).Nevertheless, we see three notable trends: (i) At fixed velocity dispersion, the mean abundances of elements in BCGs are equal or somewhat larger to those in satellite galaxies.There are some exceptions, such as the [Fe/H] and the [Ca/Fe] abundances which are larger for the satellite galaxies.On the other hand, for some elements produced in CC-SN explosions ([Si/Fe] and [K/Fe]), the values for BCGa are higher at fixed velocity dispersion.
(ii) In a number of BCGs, elements primarily produced in SN Ia explosions have larger abundances than in satellite galaxies.One can see this based clearly on the values of [Mn/Fe] and [Co/Fe].
In the case of [V/Fe], while both satellite galaxies and BGCs have similar mean and high values, satellite galaxies are more likely to show lower values.
(iii) Correlations between the elemental abundance and velocity dispersion are stronger for satellite galaxies, especially so for the overall metallicity ([Fe/H]) and the CC-SN produced [Na/Fe] and [Mg/Fe].
These trends can be explained as a combination of the differences in the galaxy mass, where BCGs are somewhat more massive than non-BCGs in our sample (Fig. 1), but also as a consequence of the specific mass assembly paths for these two galaxy types.While all galaxies have high [/Fe] values, indicating short starbursts (at early times), BCGs do not show dependance of the [/Fe], or the metallicity, on the velocity dispersion (or, therefore, galaxy mass).We note that [Fe/H] also correlates with the I-band surface brightness in the aperture of the MUSE data used for this analysis, and may thus be a result of either an averaging of the radial [Fe/H] profile due to merging or a prolonged star formation leading to a higher central surface brightness.We also note that the [Ca/H] abundance does not show any difference between BCGs and satellite galaxies, and The tentative differences between BCGs and satellite galaxies are consistent with the hypothesis that BCGs have more violent formation paths, experiencing more major mergers and accreting stars from different types of galaxies.Our data probe the central half-light regions of galaxies, regions where the episodes of star formation are the most likely, but also the regions influenced the most by similar mass mergers.Our chemical evolution model showed that the star formation occurred long time ago and it was short.However, BGCs have both the -element abundances and Fe-peak elements larger than satellite galaxies, suggesting somewhat more extended star formation histories, or perhaps a different number of bursts of star formation.Similar conclusions were reached by other studies (e.g.Thomas et al. 2010;La Barbera et al. 2014;Gallazzi et al. 2021) using much larger samples of galaxies, over a larger mass range and different levels of star formation.One of the results from these studies was that central galaxies can appear somewhat younger and with lower [/Fe] abundances than their satellites.Nevertheless, we do not consider this inconsistent with our findings due to the very different types of studied galaxies.The above-mentioned literature studies have a much larger leverage on the environment type, spanning poor groups to large clusters, while all our satellites are in similarly rich and massive clusters of galaxies, where they are typically 2nd, 3rd or 4th brightest galaxy (Krajnović et al. in prep).In that sense, and given the low number of galaxies, it is remarkable that we detect differences between M3G galaxies.Finally, the satellites in the above-mentioned studies have significantly lower masses compared to the centrals, while in our case, even though the masses of BCGs are somewhat larger than of satellites, they are generally comparable (and > 10 12 M ⊙ ).
Major mergers, when occurring between galaxies without gas, change   , or, generally, the mass -size relation, only mildly (Boylan- Kolchin et al. 2006;Bezanson et al. 2009;Hopkins et al. 2009;Naab et al. 2009).They also flatten any existing trends in metallicity or element abundances (e.g.Kobayashi 2004;Di Matteo et al. 2009).Such an "equalising" effect in terms of chemical properties can be expected between galaxies of the same type (i.e.BCGs) even if they have somewhat different masses (or   ).The differences that we observe between the satellites and BCGs suggest that, while they have similar masses and effective velocity dispersions, the satellite galaxies experienced less gas-free major mergers.This is also supported by the kinematics of M3G galaxies, as about 50% of BCGs show prolate-like rotation around the major-axis (Krajnović et al. 2018).On the other hand, only about 30% of satellite galaxies show the same type of rotation (Krajnović et al. 2018).Such prolate-like rotation is typically expected to originate in major, gas-free mergers (e.g.Ebrová & Łokas 2017;Li et al. 2018), suggesting a clear formation way for some galaxies.
The IMF differences between BGC and satellite galaxies are small, but are, nevertheless, noticeable.There is a remarkable trend of the high-mass end IMF slope with   for satellite galaxies (Fig. 6), which does not exist for BCGs, and can be traced back to the differences in elemental abundances.Noteworthy is also the distribution of satellite galaxies in  LM − HM plane (Fig. 7), where, again, satellite galaxies show an anti-correlation (at 2 level), while BCGs do not show a statistically significant correlation.Marking the galaxies with their kinematics (with prolate-like and without prolate-like rotation) enhances this difference somewhat: galaxies without prolate-like rotation behave similarly like satellites.Prolatelike rotation is a direct evidence that the last merger was gas-free and a major merger, but galaxies without prolate-like rotation could have also experienced similar type mergers, "diluting" the signal.Nevertheless, the trend suggest that the observed difference originate in different evolution paths taken by BCGs and satellite galaxies.

The shape of the IMF in the most massive galaxies
In Fig. 12 we illustrate the shape of the IMF recovered through alf and chempy models.We combine the low-mass and the high-mass slopes assuming the transition occurs at 0.5 M ⊙ .As already visible in Fig. 7, the spread in the slopes values is generally small, while it is a bit larger on the low-mass end.Both BCGs and satellite galaxies are consistent with bottom heavy low-mass IMF and a top heavy high-mass IMF, overabundant in low and high mass stars compared to the standard IMF parametrisations.This is equivalent to what Martín-Navarro (2016) invoked as a non-canonical shape of the IMF providing both the high [Mg/Fe] values and being overabundant in low mass stars (relative to Salpeter IMF) in massive ellipticals (e.g.Conroy & van Dokkum 2012b;Spiniello et al. 2012;Cappellari et al. 2012;Ferreras et al. 2013).
An alternative way of looking at Fig. 12 is that it does not show a static picture of the IMF in any of these galaxies, but could be considered as a depiction of a time-varying IMF, where the stellar mass (the x-axis) is a proxy for time: massive stars, which do not exist anymore, signify the past, while low mass stars, still burning hydrogen, relate to the present.Here the initial IMF in progenitors of our massive galaxies had an excess of high-mass stars that were born in a short (∼ 0.3 Gyr), but strong starburst, characterised by CC-SN and -element enrichment of the ISM and the subsequent stellar populations.This phase is followed by a prolonged (∼ 1 Gyr), less violent star-burst characterised with a bottom heavy IMF (Vazdekis et al. 1997;Weidner et al. 2013b;Ferreras et al. 2015).
Our results are broadly consistent with the notion that the IMF is not universal or invariant, but should be considered within a specific time and a spatial volume.What we observe is a "galaxy A composite distribution of the initial mass functions in the  () =    ∼  −  form for M3G galaxies.The IMFs are constructed using the low-mass (< 0.5 M ⊙ ) and high-mass (> 0.5 M ⊙ ) end slopes obtained from the alf analysis of the MUSE spectra and the chempy modelling of recovered elemental abundances.For comparison, we plot also the Kroupa (dashed line) and the Salpeter (solid line) IMFs.Orange lines are BCGs and light blue lines are the satellite galaxies of the M3G sample.Vertical dashed-dotted line approximately shows the boundary between stars that still exits in M3G galaxies, and stars that ended their evolutionary paths.Note that this plot is not meant to represent a fixed IMF.Instead it shows a temporal variation of the IMF, in the sense of the IGIMF theory (e.g.Jeřábková et al. 2018), where, for the massive galaxies, the high-mass end is representative of the very early age, and the low-mass end of the later star formation.
IMF" (gIMF) characteristic for the full galaxy, and different from the much more local "stellar IMF" (sIMF; Hopkins 2018).In this context, the "galaxy wide IMF" 5 (Jeřábková et al. 2018), arising from an "integrated galaxy-wide IMF" (IGIMF) model (Kroupa & Weidner 2003;Weidner & Kroupa 2005;Kroupa et al. 2013;Recchi & Kroupa 2015;Yan et al. 2017), is obtained as a summation of many sIMF.The gIMF found in the centres of present day massive galaxies starts with an IMF overabundant in high mass stars (top heavy), and then evolves into an IMF with an excess of low mass stars (bottom heavy), which are also enriched in metals.Crucial point is that the gIMF is not simultaneously bottom heavy and top heavy, but the observed excesses of the high and low mass stars happen sequentially.This evolution of the IMF is present in both galaxy types of our sample (BCGs and satellites), but the differences in the correlation of the IMF slope (specifically  HM vs.   ) suggest that their subsequent evolutionary paths are nevertheless sufficiently different, most likely in terms of the number and types of mergers.
Given that all galaxies in our sample end up as massive, quiescent and giant ellipticals, a broadly speaking argument would be as follows.During the early formation, properties of these galaxy (e.g.gas accretion, start formation, the transition between top to bottom heavy IMF, etc) depend on their total mass and fuel availability (i.e.environment).At some point, some of these galaxies (i.e.BCGs), however, end up being in the centre of the gravitational wells and experience more mergers of all kind, with an increasing fraction of gas-free mergers as time passes.A consequence of this "privileged" position is the accretion and mixing of stars formed in different 5 Different authors have different names, and we will use gIMF.
types of galaxies and under different conditions, maintaining the chemistry, elemental abundances and metallicity of the stars, but removing the signatures of the dependence with galaxy velocity dispersion (or mass).Massive galaxies outside of this "privileged" position are able to retain more of the "fossil" records from the epoch of early star-formation, recognisable in the dependancies with the velocity dispersion (i.e. HM vs.   ).

SUMMARY
We present in this work an analysis of the stellar populations in the M3G sample of 25 very massive galaxies observed with VLT/MUSE instrument.Our goal is to reconstruct the chemical evolution of these galaxies and place constraints on both the low-and high-mass IMF slopes, as an overview of the change in the star-formation.We focus on the central half-light regions of our galaxies, and spectra of very high signal-to-noise ratios.
Using advanced stellar population tool alf, we measure 19 elemental abundances, age, and low-mass end slope of the IMF of the galaxies.The elemental abundances are typically high, as expected for massive galaxies, and are broadly consistent with previous observations based on SDSS galaxies, especially when compared with the measurements obtained with the same version of the tool.Next to higher average values of some elements (e.g.[O/Fe], [Mg/Fe] and [V/Fe]), we find significant correlations between the elemental abundances and the global velocity dispersion.Notable correlations are found for [Fe/H], [N/Fe], [Na/Fe] and [Ni/Fe], while [V/Fe] shows an anti-correlation.The correlations are typically more pronounced (significant at 2 level) for satellite galaxies for [Fe/H] and [Mg/Fe], while in the case of [V/Fe] it is the BCGs that show a significant correlation.
Our stellar population models are constrained with a unimodal IMF which is free to vary between galaxies.We find that in most cases the slope of this low-mass end IMF is similar to the Salpeter, and often larger, but there are at least two galaxies which have the IMF slope consistent with that of (low-mass end) Kroupa IMF.There is no clear difference in the derived slopes between BCGs and satellites.This results is consistent with previous works, both in the sense that very massive galaxies (with high   ) have super-Salpeter IMF slopes and that there are massive galaxies with smaller, Kroupa-like IMF slopes.
We use the derived elemental abundances of selected elements ([Fe/H], [O/Fe], [Mg/Fe] and [Na/Fe]) to constrain the chemical evolution using the Chempy model.While these elements are typically well reproduced by the set of standard yields, other chemical abundance are poorly reproduced.Changing the yields does not necessarily improve the reconstruction of elemental abundances.The chemical evolution modelling results show that the very massive galaxies have similar formation time scales and star formation efficiency.
One parameter critical for the chemical evolution shows an interesting change between galaxies: the high-mass end slope of the IMF.Satellite galaxies show a clear correlation between this IMF and the velocity dispersion, while this is not true for BCGs.Having both the low-and high-mass end IMF slopes, we combine them and show that there is a general tendency that the more top heavy IMF slope (high-mass end) is paired with a more bottom heavy (lowmass end) IMF slope.This correlation is more pronounced for the satellite galaxies, but only at a 1 confidence level.Furthermore, the changes in [/Fe] and [Na/Fe] abundances can be explained by this change in slope of the high-mass end of the IMF with the global velocity dispersion.We test to what extent are the measurement of both the low-and high-mass end IMF slope susceptible to systematic errors in the models, by masking different absorption-lines, changing the wavelength range used for fitting, changing the SN yields, and modifying the the duration of the chemical evolution or their outflow feedback fraction, but the general trends in the IMF slopes remain.
Our IMF results do not necessarily imply that the IMF was simultaneously bottom heavy and top heavy, and that we derived a static shape of the IMF.We take our results as a confirmation of the galaxy wide IMF hypothesis, where the gIMF is not universal nor static, but it changes with time.The high-mass end IMF is representative of the situation in the first few hundred Myr of galaxy evolution and is characterised by an excess of massive stars, coinciding with the first intensive star burst.As the star formation changes to a more extended and less violent phase, the gIMF transits to the low-mass end IMF that we measure in nearby (massive) galaxies.It is characterised by an excess of low mass stars, having the chemical imprint of the high mass stars.Such a temporal variation of the gIMF was previously invoked to explain the high metallicities, enhanced -elements and bottom heavy IMFs in massive galaxies.Our results show that a combination of top and bottom heavy IMFs, at different epochs of formation are indeed required to explain the chemistry of low mass stars in very massive galaxies.
The small differences between the BCGs and satellite galaxies (in the elemental abundances and the high-mass IMF slope) can be understood as small variation in the formation paths of these two galaxy types.All investigated galaxies are very massive and live in dense environments, but the dependence of elemental abundance (and IMF) on the velocity dispersion, together with evidence based on stellar kinematics, suggest that BCGs experienced more gas-free similar mass mergers than the satellite galaxies.

APPENDIX A: AN EXAMPLE OF SKY SUBTRACTION METHOD
In Fig. A1 we present an example of our method to remove the remaining sky contribution from the galaxy spectra.As briefly explained in Section 3.1, we combine MUSE spectra outside the central effective radius aperture, and attempt to fit a stellar population model.This fit is typically not very well constrained as the outer spectra have essentially no information on the stellar continuum.The sky residuals are then obtained subtracting the stellar population fit from the outer spectrum.In principle, we could have also estimated the sky residuals by removing the median, or by a simple polynomial fit of the outer spectrum.For PGC003342 the difference in sky residuals obtained in these ways are far below 1%, except for two features at 5304Å and 7160Å, possibly corresponding to some Fe and Ca features, where they just approach the 1% level.They are carried over to the sky spectrum, but their influence to the final results are negligible.Presented sky residuals are typical for all other galaxies.The residual sky spectrum is subtracted from the inner spectrum and the result is shown on Fig. A1 with the thick teal line, in comparison with the original spectrum produced by the data reduction pipeline (thick black line).The improvement is obvious, even though some spectral regions remain noisy.We mask those during the alf fitting.Note also that below 4960Å we did not attempt to subtract the sky.

APPENDIX B: ELEMENT ABUNDANCES AND SYSTEMATIC UNCERTAINTIES
We first present here correlations between different parameters resulting from the nominal alf fits as presented in Section 4.1.In particular, we correlate the elemental abundances (and stellar ages) with the velocity dispersion in Fig. B1 and with the low-mass IMF slope in Fig. B2.On Fig. B1 we overplot the values from work by Gu et al. (2022), also obtained by alf, from Magellan/LDSS3 long-slit spectra of massive galaxies.These spectra have a comparable S/N (∼ 100 − 200) to our data, but their wavelength range is somewhat larger (4000-10300), and they do not cover the same one-effective radius aperture as we do.The alf set up also differs in several ways.Notably, they use an IMF model with three slopes, where the slope for the stellar masses higher than 1 M ⊙ is fixed to the Salpeter value.Fit is also divided into five spectral ranges, masking several sky lines.
There are some differences in the recovered elemental abundances, notably for [O/Fe] and partially for [Mg/Fe], while the metallicity and [Na/Fe] are in a good correspondence.Placing the observed differences between the two studies into context requires however also a broader comparison.As we show in Fig. 8, different alf setups deliver different results, with the overall largest difference when the wavelength range is limited to the blue region (< 690 nm).The most interesting are the four elemental abundances used to constrain the chemical models (see end of Section 3.2).These abundances show relatively small changes between various alf runs (∼ 0.05 dex), largest begin when fitting blue region only (∼ 0.1 dex).This is consistent with the conclusions of Gonçalves et al. (2020), who state that the results of the full spectral fitting (in their case for age, metallicity and [Mg/Fe]) can be dependant on the spectral range used.A similar results was already seen by Conroy et al. (2014), when they limited the fitting range to < 580 nm, except that in their cases [O/Fe] and [Na/Fe] showed larger discrepancies.
Given that one can expect differences in recovered parameters depending on various assumptions and the set up of the fit, we compare our results with a compilation of several literature studies which derive the metallicities and -elements abundances in diverse and independent ways.Specifically, we select works that use both the index fitting and the full spectral fitting methods.The results are shown on Fig. B3.There are a few noteworthy points to consider: • The data come from a variety of observations, using both longslits and integral-field units, and therefore covering different regions of galaxies.Walcher et al. (2015) data are based on SDSS 3 ′′ fibre spectra, ATLAS 3D (McDermid et al. 2015) and SAMI (Scott et al. 2017) data are based on one effective elliptical apertures, similar to our case.Other data sets are taken from long-slit observations, and spatially limited.
• The uncertainties are often comparable with our values, when the same velocity dispersion range ( > 250 km/s), as a proxy for galaxy mass, is considered, as already discussed in Section 5.1.
• The differences between our and literature results are similar to the differences between various literature results.
• A number of literature studies have comparable metallicity ([Fe/H]) or element abundances ([Mg/Fe]), but not both at the same time.Our metalicities are similar to those of Gu et al. (2022) and Walcher et al. (2015), while element abundances are similar to those of Thomas et al. (2005);Sánchez-Blázquez et al. (2006b) and Loubser et al. (2009).Our metallicity and abundances are similar with those of McDermid et al. (2015) and Scott et al. (2017), which are also relatively consistent with each other over a larger galaxy mass range.These similarities and differences are not related to the method (full spectral of index fitting): for example the SAMI, our and those of the ATLAS 3D abundances were obtained with two different types of full spectral fitting methods and an index fitting method, respectively, using different single stellar population models, but they still provide consistent results.
Overall, we can conclude that our results are consistent with those found in the literature, for the same galaxy mass (velocity dispersion).The comparison is essentially limited to systematic effects, which include the S/N of the data, aperture size, and perhaps most importantly, the assumptions on the stellar populations models.
In Table B1 we summarise the values of the Pearson coefficients for relations between element abundances and velocity dispersion, low mass IMF slope, as well as the low and high mass IMF slopes with the velocity dispersion.The correlation coefficients were calculated for three different samples: the full M3G sample, and for BCGs and satellites separately.The errors were estimated at a 2 level by taking the lower 2.5% and upper 97.5% limit values by bootstrapping.Here we list the most noteworthy relations, where at least one sample shows a significant correlation.An exception is the  LM −  relation which does not show a correlation, but is added for comparison with the  HM −  relation, which is significant for both the full sample and the satellites only.2017) -(S17) red line.T05 and SB06 show average values for metallicity and alpha abundance from these surveys, together with the average uncertainties (shaded regions).S17 shows the functional dependences fitted to their data.We used both full spectral and index fitting methods for comparison, where T05, SB06, L09 and M15 use index fitting methods.The label M15 SFH refers to the metallicities obtained by M15 using the full spectral fitting (the same data as for the index fitting).S17 and W15, and our work, use full spectral fitting.As in the rest of the paper, BCGs are shown with light red and satellites with light blue symbols.
newly inflowing gas has zero metallicity (Rybizki et al. 2017).Similarly, setting setting x  = 0.2 determines the fraction of enriched gas that is removed from the ISM, and therefore lost for direct star formation.While this gas will return to the ISM from corona, its chemical composition will be first diluted with the primordial gas.The results of these chemical models (not shown) are actually similar to previous tests: there are certain changes to some parameters (i.e. the mass fraction of corona), but they are no more pronounced than the changes induced by longer duration of the chemical models.Crucially, the high-mass end IMF slope  HM changes very little.For the case of no outflow feedback fraction there is a minor overall steepening (for an average of 0.02), while the correlations between  HM and   remain the identical.The new Pearson correlations coefficients are: the full sample is  = −0.57+0.31 −0.18 , for BCGs  = −0.28+0.69  −0.36 and for SAT  = −0.86+0.11 −0.07 .In the case of the fixed outflow feedback fraction, the slope actually decreases (for an average of 0.05), but the correlations remain the same.The new Pearson correlations coefficients are: the full sample is  = −0.49+0.41  −0.26 , for BCGs  = −0.12+0.66 −0.41 and for SAT  = −0.85+0.11 −0.08 .We conclude that the main chemical modelling results pertaining to the high-mass end IMF slopes are robust.

APPENDIX D: ALF SPECTRAL FITS
This paper has been typeset from a T E X/L A T E X file prepared by the author. .Chempy results as a function of velocity dispersion.We show the results of two chemical models, the nominal models with the chemical evolution lasting for 3.5 Gyr (circles) and a test model with the chemical evolution lasting 7 Gyr (squares).The youngest galaxy, PGC073000 is highlighted with a large circle and a square.From left to right, top to bottom, we show versus the stellar velocity dispersion within R eff : a) the star formation time scale , b) the best-fit SN Ia delay time, c) the amount of SN Ia exploding over a period of 15 Gyr per solar mass, d) fraction of the outflowing enriched gas deposited in the corona instead of the star forming zone (note that 0.001 is a lower limit imposed through a prior), e) the mass fraction of the corona, f) the star formation efficiency, g) the power law index of the high-mass IMF, h) the shape parameter  of the Gamma distribution, and i) the fraction of high-mass supernovae exploding as core-collapse supernova instead of as hypernova.BCGs and satellite galaxies are shown with red and blue symbols, respectively, as indicated on the legend of panel h).
Figure1.Mass -size relation for M3G galaxies (circles), divided into BCGs (orange) and satellites (light blue).Squares are galaxies from the MASSIVE survey (values taken fromMa et al. 2014), and pentagons show massive part of the ATLAS 3D Survey of nearby early-type galaxies (values taken fromCappellari et al. 2013).Sizes of galaxies from all surveys are based on the 2MASS observations.The masses of the MASSIVE galaxies are based on the 2MASS K-band absolute magnitude and calibration eq.(2) fromCappellari (2013).ATLAS 3D masses come from dynamical models.M3G masses are based on the virial estimate and taken fromKrajnović et al. (2018).The dashed lines are the lines of the constant velocity dispersion, based on theCappellari et al. (2006) parametrisation of the virial mass estimator.

FluxFigure 2 .
Figure 2. Example of the best-fit alf model to the MUSE data of PGC003342.Grey patches show spectral regions that were masked because of telluric absorption, sky emission or line emission.The data-model residual is shown in green (grey) in unmasked (masked) regions, and shifted to 0.5 the average flux value.The blue region denotes the 1- uncertainty on the data.Similar spectra for other galaxies are presented in Fig. D1.

Figure 4 .
Figure 4. Selected elemental abundances as a function of global velocity dispersion (left) and low-mass end IMF slope (right).BCGs are denoted with orange points, satellite galaxies with light blue points.We plot here only elements for which a significant correlations, either for BCGs or satellite galaxies.The Pearson correlation coefficient  (for the full sample) is indicated on each panel, with uncertainties indicating the 2 level.The correlations for all other elements are in Figs.B1 and B2.For comparison we include data from the MASSIVE sample(Gu et al. 2022) with small green symbols.

Figure 5 .
Figure5.Low-mass end slope of the IMF,  LM , versus velocity dispersion   within the effective radius aperture.BCGs are denoted with orange points, satellite galaxies with light blue points.
(1): name of the galaxy; Column (2): global velocity dispersion estimated within an elliptical aperture of a semi-majr axis radius equal to the half-light radius of the galaxy.Values in the subsequent columns are also estimated within the same aperture.Column (3): low-mass end IMF slope from the nominal alf run; Column (4): high-mass end IMF slope from the nominal chempy run; Column (5): low-mass end IMF slope when Na 8190 line is masked;Column (6): low-mass end IMF slope when NA 8190 and NaD lines are masked; Column (7): high-mass end IMF slope usingChieffi & Limongi (2004) yields.Note that for  CL04 we did not run MCMC chains to calculate the uncertainties for individual galaxies, but the errors are expected to be similar to  HM .Column (8): BGCs are marked with 1.

Figure 7 .
Figure7.High-mass end slope of the IMF derived using Chempy with the yields ofNomoto et al. (2013) versus the low-mass end slope derived using alf.BCGs and satellite galaxies are shown with orange and light blue symbols, respectively.Galaxies with prolate-like rotation (with kinematic misalignments larger than 75 • ) are shown as green pentagons.The colour of the error bars of the pentagons relate them to BCGs or satellites.Horizontal dashed lines indicate the Salpeter (higher) and Kroupa (lower) values for the low-mass end slopes, for comparison purposes.

Figure 9 .
Figure9.The relation between the low-mass end IMF slopes derived using alf by masking Na lines (-axis) and not masking (-axis).We used two different masks, only applying it to the Na 8190Å line, or both on 8190Å line and NaD doublet.Only small differences are seen between two masks and there is an overall good agreement with the unmasked values.

Figure 10 .
Figure 10.Elemental abundances from observations as a function of central velocity dispersion for a subset of light elements.All elements are shown with respect to the Fe abundance.Lines denote model predictions based on the average parameters from the Chempy fits to the observational data.Note that fits included only the elements O, Na, Mg and Fe.Models are based on the yields of Nomoto et al. (2013), Chieffi & Limongi (2004) and Limongi & Chieffi (2018) (set R).

Figure 11 .
Figure11.High mass slope of the IMF derived using the yields ofChieffi & Limongi (2004) versus the low-mass end slope derived using alf.BCGs and satellite galaxies are shown with orange and light blue symbols, respectively.Note the similarity of this figure with Fig.7, once the systematic shift to lower  HM is take into account.Horizontal dashed lines indicate the Salpeter (higher) andKroupa (lower)  values for the low-mass end slopes, for comparison purposes.
Figure12.A composite distribution of the initial mass functions in the  () =    ∼  −  form for M3G galaxies.The IMFs are constructed using the low-mass (< 0.5 M ⊙ ) and high-mass (> 0.5 M ⊙ ) end slopes obtained from the alf analysis of the MUSE spectra and the chempy modelling of recovered elemental abundances.For comparison, we plot also the Kroupa (dashed line) and the Salpeter (solid line) IMFs.Orange lines are BCGs and light blue lines are the satellite galaxies of the M3G sample.Vertical dashed-dotted line approximately shows the boundary between stars that still exits in M3G galaxies, and stars that ended their evolutionary paths.Note that this plot is not meant to represent a fixed IMF.Instead it shows a temporal variation of the IMF, in the sense of the IGIMF theory (e.g.Jeřábková et al. 2018), where, for the massive galaxies, the high-mass end is representative of the very early age, and the low-mass end of the later star formation.

Figure B1 .
Figure B1.Elemental abundances based on the analysis of MUSE spectra of M3G galaxies with alf as a function of velocity dispersion within the effective radius.BCGs are shown with orange and satellite galaxies with light blue symbols.The top left plot shows the luminosity weighted (log 10 ) age within the same aperture as the elemental abundances.The Pearson correlation coefficient  is indicated on each panel.Uncertainties on the velocity dispersion are smaller than the symbols (∼ 2 km/s).For comparison we include the data fromGu et al. (2022) with small green symbols.Note that [Fe/H], [N/Fe], [Na/Fe], [Mg/Fe], [V/Fe] and [Ni/Fe] are repeated here from Fig. 4 for completeness.

Figure B2 .Figure B3 .
Figure B2.Elemental abundances as a function of the low mass IMF slope for M3G galaxies derived with alf.BCGs are shown with orange and satellite galaxies with light blue symbols.The Pearson correlation coefficient  is indicated on each panel.The strongest correlations are found between the IMF slope and [Mg/Fe] and [Sr/Fe].While there are clear trends for BCGs for some elements (e.g.[Ti/Fe], [Ba/Fe]), there are no statistically secure differences between BCGs and satellite galaxies.Note that [N/Fe], [Mg/Fe], [Ti/Fe], [Co/Fe].[Sr/Fe] and [Ba/Fe] are repeated here from Fig. 4 for completeness.
Figure C1.Chempy results as a function of velocity dispersion.We show the results of two chemical models, the nominal models with the chemical evolution lasting for 3.5 Gyr (circles) and a test model with the chemical evolution lasting 7 Gyr (squares).The youngest galaxy, PGC073000 is highlighted with a large circle and a square.From left to right, top to bottom, we show versus the stellar velocity dispersion within R eff : a) the star formation time scale , b) the best-fit SN Ia delay time, c) the amount of SN Ia exploding over a period of 15 Gyr per solar mass, d) fraction of the outflowing enriched gas deposited in the corona instead of the star forming zone (note that 0.001 is a lower limit imposed through a prior), e) the mass fraction of the corona, f) the star formation efficiency, g) the power law index of the high-mass IMF, h) the shape parameter  of the Gamma distribution, and i) the fraction of high-mass supernovae exploding as core-collapse supernova instead of as hypernova.BCGs and satellite galaxies are shown with red and blue symbols, respectively, as indicated on the legend of panel h).

5 10 15 20 25 30 Atomic number 0.50 0.25 0.00 0.25 0.50 [X/Fe] C N O NaMg Si K Ca Ti V Cr MnFe Co Ni Cu BCGs SATs Figure 3. Abundances
of elements in M3G sample galaxies, with respect to [Fe/H] within the effective radius for BCGs (orange) and satellite galaxies (light blue), normalised to the solar abundances.Shaded regions show the limits of the distributions estimated as the 16th and 84th percentiles for each element.totalmass in star formation.Star formation in the central zone happens with a star formation efficiency SFE, which defines the mass of interstellar medium (ISM) gas required to form SFE×mass of stars.Supernovae Ia are assumed to have no mass dependence.The delay time distribution has an initial delay of  Ia Gyr, for which we assume a Gaussian prior width   Ia = 0.5 Gyr and   Ia = 0.3 Gyr.The delay time distribution has a shape  −1.12 based on Maoz et al.
Notes: Columns with "lower" and "upper" headers show the 16th and 84th percentile values of the element abundance distributions.

Table 2 .
Low-and high-mass end IMF slopes for M3G galaxies.