{'fraction_non_alphanumeric': 0.030090669363869522, 'fraction_numerical': 0.01603872412672037, 'mean_word_length': 4.788499738630423, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 6, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Implications of a Temperature Dependent IMF I: Photometric Template Fitting

A universal stellar initial mass function (IMF) should not be expected from theoretical models of star formation, but little conclusive observational evidence for a variable IMF has been uncovered. In this paper, a parameterization of the IMF is introduced into photometric template fitting of the COSMOS2015 catalog. The resulting best-fit templates suggest systematic variations in the IMF, with most galaxies exhibiting top-heavier stellar populations than in the Milky Way. At fixed redshift, only a small range of IMFs are found, with the typical IMF becoming progressively top-heavier with increasing redshift. Additionally, subpopulations of ULIRGs, quiescent- and star-forming galaxies are compared with predictions of stellar population feedback and show clear qualitative similarities to the evolution of dust temperatures.


INTRODUCTION
The total luminosity of a galaxy is dominated by the most massive stars in its stellar population, yet these individual stars cannot be directly characterized. As a result, nearly everything that has been inferred about the properties of high-redshift galaxies relies on an assumed stellar initial mass function (IMF) (Bonnell et al. 2007;Conroy et al. 2009b). Several key quantities depend particularly strongly on the shape of the assumed IMF. For example, because the stellar mass is dominated by low-mass stars, inferred stellar masses and star formation rates (SFRs) are highly sensitive to the ratio of high-mass to low-mass stars in the IMF. In photometric template fitting, the most common technique used in large, high-redshift surveys, the inferred extinction, metallicity, and other properties are also covariant with the assumed IMF (cf. Brammer et al. (2008)). Thus, it is essential not merely to assume an IMF, but to assume the correct one.
On a fundamental level, star formation involves a competition between gravity and processes which provide support against collapse, most notably thermal pressure. Therefore, the IMF must be sensitive to anything which affects that balance. This should include a wide range of properties that differ between galaxies, including the central potential, existing stellar mass, star formation history, supernova rate, cosmic ray density and galactic magnetic fields, metallicity, dust density and composi-tion, AGN activity, and the environment and merger history. All of these are known to vary both between different galaxies at fixed redshift and across different redshifts (cf. Laigle et al. (2016a)). Therefore, it should be expected that the IMF is not universal, but rather differs between galaxies, and between different epochs for the same galaxy. In particular, the IMF should depend upon the gas temperature of star-forming clouds, with higher-temperature regions producing higher average stellar masses (Low & Lynden-Bell 1976;Jermyn et al. 2018). Since observations of dust even at moderate redshifts find an increase in temperature towards higher redshift (Papadopoulos 2010a;Casey 2012;Magnelli et al. 2014;Magdis et al. 2017), this implies that galactic IMFs should systematically vary with redshift, and thus the assumed IMF must do the same.
However, at present nearly all studies have assumed a universal IMF, choosing one of several estimated Galactic IMFs and applying that to every galaxy at all redshifts. In part, this is out of necessity. It has not yet been possible to directly observe the IMF outside of the Galaxy and local satellites. Indeed, even these local measurements have produced several proposed Galactic IMFs (Salpeter 1955;Kroupa 2001;Chabrier 2003). Furthermore, inferred properties such as stellar mass and SFR are known to be highly sensitive even to which of the proposed Galactic IMFs is assumed (cf. Speagle et al. (2014)). Attempts to constrain extragalactic IMFs, or even to determine whether the same IMF fits all galaxies, have produced little consensus (Cenarro et al. 2003;Treu et al. 2010;Cappellari et al. 2012). Although results at times have suggested both top-heavier (Conroy & van Dokkum 2012) and bottomheavier (Martín-Navarro et al. 2015) IMFs at high redshift, there is not yet conclusive evidence that any specific galaxy has a non-Galactic IMF.
This limited success has primarily been caused by two distinct problems. First, there are significant degeneracies between the IMF and extinction, metallicity, star-formation history, and the age of the stellar population. Thus, it is likely impossible to determine the entire shape of the IMF from high-redshift observations alone. Second, theory has not yet produced consensus predictions for the functional forms that IMFs might take at higher redshift.
In recent work, Jermyn et al. (2018) proposed a temperature-dependent family of IMFs to account for temperature variation. This family of IMFs, combined with cosmic ray feedback, was used to predict a temperature history for the star-forming "main sequence" (Steinhardt et al. 2020). Because this family of IMFs only introduces one additional parameter, T IM F , it is possible to fit photometric catalogs with the entire family of IMFs, producing a best-fit IMF in addition to the other fit parameters.
In § 2, we describe the introduction of this additional parameter into EAZY (Brammer et al. 2008), a leading photometric template fitting code. The modified fitting routines are then applied to the COSMOS2015 catalog (Laigle et al. 2016a), the largest current multiwavelength photometric catalog, in § 3. At most redshifts, galaxies are best fit with a higher T IM F than in the Milky Way, corresponding to a top-heavier IMF. In § 4, we discuss the results of various tests to determine whether these fits are truly constraining the IMF, rather than other parameters with which the IMF is covariant or even fully degenerate. Finally, several key implications of these results are discussed in § 5. This constitutes Paper I in a series of three related papers, with the remaining implications of these new measurements reserved for Papers II and III. In Paper II, the starforming galaxies and the star-forming main sequence are discussed and Paper III focuses on quiescent galaxies and quenching.

METHODOLOGY
This section describes the parameterization of a family of IMFs and the dataset and fitting techniques used to find the best-fit IMF for galaxies in the COSMOS2015 catalog.

The COSMOS2015 Dataset
In this work, fitting is performed on the COSMOS2015 catalog (Laigle et al. 2016b), the largest and most comprehensive multi-wavelength catalog available. COS-MOS2015 includes 10 broad-band filters spanning NUV (from GALEX), u (from CFHT/MegaCam) B, V, r, i, z++, Y, J, H, and IRAC channels 1 and 2, two narrow band filters (NB711 and NB816), and 12 intermediate bands (IA427, IA464, IA484, IA505, IA527, IA574, IA624, IA679, IA709, IA738, IA767 and IA827) from the COSMOS-20 survey. Most objects are not covered by every filter, with further elaborations on the catalog detailed in Laigle et al. (2016b). 518,403 galaxies in the COSMOS2015 catalog are not flagged as saturated or within bad regions. However, to properly deconvolve the IMF reasonable uncertainties are required in the photometric fitting code, so a V-band SNR-cut of 10 is used, with 139,535 galaxies passing this cut. In Paper III, an alternative K-band cut is considered, finding that the key results are insensitive to the choice of band and merely require a sufficiently high SNR across the spectrum.

Temperature Dependence of the Stellar IMF
The parameterizations used in this work are modified versions of the Kroupa IMF (Kroupa 2001), withm = 0.08M andM = 0.5M empirically determined for the Galactic IMF. Physically, the mass scale scalesm andM are believed to be the minimum masses for fragmentation and collapse, respectively. Both processes depend upon gas temperature. Cloud collapse has long been modeled with a Jeans approximation (Jeans 1902a), and fragmentation has a similar dependence during adiabatic collapse (Low & Lynden-Bell 1976;Bonnell et al. 2007). As a simplification, here it is also assumed that as the cloud collapses its temperature remains invariant so long as it is optically thin (Steinhardt et al. 2020). Jermyn et al. (2018) used these two processes to derive a temperature-dependent family of IMFs, (2) where T g is the gas temperature in star-forming molecular clouds. They argue that physically, bothm and M scale as T 2 g , for independent reasons. The typical Galactic star-forming cloud is at T ∼ 20 K, which is set as the reference temperature T 0 (Schnee et al. 2008). At T g = T 0 = 20 K, this yields the standard Kroupa IMF. At higher temperatures, this produces a top-heavier IMF, and a lower temperatures, a bottomheavier one (see Fig. 1). Figure 1. A Kroupa initial mass function normalized by total mass with described temperature dependence at T = 10 K [red], 20 K [orange], and 40 K [blue]. Note that raising the temperature leads to a more top-heavy or more bottomlight distribution and lowering the temperature leads to a more bottom-heavy distribution.
It should be noted that the specific dependence ofm on temperature requires solving the energy balance and opacity with a variety of different assumptions (e.g., opacity laws, the ambient temperature of surrounding medium, the angular rotation of the cloud, the ratio of turbulent pressure to total pressure, the strength of magnetic fields), resulting in a wide range of scaling relations, includingm ∝ T (Hopkins 2012),m ∝ T 3/2 (Jeans 1902b),m ∝ T 2 (Steinhardt et al. 2020), andm ∝ T 5/2 (Chabrier et al. 2014). Thus, some attention should be given to the interpretation of the actual bestfit top-heaviness of the IMF, i.e. T IM F , which using the assumptions of Jermyn et al. (2018) corresponds to the gas temperature, T g . For other temperature dependencies, the translation to T g is straight-forward (Table 1). The IMF parameter will correctly describe the shape of the IMF regardless of the proper association with gas temperature.

Reference
Mass Scaling Tg − TIMF Relationship  One might expect that in practice, molecular clouds in a star-forming galaxy should exhibit a range of temperature, and that therefore the true IMF might be a composite of several different components. Further, older stellar subpopulations might have formed under different conditions, and with a different IMF, than younger ones.
In this work, it is instead assumed that every galaxy has a single IMF, and that the same IMF can be used for the entire stellar population. This is done for two reasons. The first is practical; the IMF from Eq. (2) only requires the introduction of a single additional parameter, T IM F . Thus, photometric catalogs with a limited number of bands can still provide meaningful constraints on the best-fit T IM F . This is similar to other approximations; for the same reason it is also standard in template fitting to fit the stellar population as having a single age, and to fit the galaxy with a single extinction, metallicity, etc., even though in practice galaxies are not monolithic. The other is that physical mechanisms have been proposed which could create similar conditions in star-forming regions throughout a galaxy. For example, Steinhardt et al. (2020) proposed a model in which the star-forming main sequence would be driven by feedback between a temperature-dependent IMF and cosmic ray feedback. In this model, because cosmic rays can permeate the entire galaxy, star-forming regions would all have similar gas temperatures. We evaluate the predictions of this model in § 3.3.
In principle, the methods developed in this work could also be used to test this assertion more directly. In recent works, photometry of resolved portions of galaxies have been fit independently to determine whether the same best-fit properties are consistent throughout the same galaxy (Sorba & Sawicki 2018;Greener et al. 2020;Vale Asari et al. 2020;Fetherolf et al. 2020), albeit with mixed results. The same could be done with the new templates produced here.

Introducing variable IMFs into photometric template fitting
Throughout this work, fitting is performed using the photometric template fitting code Easy and Accurate z phot from Yale (EAZY Brammer et al. 2008; source code can be accessed at http://www.astro.yale.edu/eazy/), which has been shown to produce reliable photometric redshifts when compared against other software (Dahlen et al. 2013;Hildebrandt et al. 2010).
EAZY fits photometry by searching for the closest reconstructed spectrum from a library of synthetic spectra. Specifically, the process begins by constructing a larger basis of 560 synthetic templates are derived using the Flexible Stellar Population Synthesis code (FSPS; Conroy et al. 2009a, Conroy & Gunn 2010 which span a variety of ages, metallicities, star-formation histories, dust attenuation and extinction curves (Table 2). These synthetic spectra follow FSPS convention and includes nebular continuum and line emission features derived from the Cloudy models Byler et al. (2017). Then, a smaller basis of 12 theoretical spectra is constructed to span the larger basis of 560 spectra using non-negative matrix factorisation (NNMF). The NNMF reduction takes a large number of synthetic models and computes a reduced basis of templates which best reproduce synthetic broad-band photometric calibration catalogues (see Brammer et al. (2008) for a complete discussion). Thus, the reduced basis represents the "principalcomponent" spectral templates. This analysis follows convention by utilizing the reduction coefficient already derived and utilized in the established EAZY templates (Brammer et al. 2008).
The assumption is then that linear combinations of the reduced basis spans the entire range of real galaxy spectra, yet with few enough parameters that the best fit can be constrained from a small number of photometric bands. EAZY then finds the best-fit linear combination of these 12 basis spectra to the observed photometry, while accounting for redshift and convolving with filter profiles. Best-fit galaxy properties are then calculated as a (luminosity-)weighted average of the properties corresponding to each template. 8 -50 Table 2. Range of galactic parameters spanned in the initial 560 FSPS spectra. These include the average age of the stellar population, characteristic time for star formation τ , metallicity taken from the tabular Mist isochrones, and the dust extinction Av in magnitudes.
In introducing an additional parameter, the goal is to alter this process as little as possible, so that previous experience using and interpreting EAZY fits will still be applicable. Thus, the only novel part of this approach comes in the initial template definition. The larger basis of 560 FSPS templates is calculated independently for every T IM F in a grid of temperatures. Those are then reduced to a smaller, 12-template basis as before, and then used to find a best-fit linear combination to the photometry. Comparing the residual χ 2 n for each fit allows the determination of a single best fit across all choices of T IM F . with increasing temperature [from 20K to 50 K] shown in progressively brighter hues. Increasing temperature results in larger fluxes in the UV [λ < 4000Å] as expected from the additional massive stars. Note the large variability in the young galaxy with temperature while the older galaxy is invariant to temperature-variations.

Distinguishing a change in the IMF from covariant parameters
The greatest impact of increasing T IM F lies in an increased fraction of high-mass stars (Fig. 2), and there-fore a bluer spectrum. Although a bluer spectrum can also be produced by decreasing extinction, extinction is predominantly exponential, while a change in the IMF produces a more complex shape, with different curvature. The two shapes are therefore distinct, so that a change in IMF allows a better fit.
Often, the best-fit model using an incorrect shape might be expected to be most incorrect at both extremes. However, flux at wavelengths longer than ∼ 10 4Å has only a small contribution from starlight. As a result, a change in the IMF has negligible effect, since dust and metallicity are fit independent of star formation when using EAZY. Although the true difference is far more significant on the blue end, a best-fit model with a 20K IMF will typically reduce the total χ 2 by varying the extinction (i.e., less dust than is physically present to compensate for a top-heavier IMF). Therefore, the improvement at high T IM F instead lies predominantly at both extremes of the flux most sensitive to starlight.
In order to confirm this, the residuals of (synthetic) model and observed spectra can be compared. For real spectra, the residuals are calculated by subtracting the best-fit model template (convolved over filter-profiles) from the observed photometric flux and dividing by the reported standard uncertainty in each filter. For model spectra the methodology is identical, with the "observed" photometric flux being a convolution of filterprofiles with a synthetic FSPS model spectrum. The model FSPS spectrum is that of a star-forming galaxy at T IM F = 32K and z = 1.35. The fit was performed using uncertainties in each band that are typical of the overall catalog, even though the model FSPS fluxes were used precisely, with no artificial scatter. The resulting residuals exhibit a similar shape to the actual residuals for objects at 1.2 < z < 1.5 with best-fit T IM F from 30-35K (Fig. 3). Thus, the improvements in fitting for both model and observed photometry are predominantly at the edges of the starlight-sensitive portions of the SED. Such a result is consistent with the interpretation that these objects yield superior fits for T IM F = 30 − 35K than for T IM F = 20K.

Temperature grid and outliers
Fits are performed over a temperature-grid of IMFs, spaced every 1K from 8K to 50K. Most objects are wellfit within a significantly tighter interval of temperatures (typically from 20 to 40 K). However, this analysis still presents the larger interval to emphasise the tight and singular relationship between redshift and T IM F . The lower limit is 8K, as the steepest portion of the IMF then applies for M > 0.08M . Below approximately 0.08M Figure 3. Median difference in χ between the best-fit templates at TIMF and 20K for objects at 1.2 < z < 1.5 with best-fit TIMF from 30-35K as a function of rest-frame wavelength (orange). Shading indicates lower and upper quartiles. The difference is shown as a running boxcar with width constant in log-space and set equal to the wavelength. The primary improvement in fit lies in the UV and NIR, likely due to a difference in shape between (power-law) extinction and the more complex effects of a change in TIMF . The fit residuals between a TIMF = 32K FSPS template and a bestfit reconstructed spectrum at TIMF = 20K are shown for comparison (blue), exhibiting a similar shape the objects are no longer stars in the sense that fusion can no longer take place, so to create an even bottomheavier IMF a different IMF shape would be required (as suggested in Bonnell et al. 2007). The upper limit is arbitrarily set at 50 K as no significant population of objects is fit from 45-59 K. With a 60 K upper limit, no significant population would be fit between 50-59 K. At every temperature within the grid, EAZY computes the best-fit linear combinations of templates, a best-fit photometric redshift z phot and a goodness of fit, χ 2 . The best-fit temperature is the location of the global minimum within the χ 2 -temperature landscape.
One concern with this gridded approach is that if there are many local minima and a coarse grid, it is possible that the minimum χ 2 n may not lie at the global minimum. This is a known issue more generally in photometric template fitting, and one of the reasons that EAZY chooses to fit using an eigenbasis rather than gridding in all parameters. Fortunately, for most objects in the COSMOS2015 catalog, the χ 2 n surface typically indicates only one local minimum, which is also the global minimum (e.g., Fig. 4).
For 60% of all objects in the catalog there is a single local minimum within the χ 2 -temperature landscape, which is also the global minimum. For 40% of the catalog, there are several local minima or the only minimum lies at bounds of the temperature range (at 8K or 50K). A closer examination of these cases reveal that they typically have poorly-constrained photome- try with a low SNR. With the SNR-cut instituted above around 80% have a single local minimum within the χ 2temperature landscape 1 . Furthermore, the number of objects fit at the temperature bounds is tightly correlated with the number of photometric bands used for template-fitting. Decreasing the number of photometric bands used yields a higher fraction of outliers. Conversely, objects observed in a large number of bands are often well-constrained. Objects with poorly-constrained photometry yield poorly constrained T IM F , so for the remainder of this work objects with no clear local minimum are cut from the sample as low quality. As might be expected, this quality cut affects a greater number of objects at lower redshift but a greater fraction of objects at higher redshift (Fig. 5).
A possible concern is that cutting outliers in T IM F may bias towards younger galaxies. Constraints on the IMF are stronger for younger stellar populations than for older ones, because changes in the IMF lead to larger changes in the resulting spectral energy distributions (SEDs) for younger galaxies (Fig. 2). As might be expected, in both cases, increasing temperature leads to more massive stars and thus more emission in the UV, while emission in the optical and IR remains unchanged. However, the older galaxy displays less variability with temperature because massive stars have shorter lifetimes, so most of these massive stars are no longer present.

IMF TEMPERATURES FOR THE COSMOS2015 CATALOG
Using the techniques described in the previous section, a modified version of EAZY was run to find a best-fit IMF and other galactic parameters for every object in the COSMOS2015 catalog with sufficient signal to noise in the V-band. As described in § 2.1, because constraining the IMF involves fitting for an additional parameter, doing so successfully requires a stronger quality cut than in previous COSMOS2015 analysis. Of the 518.403 galaxies objects for which Laigle et al. (2016a) were able to constrain the redshift, only 139.535 had sufficient signal to noise to constrain the IMF. This cut is evaluated in more detail in 4.3 in order to confirm that it does not introduce an additional bias.

A Redshift-IMF Relation
Perhaps the most striking result is that nearly every galaxy in the COSMOS2015 sample is best fit not using a Galactic IMF, but using one at T IM F > 20 K, or topheavier than the Milky Way (Fig. 6).
Further, there is a strong correlation between redshift and the best-fit T IM F , with the characteristic T IM F rising from ∼ 25 K at z ∼ 0 to ∼ 35 K by z ∼ 4. If T IM F is indeed an indicator of conditions in star-forming molecular clouds, it appears that these clouds have similar conditions across most galaxies at the same redshift. Further, those conditions change with redshift, implying that there is a cosmic synchronization to key aspects of galactic evolution.
These results contradict the current assumption that the IMF is universal and constant, so that a Galactic Figure 6. Left: Best fit temperature from 10-50 K versus lookback-time from a sample of 139,535 COSMOS2015 galaxies with SNR>10 in the V-band (Laigle et al. 2016a). At each redshift, the distribution is individually normalised, in order to emphasise the temperature distribution at all redshifts. With increased redshift fewer galaxies are fit at lower temperatures. Right: Boxcar smoothed mean with standard deviation of best-fit gas temperature at different lookback-times [with mean determined from objects in 2 Gyr width age bins and not including galaxies fit at the bounds of temperature range]. The mean temperature increases from ∼ 28 K to ∼ 36 K from present to 12 Gyr, while the spread decreases.
IMF can be used to analyze all galaxies at all redshifts. Thus, one significant implication is that photometric observations of most galaxies may need to be reinterpreted. This might change key inferred properties, including not just stellar masses and star formation rates but also extinction, age, metallicity, etc., which may have been incorrectly inferred for most galaxies. This is discussed in more detail in § 3.4. The specific implications for the star-forming main sequence are discussed further in Paper II, and for quiescent galaxies and quenching mechanisms in Paper III.
Finally, although gas temperatures have historically not been possible to observe for most galaxies, dusttemperatures present an independent probe of the trend seen in Fig. 6 and have been estimated for a broader ensemble of evolving galaxies. The inferred characteristic dust temperature of galaxies on the star-forming main sequence decrease from approximately 35 K to 25 K from z = 2 toward the present time (Magnelli et al. 2014;Magdis et al. 2017). If the temperatures of dust and cold gas clouds are similar, then one would expect a similar trend in gas temperature, which is consistent with the observed evolution in T IM F . It should be noted that the two cannot be compared directly, since the luminosity averaged dust temperature is dominated by hotter dust, while the gas temperature is mass-averaged as it is inferred from the IMF. However, despite the discrepancies, T IM F measurements are entirely independent of dust temperature estimates, so it is striking that the two produce consistent results.

Physical Populations
Galaxies are known to fall into subpopulations with distinct properties. Perhaps the strongest separation is between star-forming and quiescent galaxies, a classification which can easily be determined from color for most galaxies (Williams et al. 2009;Arnouts et al. 2013). Cool gas in star-forming galaxies is used to produce stars. Various mechanisms have been proposed for quiescence, either depleting the cool gas or producing conditions under which cool gas exists but cannot efficiently produce stars. It is therefore natural to wonder whether both star-forming and quiescent galaxies exhibit similar IMFs.
Using UVJ selection, the galaxies in the COS-MOS2015 catalog were divided into star-forming and quiescent populations. The distribution of best-fit T IM F for the two populations differs (see Fig. 7). The continuous range of temperatures at a given redshift seen in Fig. 6 consists of a combination of cooler quiescent galaxies and a range of slightly hotter star forming Figure 7. Box-car smoothed mean and spread of inferred temperature for galaxies in age bins of width 2 Gyr. Quiescent galaxies display Milky Way like IMFs and lie on the lower range of temperatures spanned by star-forming galaxies. Well-measured quiescent galaxies are sparse beyond 9 Gyr, with significant numbers of star-forming galaxies observed up til 12 Gyr.
galaxies. Thus, a selection based on T IM F provides a surprisingly effective separation of the two classes. The coolest galaxy at any redshift is typically quiescent.
An examination of best-fit temperatures for synthetic spectra (see § 4.3) implies that the total uncertainty on individual temperature measurements is approximately 10K. Thus, many individual objects cannot be selected as quiescent or star-forming from their best-fit T IM F alone. However, the uncertainty on the mean of a population of galaxies is naturally far better constrained given the Central Limit Theorem. For 0 < z < 1, there are 7.800 quiescent galaxies and 68.000 star forming galaxies, so the star-forming galaxies mean T IM F is more than 30σ higher than the quiescent sample. However, while the statistical uncertainty on an ensemble is small, it is difficult to conclusively show that this separation truly comes from a difference in temperature. Because the EAZY template basis only spans a limited subspace of possible galaxy spectra, an alternative explanation could be that typical, e.g., quiescent spectra lie outside that subspace, while typical star-forming spectra are better modeled by that limited basis. However, these results are consistent with the naive prediction that a less active galaxy should be colder due to less stellar light and lower cosmic rays production. This is also consistent with measured dust temperatures, which tend to increase towards higher SFR at fixed stellar mass and redshift Magnelli et al. (2014).
A smaller subpopulation are the the ultra-luminous infrared galaxies (ULIRGS), hypothesized to be caused by mergers of gas-rich galaxies (Lonsdale et al. 2006). The typical best-fit IMF of ULIRGs (Fig. 8) is nearly independent of redshift within the range of lookback times investigated from 6 to 12.5 Gyr. In that respect, ULIRGs behave differently than typical starforming galaxies, which have hotter gas towards high redshift. This might reinforce the belief that ULIRGs are governed by a distinct process from those which drive the star-forming main sequence, and thus typical starforming galaxies. This redshift-independence of temperature has also been measured in dust-temperatures of ULIRGs (Béthermin et al. 2015), again suggesting that gas and dust temperatures are correlated, either because they are in equilibrium or due to feedback processes. Figure 8. Box-car smoothed mean and spread of inferred temperature for ULIRGs with age bins of width 2 Gyr. The average temperature is consistent with being independent for lookback-times ranging from 6 to 12.5 Gyr.

Comparing Temperature Evolution with Theory
Although previous template fitting techniques have generally assumed a Universal IMF, there have been several solid theoretical reasons to expect that this might not be the case. It is therefore natural to consider whether the best-fit T IM F are consistent with what might have been expected theoretically.
Within the Milky Way, the typical temperature of cold gas in molecular clouds is set primarily by stellar radiation. Different molecular clouds may have different temperatures, depending upon their proximity to luminous stars. The most common case will be temperatures set by low-mass stars, as the Galactic SFR is relatively low, so the probability of a high-mass star in close proximity to a molecular cloud is small. However at higher SFR, as was the case for much of the Milky Way's history, the dominant contribution may instead come from cosmic rays, which permeate the entire galaxy and therefore can heat clouds uniformly (Papadopoulos et al. 2010;Papadopoulos 2010b). At z > 6, beyond the redshift domain discussed in this work, the CMB temperature also exceeds that of typical Galactic molecular clouds (Jermyn et al. 2018). In both cases, clouds would have similar temperatures throughout the galaxy. In this work, a single characteristic temperature is fit due to the limited information available, so that the resulting bestfit parameters follow the standard convention of treating galaxies as monolithic, with one SFR, T IM F , extinction curve, metallicity, etc.
In the Steinhardt et al. (2020) model of galaxy evolution the feedback processes between cosmic ray heating from the deaths of massive stars and the temperature dependence of the IMF result in an attractor solution with a common T IM F . They predict that T IM F should be nearly constant for star-forming galaxies at fixed redshift, and declining towards lower redshift. Furthermore, quiescent galaxies should have lower T IM F than star-forming galaxies, but should not continue to drop past quiescence, as the dominant heating source at low SFR is instead starlight from low-mass stars. These predictions are consistent, at least qualitatively, with the results found here, but there are sufficient free parameters that Steinhardt et al. (2020) were unable to make robust quantitative predictions with their model.
There are also several other observations which hint at the same conclusion. Although absolute calibration is difficult, so that inferred dust temperatures vary quantitatively between models, a corresponding attractor dust temperature and decrease in dust temperatures toward low redshift is seen on the star-forming main sequence (Casey 2012;Magnelli et al. 2014;Béthermin et al. 2015;Magdis et al. 2017;Cortzen et al. 2020). This, too, would suggest that that the IMF typically becomes more top-heavy at larger redshifts.
Notably, this is one of many other tight correlations known in galaxy evolution including the star-forming main sequence, a tight correlation between star formation rates and the existing stellar mass at fixed redshift (cf. Peng et al. (2010); Speagle et al. (2014);Tomczak et al. (2016)). Such a correlation might imply a universality to star formation with a common evolutionary track , but alternatively could be explained by averaging over a large number of shorter bursts (Kelson 2014). The consistency of T IM F is clearly consistent with the former explanation. For the latter, it would appear to be inconsistent with bursts on a galaxy-wide scale, for which T IM F would vary significantly. However, a galaxy composed of several subpopulations, which are individually either starbursting or not at any given time, could still produce a nearlyconstant T IM F when the entire galaxy is fit simultane- Figure 9. Best fit stellar mass (M) and star formation rate (SFR) from EAZY template-fitting for individual galaxies as a function of TIMF . The dashed line is the median with shading indicating first and third quartiles. Any observed galaxy will only have a single best-fit IMF, but the stellar mass and SFR may be biased by an erroneous choice of IMF. Stellar Mass decrease monotonically with increasing temperature as we would expect from a top-heavier IMF, while SFR increases slightly beyond 40 K.
ously and treated as monolithic. Testing this hypothesis will require fitting resolved subpopulations of nearby star-forming galaxies in order to determine whether they all have the same T IM F .

Reinterpreting Galactic Properties
Nearly every galactic property inferred from photometric template fitting is sensitive to the low-mass stellar population which is not luminous enough to be seen, and therefore depends upon the assumed IMF. Thus, changing the IMF will also change these inferred properties. Following the standard EAZY parameter reconstruction (Brammer et al. 2008), the galactic properties of individual templates can be propagated by a luminosity-weighted sum of the basis templates into the properties of any photometric observations. A first approximation of this effect can be seen in Fig. 9, where the best-fit parameters are shown for identical sets of photometric colors for different T IM F .
Perhaps the most straightforward changes are those in stellar mass and IMF for small variations in the Galactic IMF (T IM F ∼ 20 K). Bottom-heavier IMFs (lower T IM F ) produce larger inferred M * and SFR, as more stellar mass is 'hidden' in small stars. The luminosity and mass are respectively dominated by the upper and lower end of the stellar mass distribution, so a bottomheavy IMF decreases the luminosity relative to the mass. Thus, for the same measured luminosity, a greater mass is inferred. Correspondingly, a top-heavier IMF will have lower stellar masses and SFR than at T IM F = 20 K. Given the evolution evident in Fig. 6, this implies that previous analyses may have systematically overestimated stellar masses, particularly towards high redshift, with the true stellar masses being around 40 % less than predicted for a typical z = 2 galaxy. Naturally, the exact bias in inferred properties will depend on the parameterization of a variable IMF.
Although at T IM F < 30 K, the trends for stellar mass and SFR as a function of T IM F are qualitatively similar, above 30 K they diverge. The SFR is no longer monotonic, and overall exhibits a weaker temperature dependence than M * . Notably, this means that varying the IMF will change the specific star formation rate (sSFR). The implications of this for the star-forming main sequence are described in detail in Paper II.
The counterinuitive increase in best-fit SFR towards higher T IM F is due to covariances with other galactic properties which are fit simultaneously. A higher T IM F produces more high-mass stars, and thus a bluer stellar population. Since the spectral tilt of the underlying photometry has not changed, instead extinction is increased to compensate, and the best-fit combination becomes both higher-SFR and higher-extinction than before. Of course, these covariances also present a significant concern. If these properties were fully degenerate, then it would never be possible to determine the true IMF, SFR, or extinction. Fortunately, they are only partly degnerate, as shown in § 4.2.
Finally, note that although there is a strong dependence between T IM F and other parameters, there is always a single, best-fit combination of T IM F , M * , SFR, etc. The dependence shown in Fig. 9 is useful to estimate the extent to which properties calculated assuming a Galactic IMF might be incorrect. For many highredshift galaxies, a Galactic IMF can be rejected, although for some galaxies, particularly a lower redshift, it cannot be.

ROBUSTNESS OF IMF MEASUREMENTS
Photometric template fitting uses limited information to fit many parameters with high covariance. As a result, quantities inferred for individual galaxies are rarely tightly constrained. The introduction of an additional parameter, T IM F , only adds to the difficulty and weak-ens the constraints. A stronger quality cut is needed to select objects for which T IM F can be measured than in the previous COSMOS2015 catalog, with fewer parameters. In this section, the robustness of these measurements is evaluated in more detail.

Significance and Uncertainty Estimation
The χ 2 landscapes produced by photometric template fitting can be remarkably complex, characterised by a large number of local minima, which are often wide but shallow (Speagle et al. 2016). In addition, because the set of physical models spanned by templates is only expected to approximate the true spectrum, rather than contain it, the errors will not be chi-squared distributed, and therefore the likelihood distribution for T IM F cannot be inferred using a chi-squared test.
Instead, here uncertainties are estimated using a Monte Carlo simulation.
Each band has a wellquantified uncertainty on the flux within that filter. Drawing repeatedly from that distribution and finding best-fit parameters produces a resulting distribution of T IM F indicates a typical T IM F uncertainty of 4 K. The standard deviation in T IM F observed at fixed redshift (see Fig. 6) ranges from 8 K (at z ≈ 0) to 5 K (at z ≈ 3). Importantly, the uncertainty in T IM F for individual galaxies is smaller than the spread around z ∼ 0 (see Fig. 6), suggesting that galaxies in the local Universe are not consistent with a single T IM F . However, at high redshift the spread in temperature is consistent with being derived from a population at a single temperature. One possible explanation, described in Paper II, is that star-forming and quiescent galaxies typically lie at different T IM F . If so, a high-redshift population composed nearly exclusively of star-forming galaxies would have a single temperature but one in the local Universe, including significant fractions of both star-forming and quiescent galaxies, would not.

Degeneracies and Covariances
A significant concern in constraining the IMF is the possibility of degeneracies in the spectrum with other galactic properties. In particular, there is in principle a complete denegeracy between the IMF and the star formation history (SFH). An identical observed stellar population can be produced from different IMFs by varying the star formation history. This idea has also been used, assuming a Galactic IMF, to determine extragalactic SFHs (Panter et al. 2003;Sánchez et al. 2019). The most recent SFR is determined from the highest-mass objects, then that population subtracted off and the next-most-recent SFR is determined from the highestmass residual, and so on. However, such a solution is possible for a wide range of IMFs. In principle, a particularly incorrect IMF may require a negative SFR at some times, but in practice, many reasonable choices of IMF produce non-negative results.
Here, this degeneracy is broken by assuming a shape for the SFH. In the standard EAZY methodology, templates are generated using a delayed-tau SFH, with a range of different ages for the stellar population (Brammer et al. 2008). For this single shape of the SFH, it is then possible to uniquely determine a best-fit IMF. Changing the assumed shape of the SFH will, of course, change the resulting best-fit parameters. As an extreme example, one can generate templates using a history consisting only of a single starburst, producing a single stellar population. However, applying these SSP-templates to observed photometry yields a best-fit T IM F similar to those with the delayed-tau model (see Fig. 10). Thus, the inferred IMF remains relatively robust regardless of varying assumptions of SFH. Figure 10. Inferred TIMF for 1,000 objects given a basis of 12 representative SSP templates and the EAZY templatebasis of composite delayed τ -model SFH history. Regardless of assumption in star formation history similar temperatures are inferred with a correlation coefficient of 0.5 for objects found with a χ 2 -minima within the temperature range.
An additional concern are the covariances between T IM F and properties such as metallicity and dust. For example, to first approximation, a higher T IM F produces a bluer spectrum for the same physical properties, much as lower extinction would. To investigate this, the same Monte Carlo simulation described in the previous section was used, comparing the distributions of pairs of physical properties (see Fig. 11). As expected, there is large degeneracy between T IM F and M * , as the IMF is a near-direct scaling between the bright end of the stellar mass function and the remainder. Furthermore, an increasingly top-heavier IMF is offset by increasing extinction. Indeed, this trade-off between extinction and IMF could produce the characteristic shape of residuals discussed in § 2.4. Finally, there are no significant degeneracies between T IM F and metallicity. However, metallicity is derived from the strength of lines and is poorly constrained from photometry. Therefore, the small statistical variation is not representative of the systematic uncertainties associated with constraining Z for template-fitting (Mitchell et al. 2013).
Additionally, the co-variance between T IM F and the inferred z phot can be estimated by examining the 23.000 objects with additional spectroscopic redshift in COS-MOS2015 catalog (Laigle et al. 2016b). Here, the scatter between spectroscopic redshift and photometric redshift, z phot − z spec , is unchanged with a variable T IM F . Thus, the improvements in the fit does not decrease the accuracy of the photometric redshift. Physically, the photometric redshift in template fitting is typically dominated by breaks in the spectrum, while T IM F is primarily fit from the continuum.
The wide range of best-fit stellar initial mass functions can also be validated by inspecting the spectral energy distributions.
A composite observed SED has been constructed combining the COS-MOS2015 photometry with known spectroscopic redshifts using the EAZY-py implementation for composite SED (https://github.com/gbrammer/eazy-py, Brammer et al. (2008)). The observed photometry is (using filter-profiles and the spectroscopic redshifts) converted to an estimated rest-frame flux. Thus, using the multitude of bands and the wide range of redshifts one can sample a galaxy-populations rest-frame spectrum. These observed SED are not corrected for dust attenuation, which is however included in the fitting procedure for T IM F . As they do not account for dust, composite SED have limited interpretability, but are useful for determining if separations in parameters-space leads to a separation of objects in color-space. The stacked photometry can be seen for galaxies split in two ranges of best-fit temperatures in Fig. 12. Importantly, these populations show different spectral shapes, with more flux in ultraviolet (UV) and bluer part of the optical spectrum for hotter galaxies. Thus, stratifying objects in T IM F leads to a separation of objects in the specific regions of color-space most sensitive to a change in IMF. Frequency of best fit TIMF versus metallicity, extinction coefficient and stellar mass for a single object in COSMOS2015 with 10,000 unique perturbations of observed photometric flux. The perturbed fluxes are generated by sampling from a normal distribution with mean observed flux and standard deviation flux error, as reported in the catalog. TIMF displays respectively negative and positive covariance with the inferred mass and dust extinction, with no clear correlation to metallicity. Figure 12. Composite SED normalised at 5600Å for galaxies with best-fit temperature in 10-20 K range [red] and 35-45 K range [blue] with standard deviation in shading. The central line indicates median SED with shading indicating 16 and 84 percentiles. The cooler galaxies display systematically less flux than the hotter galaxies at all wavelengths of the ultraviolet [λ < 4000Å]. Cooler galaxies look more quiescent than the hotter galaxies, which share spectral features with more star-forming galaxies.

Quality Cuts
A final concern is the possibility that the tight relationship between best-fit T IM F and redshift arises from the cuts applied to the sample, rather than from a true astrophysical relationship. However, this concern may be alleviated by applying the same approach to synthetic photometry of spectra with known T IM F . Here photometry has been constructed using the EAZY tem-plate basis, with scatter characteristic of galaxies at different signal-to-noise from the COSMOS2015 catalog. A best-fit IMF temperature (Fig. 13) is then determined. For the majority of objects in the COSMOS2015 catalog, the photometry is not constrained well enough to determine T IM F , but for the highest-quality objects, the input T IM F can be recovered.
The spread observed between input and inferred temperatures of synthetic photometry is typically about 10 K for the 95th percentile; This is significantly larger than the 5 K spread suggested by the uncertainty from the Monte Carlo simulations (see § 4.1) and the physical spread observed at fixed redshift (see § 6). This reflects a generic issue present in template-fitting codes when the observed galaxies are outside the synthetic models parameter-space. Any fitting would in such cases typically happen at the bound of the parameter-space, which would artificially suppress the variance (Brammer et al. 2008). In future, this could be alleviated with a deeper understanding of galactic spectra and improvements to FSPS modeling, but this analysis remains agnostic, since the T IM F for an ensemble of galaxies is highly significant with either metric of uncertainty.
Additionally, the average best-fit T IM F is consistent with no systematic bias for increasing uncertainties on flux-measurements. However, the fraction of objects best fit at the bounds of the temperature range monotonically increases with increased noise (see Fig. 14). Therefore, the well-measured galaxies are important for constraining the IMF, but not for the overall trend with . For each panel 10,000 synthetic photometry are generated with a known input temperatures ranging from 10 to 50 K. The photometry is scattered with Gaussian uncertainties in each filter. The uncertainties are drawn corresponding to the representative SNR for each band within COSMOS2015 [from left to right the panels show the SNR from the median, 80th, 95th and 99th percentile]. With increasing SNR, the inferred temperature increasingly matches the underlying input temperature. However, even for the 99th percentile a systematic bias between input and inferred exists, where especially bottom-heavier IMFs remain hard to distinguish.
redshift seen in § 3. The inclusion of temperatures outliers does not influence the trends presented in this paper; The redshift evolution and the reinterpretation of galactic properties displays the distribution and dependence of T IM F , while the differences between ULIRGs, star-forming and quiescent galaxies remains valid with or without outliers. Typically, objects fit at the bounds are ill-constrained because of limited coverage in optical and UV. Thus, because the inferred temperature may not represent the underlying IMF, we recommend caution when interpreting these outlying T IM F .

DISCUSSION
This work introduces an additional parameter, T IM F , into photometric template fitting. This allows a selection of the stellar initial mass function (IMF) from a theoretical family of possible IMFs, with a higher T IM F producing a top-heavier stellar mass distribution. A Milky Way-like IMF (here, a Kroupa IMF), which previously has typically been assumed to be Universal, is part of this family and corresponds to T IM F = 20 K.
Fitting the COSMOS2015 catalog with this additional parameter produces three key results: • Top-Heavier IMFs: Most galaxies are best-fit with T IM F > 20 K, typically ranging from 25-40 K. This means that most galaxies should have Figure 14. Fraction of objects best fit at the bounds of temperature range as a function of SNR in V band . As noise decreases the fraction of outliers decreases monotonically. Note in paper III, it is shown that a SNR-cut in K band yields equivalent results. Objects still best fit at the temperature bounds at high SNR typically have little coverage within UV and optical bands.
top-heavier stellar populations than previously expected. One likely consequence is that stellar masses(M * ) and star-formation rates (SFR) have previously been overestimated, with greater overestimation associated with galaxies with higher

Limitations From Assuming a Monolithic IMF
In order to introduce only a single additional parameter, in this work galaxies are modeled as having a single IMF for their entire stellar population. This assumption is very likely unphysical. If T IM F evolves with redshift, as found here, then the stellar population in a galaxy would be comprised of components produced from different IMFs. Further, unless the dominant contribution to galactic heating truly permeates the entire galaxy, different molecular clouds should have different gas temperatures. This issue might be less critical, since the coolest clouds should be the sites of the bulk of the star formation.
In some respects, the same is true of other parameters in photometric template fitting. Galaxies in the COSMOS2015 catalog are described as having a single, luminosity-weighted, age, extinction, and metallicity (or, in some cases, probability distributions for that single parameter; Laigle et al. 2016a). Fitting portions of galaxies which can be resolved into multiple components often results in these parameters varying between those components (Greener et al. 2020;Fetherolf et al. 2020).
EAZY fits galaxy photometry as a linear combination of a basis which includes components of different age, extinction, and metallicity. Thus, at least in principle, a galaxy with distinct subpopulations could be well-fit with this technique. However, the 12 basis vectors are insufficient to span the full space of subpopulations.
Introducing the additional parameter of T IM F would not have allowed 12 eigenspectra to span the full space spanned by the standard EAZY template basis. Expanding the size of the basis would have required too many additional measurements to constrain the best fit, so that high-redshift objects could not be modeled. As a result, instead T IM F is treated differently than these other parameters, fit from a grid of single values. A grid-like approach has been used successfully in other template-fitting codes, which model the full stellar population with a single set of parameters. However, this means that stellar populations produced from multiple IMFs will be fit with a single, luminosity-averaged, IMF. As a galaxy becomes quiescent, even without new star formation the best-fit T IM F could vary, as the luminosity-averaged IMF transitions from dominated by the youngest population to older ones.

Adding Additional Parameters
One solution would be adding a second parameter, corresponding to an evolutionary history of the IMF. It would be straightforward to this parameter to span a reasonable range of shapes for that history. However, it should be noted that the addition of one parameter required a significantly more stringent quality cut (Fig. 5) than in the full COSMOS2015 catalog. At present, only a small fraction of objects in COSMOS2015 can be wellmodeled, and other photometric surveys typically have fewer available bands. Adding a second new parameter would have a similar effect, and it might become impossible to fit any objects at, e.g., z > 2.
If sufficient information is available to add a parameter, perhaps a more useful choice would be the dust extinction curve. As with the IMF, currently a single extinction curve is typically assumed for fitting an entire photometric catalog, with only the normalization determined by the fit. However, there is a strong argument for variability with measurements of local galaxies finding a range of different extinction curves (Gordon et al. 2003;Salim & Narayanan 2020).
There is also the possibility of a strong covariance between the best-fit extinction curve and the best-fit IMF shape. Thus, modifying template fitting codes to investigate this possibility would provide an important test of the robustness of the results in this work. This likely requires a substantial rewrite to either EAZY or a different template-fitting code, and is thus beyond the scope of this work.
Another significant parameter not currently included is the shape of the star-formation history (SFH). Although there is theoretically a complete degeneracy between the SFH and the IMF, in practice different reasonable choices of SFH produced similar T IM F (Fig. 10).
Ultimately, despite the limitations of adding only a single parameter to describe the range of possible stellar initial mass functions, the best-fit IMF is well constrained from ∼ 30-band, multi-wavelength photometry. For most galaxies, this produces an IMF different than that of our own Galaxy. Because nearly every quantity derived from photometric template fitting is sensitive to the IMF, finding ways to include a range of IMFs should be considered necessary for future template fitting pipelines. The resulting implications for models of star-forming galaxies (Paper II) and quiescent galaxies (Paper III) are discussed in companion papers.