High-z Dusty Star-forming Galaxies: A Top-heavy Initial Mass Function?

Recent estimates point to abundances of z > 4 submillimeter galaxies far above model predictions. The matter is still debated. According to some analyses the excess may be substantially lower than initially thought and perhaps accounted for by flux boosting and source blending. However, there is no general agreement on this conclusion. An excess of z > 6 dusty galaxies has also been reported albeit with poor statistics. On the other hand, evidence of a top-heavy initial mass function (IMF) in high-z starburst galaxies has been reported in the past decades. This would translate into a higher submillimeter luminosity of dusty galaxies at fixed star formation rate, i.e., into a higher abundance of bright high-z submillimeter galaxies than expected for a universal Chabrier IMF. Exploiting our physical model for high-z protospheroidal galaxies, we find that part of the excess can be understood in terms of an IMF somewhat top-heavier than Chabrier. Such an IMF is consistent with that recently proposed to account for the low 13C/18O abundance ratio in four dusty starburst galaxies at z = 2–3. However, extreme top-heavy IMFs are inconsistent with the submillimeter counts at z > 4.


Introduction
The opening of the submillimeter window on galaxy formation and evolution, thanks to the 850 μm surveys with the Submillimeter Common-User Bolometer Array (SCUBA) on the James Clerk Maxwell Telescope (e.g., Smail et al. 1997; Barger et al. 1998;Hughes et al. 1998), has revolutionized the field. SCUBA blank-field pointings demonstrated that the most active star formation phases of high-z galaxies are heavily dustenshrouded and therefore largely missed even by the deepest optical surveys (for a review, see Casey et al. 2014). The results of these surveys have strongly shaken up the leading galaxy formation paradigm of the time which underpredicted the submillimeter counts by one order of magnitude or more (e.g., Kaviani et al. 2003;Baugh et al. 2005).
Part of the discrepancy was later shown to be attributable to an overestimate of the observed counts (Coppin et al. 2006;Geach et al. 2017), mainly due to source blending within the relatively large SCUBA beam and to insufficient corrections for flux boosting (Eddington bias;Eddington 1913). Nevertheless, submillimeter data are still challenging for semianalytic galaxy formation models (SAMs). In particular, SAMs generally underpredict the bright end of the star formation rate (SFR) function at  z 2 (Niemi et al. 2012;Gruppioni et al. 2015). A further challenge came from the first searches for submillimeter selected z4 galaxies (Dowell et al. 2014;Asboth et al. 2016;Ivison et al. 2016): the inferred surface densities at flux densities larger than several tens of mJy at 500 μm were found to be in excess of model predictions by about one order of magnitude. More recent studies (Donevski et al. 2018;Duivenvoorden et al. 2018) suggest that the discrepancy with models may be accounted for by flux boosting due to instrumental noise and confusion (including the contribution from clustering), plus source blending. However, the issue is not settled yet (see Section 2).
On the other hand, a key ingredient to understand the star formation process is the initial mass function (IMF), i.e., the relative number of newly formed stars as a function of their mass. It is still being debated whether or not the IMF is universal or varies with conditions in star-forming regions. Bastian et al. (2010), in their comprehensive review, concluded that the data available at the time did not provide "clear evidence that the IMF varies strongly and systematically as a function of initial conditions after the first few generations of stars." On the theoretical side, Padoan et al. (1997) presented a model whereby the IMF depends on the average physical parameters (e.g., temperature, density, and velocity dispersion) of the star formation sites. 5 The model predicts a larger fraction of massive stars (i.e., a top-heavier IMF) in starbursts. In starburst galaxies the cosmic-ray energy densities generated by the deaths of the most massive stars may raise the gas temperatures even in the coolest star-forming clouds, thus increasing the Jeans mass, hence the characteristic stellar mass, i.e., the IMF shape (Papadopoulos et al. 2011). Chiosi et al. (1998 argued that this model alleviates or solves some difficulties faced by models for the formation and evolution of elliptical galaxies. Recently, observational evidence of a top-heavy IMF (i.e., of an IMF with more massive stars than expected from canonical IMFs, such as that by Chabrier 2003 or by Kroupa et al. 2013) have been reported to account for the low 13 C/ 18 O abundance ratio found in both local ultraluminous infrared galaxies (Sliwa et al. 2017;Brown & Wilson 2019) and in four gravitational lensed submillimeter galaxies at z∼2-3 (Zhang et al. 2018).
A greater abundance of massive stars at high-z translates into a higher submillimeter luminosity of dusty galaxies at fixed SFR, i.e., into a greater abundance of submillimeter bright galaxies than expected under the assumption of a universal IMF. How do the counts of high-z galaxies implied by the topheavier IMFs proposed by Zhang et al. (2018) compare with observations?
To address this question we need a physical model, linking the IR luminosity of galaxies to their star formation history, taking also into account the effect of active galactic nuclei (AGNs) growing in their centers. The physical model by Cai et al. (2013), briefly described in Section 3, allows us to carry out this investigation. The baseline model adopts a Chabrier (2003) IMF, assumed to be universal. This IMF is described by a broken power law: µ where m is the stellar mass and dN is the number of stars within d log m.
An extremely top-heavy IMF ( = dN d m log constant) was previously proposed by Baugh et al. (2005) to bring the counts of submillimeter galaxies predicted by their semianalytic model into agreement with observations. Less top-heavy IMFs have been invoked by Lacey et al. (2016) and Cowley et al. (2019) to predict submillimeter counts. However other models (e.g., Granato et al. 2004;Cai et al. 2013) accounted for the submillimeter counts adopting a standard (Chabrier 2003) IMF. None of these models considered the counts of submillimeter selected galaxies at z4 or z6.
In Section 4 we investigate the effect of a top-heavier than Chabrier IMF on several observables. In Section 5 we summarize and discuss our main conclusions. Throughout this paper we adopt a flat ΛCDM cosmology with the latest values of the parameters derived from Planck Cosmic Microwave Background (CMB) power spectra: 2. Dusty Star-forming Galaxies at z>4

Summary of Relevant Observations
The submillimeter extragalactic surveys with the Herschel Space Observatory at wavelengths extending over the peak of the dust emission spectrum up to high redshifts, have made it possible to select candidate high-z galaxies. The first searches for z4 galaxies were carried out by Dowell et al. (2014), Asboth et al. (2016), and Ivison et al. (2016). They used Herschel Spectral and Photometric Imaging Receiver (SPIRE) flux densities at 250, 350, and 500 μm to search for sources with red colors, indicative of high z. Since the spectral energy distribution (SED) of local ultraluminous dusty star-forming galaxies typically peaks at rest-frame wavelengths λ∼100 μm, the SEDs of z4 dusty star-forming galaxies are expected to typically peak at l m  500 m in the observer's frame (500 μm risers).
On this basis, Dowell et al. (2014) Oliver et al. 2012), totaling an area of 21 deg 2 . Follow-up spectroscopy of the first 5 of the 38 sources selected in this way showed that 4 have z in the range 4.3-6.3 and the fifth has z=3.4, confirming their method is efficient at selecting high-z dusty galaxies. After correcting for completeness and purity, the surface density of risers with m  S 30 500 m mJy was found to be  . Submillimeter images with SCUBA-2 and with LABOCA were obtained for a suitably chosen subset of these sources. Photometric redshifts were estimated comparing the observed submillimeter SEDs with those of three well sampled SED templates, found to be representatives of high-z dusty starforming galaxies. In this way, the redshift distribution was derived. The fraction of z>4 galaxies was found to be 32±5%. Their surface density is  -2.1 0.6 deg 2 , somewhat lower than, but consistent with the estimate by Dowell et al. (2014) and substantially higher than the Cai et al. (2013) and the Béthermin et al. (2017) predictions.
However, as suggested by more recent studies (Donevski et al. 2018;Duivenvoorden et al. 2018), the discrepancy with models may be accounted for by flux boosting and source blending. Donevski et al. (2018) Asboth et al. (2016) and estimated their redshifts by SED fitting. Both groups used model-based simulations to correct the raw counts for the effects of flux boosting and of source blending. On this footing they found consistency with model predictions (see Figure 1). Jin et al. (2018) have published a "super-deblended" farinfrared to (sub)millimeter photometric catalog in the Cosmic Evolution Survey (COSMOS) field with spectroscopic or photometric redshifts. We have extracted from their catalog z>4 sources with > m S 20 500 m mJy and signal-to-noise ratio 5 and derived their differential number counts after applying the correction for flux boosting (see Equation (1)). Again, the corrected counts are consistent with model predictions.

Potential Observational Biases
A strong effect of flux boosting on the counts of z4 galaxies is indeed expected as a consequence of their steep slope. As shown by Hogg & Turner (1998), in the case of Gaussian noise, the ratio of the maximum likelihood true flux, S ML , to the observed flux, S o , is related to the slope of the integral counts, q, and to the signal-to-noise ratio, r, by The observationally determined counts shown in Figure 1 consistently indicate, for > z 4 galaxies,q 3.6, implying that there is no finite maximum likelihood value for r<4.3. For r=4.5, close to the critical value r=4.3, and r=5, usually taken as the threshold for reliable detections, we have and 0.757, respectively. This bias in the measured flux density translates to a bias in the integral source counts of , of factors of 4.7 and 2.7, respectively, for the two values of r. Note that these are just examples to illustrate that noise has indeed a strong impact on observed counts. This simple estimate does not account for the non-Gaussian, spatially correlated fluctuations due to clustering as well as for the complexity of source selection. Both Donevski et al. (2018) and Duivenvoorden et al. (2018) have extensively used simulations from Béthermin et al. (2017) to evaluate their selection method. To deal with these effects, Béthermin et al. (2017) have produced a 2 deg 2 simulation of the extragalactic sky from 24 μm to 1.3 mm, called Simulated Infrared Dusty Extragalactic Sky (SIDES). The simulation uses an updated version of the empirical model by Béthermin et al. (2012Béthermin et al. ( , 2015 and includes clustering based on the spatial distribution of dark matter halos yielded by a state-of-the-art Nbody simulation; halos are populated with galaxies using abundance matching. Since a 2 deg 2 cone is too small to provide enough statistics on the rare z>4 sources selected by most of the studies mentioned above, Béthermin et al. (2017) built also a 274 deg 2 simulation, which however does not contain clustering. The simulation of the full source extraction showed that the combination of noise, source blending, and weak lensing can dramatically boost the number density of red sources selected with the criteria by Asboth et al. (2016), compared to their real number density. This conclusion contradicts that by Asboth et al. (2016) who also used simulations to support their surface density estimates. Béthermin et al. (2017) explained the contradiction claiming that they have identified a potentially incorrect assumption in the simulation by Asboth et al. (2016). Béthermin et al. (2017) attributed a substantial contribution to the number density overestimate to an overestimate of source flux densities at 500 μm due to multiple sources within the relatively large Herschel beam (source blending) at this wavelength; according to their simulation only ∼60% of the measured flux density comes from the brightest source within the beam. Accounting for these effects can reconcile observations with model predictions. Béthermin et al. (2017), however, acknowledged that correcting the observed counts for all the observational effects is a complex task and that their approach has limitations.
The issue is not settled, however. None of the galaxies, drawn from the sample of z4 candidates selected by Ivison et al. (2016), observed by Fudamoto et al. (2017) in 1.3-or 3mm continuum using the Northern Extended Millimeter Array and Atacama Large Millimeter Array (ALMA) were found to be affected by blending, as would be expected from simulations of Béthermin et al. (2017); multiplicity may affect at most the ∼20% of their targets that remained undetected in continuum.
Further hints of a high-z excess come from searches of z>6 dusty galaxies. Riechers et al. (2013) Figure 2).

Model
While at z1.5 star formation mostly occurs in late-type galaxies (spiral, irregular, and disky starburst galaxies), at the high redshifts of interest here the observed IR/submillimeter galaxy luminosity functions are dominated by protospheroids, i.e., by ellipticals and bulges of disk galaxies in the process of forming the bulk of their stars (e.g., Granato et al. 2004;Thomas et al. 2010;Silk & Mamon 2012). The Cai et al. (2013) model provides a physically grounded description of their formation and evolution, building on the work by Granato et al. (2004) and Lapi et al. (2006Lapi et al. ( , 2011. The function. For the latter, the analytical approximation by Sheth & Tormen (1999) was used. High-resolution N-body simulations (e.g., Wang et al. 2011) have identified two phases of the halo evolution. A first fast collapse phase, including major mergers, forms the bulk of its mass. The following phase produces a slow growth of the halo outskirts by many minor mergers and diffuse accretion. This second stage has little effect on the inner potential well where the visible galaxy resides. Based on these results, the model assumes that the star formation and the AGN growth are basically regulated by in situ processes, as confirmed by subsequent high-resolution simulations (e.g., Pillepich et al. 2015). The star formation history of protospheroids with any halo mass and formation redshift and the accretion history onto the central supermassive black holes are computed by solving a set of equations describing the evolution of gas phases and of the active nucleus, including the effect of cooling, condensation into stars, radiation drag, accretion, and feedback from supernovae (SNe) and from the AGN.
These equations yield the SFR of each galaxy and the bolometric luminosity of the AGN as a function of halo mass, formation redshift, and galactic age. The SFR was converted to the total infrared (IR; 8-1000 μm) luminosity using the standard calibration (e.g., Kennicutt & Evans 2012). The epoch-dependent bolometric luminosity functions of galaxies, of AGNs, and of objects as a whole (galaxy plus AGN) are computed coupling luminosities with the halo formation rate.
To derive monochromatic luminosity functions, the SED of the strongly lensed z;2.3 galaxy SMMJ2135-0102 (Swinbank et al. 2010), modeled using the GRASIL software package (Silva et al. 1998) was adopted for the dusty galaxy component. This SED was found to be representative of high-z dusty galaxies (González-Nuevo et al. 2012;Ivison et al. 2016). For AGNs two main phases of the coevolution with the galaxy component were considered. For the first phase, when the black hole growth is enshrouded by the dusty interstellar medium (ISM) of the host galaxy, a heavily absorbed AGN SED, taken from the library by Granato & Danese (1994), was adopted. For the second phase, when the AGN shines after having swept out the galaxy ISM, we adopted the mean QSO SED by Richards et al. (2006), extended to submillimeter wavelengths assuming a graybody emission with dust temperature T dust =80 K and emissivity index β=1.8.
After determining the values of the parameters from a comparison with a set of multifrequency luminosity functions (see Cai et al. 2013 for details) the model successfully reproduced a broad variety of data: multifrequency and multiepoch luminosity functions of galaxies and AGNs beyond those used for the model calibration, counts of unlensed and gravitationally lensed sources, redshift distributions, and clustering properties (Xia et al. 2012;Cai et al. 2013Cai et al. , 2014Carniani et al. 2015;Bonato et al. 2017). In particular the model accounts for the redshift-dependent IR luminosity function determined by Gruppioni et al. (2013) up to z;4 (Bonato et al. 2014) and fits the most recent millimeter and submillimeter counts and redshift distributions De Zotti et al. 2019;Gralla et al. 2019) as well as the SED and the power spectrum of the cosmic infrared background (Cai et al. 2013).

A Top-heavier IMF?
The bolometric luminosity of a young stellar population is dominated by massive, short-lived, UV-producing stars (e.g., Kennicutt & Robert 1998;Madau et al. 1998). Hence it is approximately proportional to the mass fraction in massive (m8 M e ) stars. For the heavily dust obscured galaxies of interest here, L L IR bol  . The calibration of the relationship between L IR and SFR for complete dust obscuration adopted by Cai et al. (2013)  The mass fractions in the 8-100 M e range are of 23%, 33%, and 44% for the "Chabrier," "top-heavy," and "Ballero" IMF, respectively. Thus the bolometric luminosity at fixed SFR increases, compared to the "Chabrier" case, by factors of ;1.4 and of ;1.9 for the "top-heavy" and "Ballero" IMF, respectively. In the extreme case of the flat IMF (i.e., = dN d m log constant) used by Baugh et al. (2005) for starbursts, the mass fraction in  M 8  stars increases by a factor of ;4. Zhang et al. (2018) mention a factor of 6-7 increase in the mass fraction of massive stars moving from the "Kroupa" (Kroupa et al. 2013) to the "Ballero" IMF. This may look in conflict with the substantially lower factor we find for the "Chabrier" case since the "Kroupa" IMF is generally considered to give a relation between L IR and SFR almost identical to the "Chabrier" IMF (e.g., Chomiuk & Povich 2011). However, Zhang et al. (2018) refer to a version of the Kroupa IMF with a much lower fraction of massive stars (8.6%) 6 than that (23%) yielded by the IMF adopted by Cai et al. (2013).
The IMF impinges also on other quantities: the SN rate and the fraction, , of the mass of a stellar generation that returns to the ISM (restitution fraction). Both factors affect the evolution of the metallicity and of the SFR. An increase of  makes more material available for star formation while feedback effects are obviously stronger for higher SN rates, with the effect of decreasing the SFR.
SNe dominate feedback effects not only for less massive galaxies but also for most of the intense star formation phase of massive galaxies. The AGN feedback takes over in the latter objects only during the last few e-folding times of the central black hole growth. The fraction of stars per stellar generation that explode as SNe, i.e., the ratio of the number of stars more massive than 8 M e to the total number of stars, is ;1% for the "Chabrier" IMF and increases to ;1.4%, 2.2%, and 37% for the "top-heavy," "Ballero," and "Baugh" IMFs, respectively. The corresponding parameter defined in Cai et al. (2013) is the number of SNe per unit solar mass of condensed stars: N SN ; 1.2×10 −2 /M ☉ , 1.6×10 −2 /M ☉ , 2.0×10 −2 /M ☉ , and 2.5×10 −2 /M ☉ , for the "Chabrier," "top-heavy," "Ballero," and "Baugh" IMFs, respectively.
The restitution fraction, , affects the SFR and the chemical evolution. It depends on the minimum mass, M min , of stars with lifetimes not exceeding the age of the universe at the considered redshift. We set M min =2 M e since the lifetime of a star of that mass with solar metallicity is slightly lower than the age of the universe at = z 4 (;1.5 Gyr) for our choice of cosmological parameters. For = M M 2 min  , =  0.393, 0.466, 0.56, and 0.83 for the "Chabrier," "topheavy," "Ballero," and "Baugh" IMFs, respectively. After having checked that the results are very weakly sensitive to the choice for M min , we have neglected its dependence on metallicity and used the same value for protospheroids at all redshifts. Note, however, that at  z 1.5 the dominant starforming population are late-type galaxies, for which we have kept the same model used by Cai et al. (2013).
For each of the considered IMFs we have rerun the code solving the set of equations written down by Cai et al. (2013).
To compute observable quantities such as number counts or luminosity functions, the output of these calculations, i.e., the SFR and the bolometric luminosity of the central AGN as a function of halo mass, formation redshift, and galactic age, need to be coupled with the SEDs of source populations.
We have kept the AGN SEDs used by Cai et al. (2013). In the case of dust luminosities powered by star formation, we need to take into account that higher luminosities entail, on average, higher dust temperatures, i.e., warmer SEDs (Schreiber et al. 2018;Liang et al. 2019). Since a sophisticated treatment is not warranted in this exploratory study we have taken this into account by assuming that the effective dust temperature, T d , hence the peak frequency of the dust emission, ν peak , scales as L IR 1 6 (e.g., Liang et al. 2019). Then, for different IMFs, n µ  k ; peak ,IR 1 6 the full continuum SED was shifted in frequency by the corresponding factor. Figure 3 shows the SEDs associated to the considered IMFs. This approach assumes, for each choice of the IMF, the same effective SED for all galaxies. This is clearly an approximation, necessary to keep the model manageable and justified by its capability of accounting for a broad variety of data. Figure 3. Galaxy SEDs, with total IR luminosity normalized to unit, associated to the "Chabrier," "top-heavy," "Ballero," and "Baugh" IMFs of protospheroidal galaxies. 6 The fraction of massive stars (8-100 M ☉ ) reported in Table 1 of Zhang et al. (2018) for the "Kroupa" IMF (6.9%) is a typo which does not affect their chemical evolution models (Z.-Y. Zhang 2020, private communication). With the corrected fraction of massive stars, 8.6%, the increase from their "Kroupa" to the "Ballero" IMF is a factor of ;5.  and =  0.54. These values differ from the values given above for the "Chabrier" IMF because of a slightly different way of computing N SN and of the choice of a lower M min (1 M e instead of 2 M e ) motivated by the fact that the focus was on lower redshifts. We have rerun the model with the new values and used the results of the updated model. The differences in the derived IR luminosity functions for the "Chabrier" IMF turned out to be irrelevant. Figure 1 shows that the observed raw 500 μm counts of z4 galaxies are consistently and substantially above the predictions of the baseline Cai et al. (2013) model, which adopted a universal Chabrier (2003) IMF. Taken at face value they strongly favor a top-heavier IMF such as the "top-heavy" or the "Ballero"; the extreme "Baugh" IMF is, however, inconsistent with the data. A similar indication comes from the estimate at z>6 (Figure 2). The correction factors derived from the Béthermin et al. (2017) simulations bring the z4 data close to agreement with the baseline model. However, several of the corrected data points are still consistent, within the uncertainties, with a top-heavy IMF; the main exception to that are the counts by Donevski et al. (2018). High-resolution follow-up of the two z>6 galaxies did not show indications of blending, so that no substantial correction appears to be necessary. However, the statistical uncertainty of this data point is large. Figure 4 compares the observed total IR (8-1000 μm) luminosity functions at z=1.8 to z=5.5 with results for the baseline Cai et al. (2013) Gruppioni et al. (2013) and by Magnelli et al. (2013) were drawn from deep surveys with the Herschel Photodetector Array Camera and Spectrometer which has a better angular resolution than Herschel/SPIRE. These data are therefore much less affected by confusion and blending problems. Wang et al. (2019) used Herschel/SPIRE data in the COSMOS field but exploited the rich, deep multifrequency data available in this field and the SIDES empirical model to deblend and deconfuse the Herschel/SPIRE images. In this way they were able to extend the submillimeter source counts by more than a factor of 10 below the confusion limit and to extend estimates of the IR luminosity function to fainter luminosities and out to higher redshifts (up to z∼6) than Gruppioni et al. (2013). The IR luminosity functions computed by Hatsukade et al. (2018) and by Magnelli et al. (2019) are based on high-resolution ALMA surveys and are therefore unaffected by confusion and blending.
The baseline model by Cai et al. (2013) yields results in remarkable agreement with the luminosity function data, especially at z;1.8 and ;2.2. An excess shows up especially at the highest redshifts (z 3.5  , 4.5, and 5.5). The significance of the discrepancy is hard to assess, however, because the published errors are purely Poisson while much larger uncertainties may be associated to the source selection, to photometric redshifts, and to the deblended photometry. Errors preferentially shift sources from more populated to less populated redshift/luminosity bins (an effect analogous to the Eddington bias), implying that the number of rare highest luminosity/redshift sources may be easily overestimated. However, taken at face value, these data support the view of an excess of high-z submillimeter galaxies that may require a top-heavier IMF, such as the "top-heavy" or the "Ballero," to be accounted for. Again, the "Baugh" IMF implies source densities well in excess of observations. ) used by Baugh et al. (2005).