Thermal Regulation and the Star-Forming Main Sequence
Charles L. Steinhardt Adam S. Jermyn Jackie Lodman


firstpage–lastpage
*
    Mark C. Wyatt
    December 30, 2023
=====================



We argue that the interplay between cosmic rays, the initial mass function, and star formation plays a crucial role in regulating the star-forming "main sequence".  To explore these phenomena we develop a toy model for galaxy evolution in which star formation is regulated by a combination of a temperature-dependent initial mass function and heating due to starlight, cosmic rays and, at very high redshift, the cosmic microwave background. This produces an attractor, near-equilibrium solution which is consistent with observations of the star-forming main sequence over a broad redshift range. Additional solutions to the same equations may correspond to other observed phases of galaxy evolution including quiescent galaxies. This model makes several falsifiable predictions, including higher metallicities and dust masses than anticipated at high redshift and isotopic abundances in the Milky Way. It also predicts that stellar mass-to-light ratios are lower than produced using a Milky Way-derived IMF, so that inferences of stellar masses and star formation rates for high redshift galaxies are overestimated. In some cases, this may also transform inferred dark matter profiles from core-like to cusp-like.





galaxies: star formation - galaxies: stellar content - galaxies: luminosity function, mass function - cosmology: cosmological parameters - cosmology: cosmic microwave background radiation - (ISM:) cosmic rays








§ INTRODUCTION


One of the most remarkable recent discoveries in galaxy evolution has been the surprising similarity in the properties of evolving galaxies on macroscopic scales, to the point that it has been possible to observationally describe the properties of a `normal' evolving galaxy.  It was already known that elliptical galaxies follow a series of scaling relations between radius, velocity dispersion, and surface brightness known as the `fundamental plane' <cit.>.  More recently, star-forming galaxies at constant stellar mass and redshift have been shown to exhibit only a narrow range of star formation rates (SFR), a result that has come to be known as the star-forming main sequence <cit.>.  There is a corresponding similarity for quasar accretion, relating luminosity to virial black hole mass and redshift <cit.>.

This similarity is particularly surprising because, in modelling star formation and quasar accretion, it would seem that every process which affects gas availability, density, composition, dynamics, temperature, and pressure must be considered.  
Furthermore, spatial degrees of freedom could well be relevant, even though there have been successful modelling efforts which neglected these <cit.>.  

Although intuitively this complexity might be expected to lead to a chaotic solution in which every galaxy has an independent evolution driven by the local environment, attractor solutions (i.e., solutions in which different initial conditions lead to similar behaviour) could also exist <cit.>.  Main sequence stars are locally chaotic, but feedback driving stars towards hydrostatic and thermal equilibrium results in an attractor solution for their macroscopic behaviour <cit.>.  Thus, it can be straightforward to predict the lifetime of the Sun although it is difficult to predict individual solar flares <cit.>.

We are therefore motivated to search for a similar equilibrium condition that might drive evolving galaxies toward a common path. 
In recent work, the authors considered the evolution of the stellar initial mass function (IMF) with temperature in star-forming clouds <cit.>, showing that a galactic IMF will be bottom-lighter than the Milky Way IMF when the temperature in star-forming regions exceeds the ∼ 20 K of Milky Way star-forming clouds.  Indeed, both dust <cit.> and gas <cit.> temperatures in evolving galaxies are measured at > 20 K, implying that the IMF should be bottom-lighter than the Milky Way-derived IMFs commonly used in photometric template fitting of their properties <cit.>.  An increase in temperature also means that fewer molecular clouds can collapse and thus decreases the overall SFR, a process that has been at the heart of models describing feedback between the evolution of the central supermassive black hole and its host galaxy <cit.>.

Further, unlike the Milky Way, at high SFR, the dominant contribution to gas temperature likely comes from cosmic rays generated in the deaths of the most massive main sequence stars <cit.>.   Thus, cosmic rays provide a feedback loop: an increase in cosmic ray generation will decrease SFR and provide a bottom-lighter IMF, which in turn changes the rate which new cosmic rays are produced.  Depending upon the precise temperature-dependence of the SFR and IMF, both equilibrium and runaway solutions might be possible, as are more complex chaotic solutions.  

Considering cosmic rays as the primary mechanism for driving temperature is natural for two additional reasons.  Unlike photons and most other products of main sequence stars and stellar death, cosmic rays permeate gas throughout the entire galaxy, contributing to temperatures even in the coolest star-forming clouds.  Moreover, despite their likely importance, cosmic rays have proven difficult to include in current numerical simulations, so that a closer consideration of their properties is natural.  Indeed, their complexity has led to several other proposed cosmic ray feedback mechanisms, most notably cosmic winds <cit.>

We outline the way in which star formation, the initial mass function, and cosmic rays interact in  <ref>.
In  <ref>, we show that this interaction results in stable attractor solutions at fixed gas mass. We then combine these interactions in  <ref> with a basic numerical model of galaxy evolution which tracks the evolution of a halo, infalling gas, and a stellar population.  In some regimes the attractor is dominated by cosmic ray heating, and the properties of this track are compared with the star-forming main sequence in  <ref>.  Whether this idea can be developed into a reasonable toy model for galaxy evolution is considered in  <ref>.



§ MOTIVATION


We are motivated by three main considerations, each of which we will describe briefly below.  In combination, these suggest the development of the model proposed in this work.



 §.§ Temperature Dependence of Star Formation


One of the most important processes governing star formation is the competition between gravity, driving the production of compact objects, and thermodynamics, which drives gas to expand rather than contract[Magnetic pressure and radiative effects also likely play a role, but our focus here is on the balance between temperature and gravity.].  Therefore, a change in gas temperature will necessarily affect the balance between these two forces, changing the nature of galactic star formation.  This balance is often expressed in the framework of calculating a Jeans mass <cit.>, although there are several other effects which cannot be neglected in a full treatment.  

An increase in gas temperature will increase the Jeans mass, and therefore fewer molecular clouds will collapse.  Those which do collapse will also collapse more slowly.  As a result, galactic star formation rates (SFR) will decrease with increasing gas temperature.

At the same time, the clouds which do collapse at higher temperatures will also produce a different stellar mass distribution, since the fragmentation process is temperature-dependent <cit.>.  In <cit.>, we provide a model of this temperature dependence and argue that an increase in stellar mass temperature produces a bottom-lighter (or, equivalently, top-heavier) stellar initial mass function.  At increased temperature, galaxies will have a lower star formation rate, but the stars produced will be more massive, and in some cases this could even lead to an increase rather than a decrease in the production of O stars despite a lower overall SFR.



 §.§ High-redshift Galactic Temperatures


Although gas in the Milky Way is typically around 20 K <cit.>, observations of higher-redshift star-forming galaxies indicate that their gas is likely at higher temperatures.  Although it is difficult to measure the gas temperature directly, dust temperatures are more approachable.  <cit.> find that dust temperatures along the star-forming main sequence range from 25 K to 40 K, with increasing temperature towards higher redshifts.  In addition, higher SFR at fixed redshift leads to a temperature increase, as might be expected due to the combination of higher luminosity and higher supernova rates.

It should be stressed that reasoning from luminosity-averaged dust temperatures to gas temperatures specifically in star-forming clouds introduces several major assumptions which are poorly justified.  It is very likely that star-forming clouds take on a variety of temperatures, just as observations of resolved galaxies indicate that star formation at any given time is typically concentrated only in some regions of star-forming galaxies (cf. <cit.>).  Nevertheless, the most reasonable interpretation of the available observational evidence is that star-forming galaxies sustain gas temperatures higher than those in the Milky Way.  

At z ≳ 6, the CMB temperature drives gas clouds to T > 20 K, and even if star-forming galaxies are hotter than the Milky Way, at some sufficiently high redshift, the CMB will dictate gas temperatures.  We consider this regime in building initial conditions for the model developed in  <ref>, but there are currently no measurements of gas or dust temperatures at these extreme redshifts.



 §.§ Cosmic Ray Temperature Driving


In the Milky Way, the ∼ 20 K typical temperature of cold gas in molecular clouds <cit.> is driven primarily by stellar radiation.  Cosmic rays alone would only drive the temperature to ∼ 5 K <cit.>, so their contribution is negligible.  

However, at higher SFRs a corresponding increase in the supernova rate will result in a higher cosmic ray density.  A galaxy with an SFR of 100 M_⊙/yr will produce cosmic rays at 100 times the Milky Way rate.  The corresponding temperature increase is complex, with simulations estimating an average temperature of[ To obtain this form, it was assumed that the cosmic ray ionization rate per H_2 molecule is directly proportional to ρ_ CR in equation (12) of <cit.>. Following their estimates for the former, 5× 10^-17 s^-1 was assumed for the Milky Way.  For similar reasons, an H_2 density of 10^4 cm^-3 Milky Way ISM temperature of T_ Milky Way≈ 20  K were used.]

    T̅_ CR≈ 0.32((5n̅^1/2ρ̅_ CR + 0.28 n̅^3)^1/2 - 0.53n̅^3/2)^2/3,

<cit.>, where T_ CR is the contribution of cosmic rays to the temperature, n is the density of H_2 and ρ_CR is the average cosmic ray density.
The over-bar indicates normalisation with respect to the Milky Way, i.e., all over-barred temperatures are divided by the Milky Way ISM temperature of 20 K.
At values that are likely typical for z ∼ 1-6 star-forming galaxies, this will lead to a temperature dominated by cosmic rays rather than starlight or the CMB.

Critically, cosmic rays are predominantly produced by supernova shocks <cit.>.
Therefore, the details of the stellar population determine the rate at which cosmic rays are generated, while the bulk properties of the galaxy determine their transport.  Thus, the rate at which new cosmic rays are generated is in turn related to the current ρ_ CR.  



§ A COSMIC RAY-DRIVEN MODEL

 

These three components suggest a model of the galactic star formation in which the CMB, cosmic rays and stellar radiation conspire to set the mean temperature gas clouds.  That temperature, in turn, determines the rate of star formation and the initial mass function, which both serve to set the density of cosmic rays and the intensity of stellar radiation.  These effects combine to form a feedback loop, which is depicted diagrammatically in Figure <ref>.



This loop is of particular interest because it directly relates the star formation rate and distribution at any point in time to the historic stellar population.
It is also comparatively simple yet potentially results in sufficiently complex behaviours as to be physically interesting.
Our aim is therefore to model and understand the phenomenology of this cycle.

Although the various components of this cycle are complex, there are relatively few observationally well-constrained properties of evolving galaxies.  At higher redshifts, these are generally averaged over the entire galaxy, even though higher-resolution and lower-redshift observations find that galaxies often have complex structure with different, e.g., star formation rates in different regions <cit.>.  

As a result of these observational constraints, it is only possible to produce a meaningful model if there are very few parameters.  Thus, here we will treat the cycle as acting in a spatially homogeneous galaxy.  We further neglect all temporal variation in metallicity, and treat cosmic rays as being absorbed by the interstellar medium on a fixed time-scale.  Even with these considerable simplifications we require several additional ingredients.



 §.§ IMF


We obtain the temperature-dependent IMF from <cit.>, who find

    ξ(M,T̅) = dN/dM(T̅)∝
    M^-0.3,     M < 0.08 M_⊙ f(T̅) 
    
    M^-1.3,     0.08 M_⊙ f(T̅) < m < 0.50 M_⊙ f(T̅)
    
    M^-2.3,     0.50 M_⊙ f(T̅) < m,

where

    f(T̅) = T̅^2.

The proportionality factors in equation (<ref>) are such that ξ is continuous and 

    ∫_0^∞ Mξ(M, T̅) dM = 1 M_⊙.

Hence ξ dM is the number of stars in a mass bin of width dM given 1 M_⊙ of star formation.
Note that this choice requires that the proportionality factors depend on T̅.



 §.§ SFR


Next we must determine the SFR as a function of temperature.
Note that the mass distribution of clumps in the Milky way is 

    n(M) ∝ M^-1.60 ± 0.06

where clumps are massive enough (10M_⊙ up to 3000M_⊙) to allow measurement <cit.>.  We take this distribution to hold for all M.  It is further assumed (cf. <cit.> that there is negligible temperature dependence to this mass distribution.

Because Jeans timescale for collapse is independent of mass <cit.>, this implies that the star formation rate at a given temperature is proportional to the number of clumps that exist above the direct collapse Jeans mass, M_J ∝ T^3/2. Thus, we find that

    SFR∝∫_M_J^∞M^-1.60 ± 0.06 dM∝ M_ gas T^-0.90 ± 0.09,

where M_ gas is the total gas mass of the galaxy.
Hence, scaling from the Milky Way, we obtain

    S̅F̅R̅ = M̅_ gasT̅^-0.90± 0.09.




 §.§ Cosmic Rays


Once formed, stars of mass M evolve and ultimately die over a timescale of

    τ≈ 10^10 yr(M/M_⊙)^-5/2.

The most massive of these new stars generate cosmic rays when they explode.
The overall cosmic ray generation rate in the Milky Way is well-approximated by assuming that each supernova converts an approximately fixed fraction q ≈ 1% of its energy into cosmic rays <cit.>.
This is in good agreement with theoretical models of supernova shocks <cit.> and so we assume that the cosmic ray generation from a single star is

    CRR(M) = q E_SN(M),

where CRR is the cosmic ray generation rate and E_SN is the supernova energy associated with a star of mass M.

The majority of the galactic supernova energy budget arises from core-collapse and pair instability supernovae, with of order E_ CC/PI≈ 10^53erg emitted in each such event, somewhat independent of the progenitor mass <cit.>.
By contrast, each Type Ia event emits of order 10^51erg <cit.>, and these do not make up a large enough fraction of supernovae to counterbalance their lower energy output <cit.>.
As a result, the cosmic ray production is proportional to the number of stars which can undergo as supernova other than Type Ia, which is just those with masses greater than 8M_⊙ <cit.>.
Putting these pieces together we may write the total cosmic ray production as

    CRR(t) = q E_ CC/PI/M_⊙∫_8M_⊙ SFR(t - τ(M)) ξ(M,T̅(t - τ(M))) dM,

or

    C̅R̅R̅(t) = ∫_8M_⊙S̅F̅R̅(t - τ(M)) ξ(M,T̅(t - τ(M))) dM.


In addition to being produced, cosmic rays are lost to a variety of processes.
Let τ_ CR be the characteristic loss time-scale.
Then

    d ρ_ CR/dt = CRR(t) - ρ_ CR/τ_ CR.

If we treat the Milky Way as being in cosmic ray steady-state and suppose that τ_ CR is fixed across galaxies then

    CRR_ Milky Way = ρ_ CR, Milky Way/τ_ CR,

so

    d ρ̅_ CR/d(τ_ CR t) = C̅R̅R̅(t) - ρ̅_ CR,

which is solved by

    ρ̅_ CR = ∫_0^tC̅R̅R̅(t')e^-τ_ CR (t-t')  dt'.




 §.§ Gas Inflow


Although our toy model treats the mass distribution of cool gas clumps as fixed ( <ref>) this is only an approximation and so it is important to consider when and how it breaks down.
We begin with a Boltzmann-like transport equation for these clumps:

    m∂ n(m)/∂ t = S(m) - ∂ F/∂ m,

where S(m) is a source of clumps of mass m and F is the rate at which clumps of mass m fragment into clumps of mass m-dm.
In writing equation (<ref>) in this form we approximate fragmentation as a differential process in which N clumps of mass m become N + N dm/m clumps of mass m - dm.
This is, however, consistent with the finding of simulations <cit.> and is conceptually similar to the energy flow assumption underpinning the highly-successful Kolmogorov cascade <cit.> so it is likely not wholly inappropriate.
Note that equation (<ref>) is only valid for m ≫ M_J.
Near and below M_J we expect clumps to either cease fragmenting or else form stars.
We therefore take the star formation rate to be of order F(M_J).

The source term in equation (<ref>) is just gas inflow.
Because fresh gas tends to have large-scale structure we see that S(m) has most of its support at large masses m > M_I, where M_I is just the scale of the typical clump in the inflow.
We therefore have an equation with a source at large m and a sink at m > M_J.
In steady-state, we would then have

    F(M_J < m < M_I) = ∫_M_I^∞ S(m') dm'.

If fragmentation is driven by gravitational interactions then above M_J this process has no characteristic mass scale and hence is self-similar.
We therefore expect that

    F(m) ≈1/t_c n(m) (m/M_J)^α,

where α is a dimensionless constant, t_c = (G ρ)^-1/2 is the gravitational free-fall time and ρ is the mass density of the clump.
Because F(m) is a constant in the intermediate regime we see immediately that n(m) ∝ m^-α, which is consistent with observations of a power-law distribution of clumps.

When the system is not in steady-state we might imagine perturbing n(m) → n(m) + dn(m) via a perturbation in S(m).
The resulting perturbation propagates in m so that

    ∂(dn)/∂ t = -1/t_c∂ (dn m^α)/∂ m.

If ρ is approximately constant in m then t_c is as well, so the perturbation propagates with characteristic time t_c.
Hence our assumption in ( <ref>) that the distribution is in steady-state is good only so long as S(m) varies slowly relative to the scale t_c.
For typical densities t_c ≈ 10^6  yr, which is quite short relative to the galaxy's dynamical timescale over which infall rates should vary, so our original assumption is justified.



 §.§ Temperature


The final ingredient in this model is a relationship between the stellar population, the cosmic ray density, the CMB and the temperature.
Equation (<ref>) gives the contribution of the cosmic ray density.
The temperature of the CMB is given by

    T̅_ CMB = 2.7 (1+z)  K,

where z is redshift.
Finally the contribution due to stellar radiation is

    T̅_ rad≈(L̅_⋆/M̅_⋆)^1/4,

where L_⋆ is the total luminosity of stars in the galaxy and M_⋆ is the total stellar mass.
To combine these we let

    T̅ = max(T̅_ CR, T̅_ rad, T̅_ CMB),

where T_ CR is the temperature due to cosmic rays, T_ rad is that due to stellar radiation and
This is a somewhat crude means of combining the two contributions, but it is accurate in regimes in which one effect dominates and provides a reasonable approximation even when this is not the case.



 §.§ Summary


Equations (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>) provide the quantitative specification of the feedback loop depicted in Figure <ref>.
Indeed together they comprise a minimal model of this feedback loop, neglecting spatial, chemical and magnetic degrees of freedom.  We might therefore hope that they could be understood analytically, searching for solutions with physical meaning.  As a system of multiple integrodifferential equations, these equations have multiple solutions and plenty of scope to support complex behaviour.  We consider the most physically interpretable of these in the following sections.



§ EQUILIBRIUM

 

The observed similarities between the growth of galaxies in different environments should motivate a search for equilibrium solutions.  As with the evolution of a main sequence star, a formal equilibrium is not required, but rather an attractor solution which slowly evolves as material is processed into stars would suffice. 

Thus, we will first search for solutions which would be in equilibrium holding gas availability fixed, then consider their broader time evolution on scales long enough that this approximation becomes invalid.  For any equilibrium solution, all of T, ρ̅_ CR, and SFR are constant.

Holding ρ̅_ CR constant requires that

    CRR = ρ̅_ CR.

Inserting equation (<ref>) we find

    ρ_ CR = ∫_8M_⊙ SFRξ(M,T̅) dM,

where we have dropped the time dependence because we require the solution to be an equilibrium.
This further allows the SFR to be pulled out of the integral, so

    ρ̅_ CR = S̅F̅R̅∫_8M_⊙ξ(M,T̅) dM.


When equation (<ref>) is dominated by the CMB or stellar radiation equation (<ref>) is straightforward to solve because ρ̅_ CR does not affect any of the quantities on the right-hand side and so is a free parameter.
In this case equation (<ref>) is trivially satisfied, so the equilibrium condition reduces to

    T̅ =  constant.

In the CMB-dominated case the temperature evolves according to the CMB and the stellar population is irrelevant.
In the stellar radiation-dominated case the equilibrium condition becomes

    T = T̅_ rad≈(L̅_⋆/M̅_⋆)^1/4.

In equilibrium the IMF is constant, so there are ξ(m,T)  SFRτ(m) stars of mass m present.
The contribution of stars of mass m to L_⋆ is therefore proportional to τ(m) L(m), or its lifetime energy output.
The specific energy output is approximately constant with stellar mass[It actually rises with stellar mass, but does so quite slowly and with finite non-zero asymptotic limits at both high and low mass <cit.>.] because τ∝ m^-2.5 while L ∝ m^3.5, so we expect a single equilibrium temperature in this limit, with temperature set by the spatial extent of the galaxy and the typical specific energy output of stars.
This solution is an attractor because the radiation temperature is set directly by the stellar population, and to leading order the stellar population has fixed specific energy output.

Finally, when cosmic rays dominate equation (<ref>) we may expand equation (<ref>) and approximate it by

    T̅_ CR≈1/2ρ̅_CR^1/3.

Hence we find

    T̅^3    = 1/8ρ̅_ CR
       = S̅F̅R̅/8∫_8M_⊙ξ(M,T̅) dM.

As T̅→ 0 the left-hand side vanishes, but the right-hand side  may be shown to scale as T^-0.3.
As T̅→∞ the left-hand side diverges, while the right-hand side falls as T^-4.3.
By the racetrack theorem there must therefore be at least one crossover point between these limits, so there exists an equilibrium.
Indeed there is precisely one equilibrium, because the left-hand side is monotonically increasing while the right-hand side is monotonically decreasing.

The same argument shows that this equilibrium is also an attractor.  An input temperature above the main sequence solution results in less cosmic ray production than is required to maintain that temperature, and similarly a temperature below the main sequence solution will produce a temperature increase.




In the next section we provide numerical evidence for the existence of unique equilibria dominated by different temperature contributions and show that, excluding large fluctuations and crossovers between equilibria, they indeed appear stable.










§ NUMERICAL SOLUTIONS



The arguments in  <ref> demonstrate that if the gas supply is held constant and the only non-negligible contribution to temperature is from cosmic rays, the system will tend towards an equilibrium solution with constant T and SFR.  However, in practice neither of these is true.  

There is certainly long-term evolution in the gas supply, from a combination of processes such as star formation which use up available gas and mergers which increase the supply.  Therefore, the equilibrium towards which the system tends should be evolving.  Moreover, if the gas supply were to change on timescales shorter than the attractor timescale, it would be possible that the system would not tend towards a common evolution at all.  Similarly, the assumption that the gas temperature is dominated by the cosmic ray contribution might hold when a galaxy has a very high SFR, but should not apply to all galaxies or at all times.  These effects cannot be easily modelled with an analytical solution, but rather require solving the system numerically.

The full complexities of the problem are well beyond the simple model in these equations, and therefore although we track several relevant observables, the solutions presented here do not comprise a full simulation.  Rather, we merely examine the dynamics of solutions of our model in the presence of an evolving gas supply for star formation, an optional halo to supply new gas, and a mass-resolved stellar population.
These numerical studies allow us to both confirm the presence of attractor solutions and examine the long-term evolution of our model.  Although written in integrodifferential form,  the equations may also be cast as a set of ordinary differential equations specifying T, ρ_ CR and the number of stars at each mass and time.  This discretised form is more convenient for numerical analysis and so is the one we use here.  



 §.§ Sample Galactic History


As a sample solution, consider the evolution of a galaxy which will eventually attain a stellar mass of 10^11, similar to the Milky Way, from very high redshift until the present.
At early times the temperature will be dominated by the CMB.
Because this has an explicit redshift dependence, its include means the model is no longer a state machine, but now depends directly on cosmic epoch.
As the CMB cools cosmic rays and stellar radiation ought to come to dominate.
The result is an evolution that has three distinct phases (Fig. <ref>). The corresponding evolution of several other parameters is shown in Fig. <ref>.



At the highest redshifts, the temperature is dominated by the CMB <cit.>.  The IMF at these temperatures is very top-heavy, so that few of these stars are present still at low-redshifts.
However, this phase results in a standardisation of galaxies: two halos containing the same baryonic mass but collapsing at different times will both enter this CMB-dominated phase and be attracted towards similar star formation rates and IMFs.  Because relatively few low-mass stars are made, this phase acts to destroy the history of when a galaxy collapsed.  Thus, by the time that the CMB has cooled to the point that it no longer dominates galactic temperatures, the two galaxies have ended up in very similar states despite their different origins.

As the CMB temperature declines, at some point the temperature becomes dominated by cosmic rays from supernovae.  This is the near-equilibrium, cosmic ray attractor solution described in  <ref>.  With a fixed total gas supply, gas availability declines as stars are formed.  This is true even though much of the mass from supernova explosions produces metal-rich regions which can be cooled to form new stars, since the current SFR always exceeds the supernova rate from the relatively rare, high-mass stars formed in the recent past.  So, this solution produces a slow decline in SFR, along with a corresponding slow decrease in temperature and increasingly bottom-heavy IMF.  Most of even the low-mass stars in the galaxy are therefore formed at higher temperatures than at present, so that using a current Milky Way-derived IMF will overestimate the stellar mass even in the local universe.



In practice, there will additionally be gas inflow from mergers, likely initially dominated by hot gas but some of which will eventually become cool gas that can be processed into stars.  

The inflow of cool or cooling gas from minor mergers is poorly constrained by current simulations <cit.>.  Halo merger rates are more robustly predicted (e.g., <cit.>).  However, even if the dark matter infall rate were well known, the baryon content, and in particular the gas fraction and stellar population of infalling material are unknown.  In Fig. <ref>, evolution for several choices of infalling gas rates are illustrated, including a scenario in which infalling baryons are given the same stellar population and gas fraction as the current state of the larger galaxy.

Depending upon the gas infall rate, a family of attractor solutions can be produced.  As a result, it will be possible to at least broadly match observations of the star-forming main sequence with the correct choice of gas infall rate.  It should be noted that this is not necessarily a feature; the ability to always choose a gas infall rate to match observations means that doing so is not properly a prediction of this model.  To use this as a test, it would be necessary to also have not just measurements of overall infall rates, but also a classification of infalling material to determine how much additional gas is available and how infalling stars add to the stellar population.  However, this model does make several other predictions which might be falsifiable, as discussed in  <ref>.

Finally, at low redshift the SFR drops to the point that starlight provides a greater contribution to the temperature than cosmic rays.  Like with CMB domination, under these conditions the SFR will be suppressed below the equilibrium solution produced by cosmic rays alone.  Therefore, this might provide a mechanism by which galaxies with a substantial, old stellar population suppress star formation and remain quiescent.  This is also the solution which best describes the current state of the Milky Way.



 §.§ Convergence and Rounding


To ensure convergence the simulations were run at a variety of time-steps.  For steps of dt = 3× 10^5 yr and below the results were nearly independent of step size, while above 10^6 yr the results were strongly dependent on step size.  This reflects the fact that the lifetimes of the massive stars which contribute the most to the cosmic ray density are of order 10^6 yr and so are under-resolved for step sizes longer than this.

Rounding does not present an issue for our purposes.  The equations have an attractor solution and do not depend sensitively on the exact state.  Therefore, it should be expected that small rounding errors will not build up over time, which is supported by our finding identical long-term numerical behaviour for slightly different initial conditions.  Over the course of ∼ 10^5 steps we incur some rounding error at the relative level of 10^-15 per step. The populations in bins vary by a few orders of magnitude over the course of the simulation, so absent amplification by the time evolution the maximum rounding-induced drift in any bin is of order one part in 10^7.



§ OBSERVATIONAL EVIDENCE AND TESTS



The model developed here accomplishes our primary goal of producing an attractor solution for star formation rates characteristic of the star-forming main sequence.  However, because the SFR and several other details depend upon the gas inflow rate, which is poorly constrained, for the correct choice of parameters this model can be made consistent with a wide range of possible observed star-forming main sequences.  Thus, the qualitative existence of the star-forming main sequence is robust, but its quantitative details do not present falsifiable predictions.  

Further, because this model is intended to be illustrative rather than complete, there are several important effects which have not been included but which will have a non-negligible impact on quantitative details.  However, beyond the existence of a main sequence, the model makes several other predictions which are at least qualitatively robust.  These can be compared with current and near-future observational constraints in an attempt to evaluate whether this is a plausible attractor mechanism for galaxy evolution.



 §.§ Gas (and dust) temperatures


A generic feature of our model is that the gas temperature ought to be a decreasing function of redshift. At early times this is because the CMB dominates the gas temperature, at intermediate times because the cosmic ray density is falling, and at late times because the stellar population ages and becomes more bottom-heavy.

While the gas temperature is difficult to directly observe, the dust temperature has been observed <cit.> and agrees with this qualitative picture.  At any fixed redshift, dust temperatures along the star-forming main sequence are nearly constant.  At fixed mass, galaxies with star formation rates increasing (decreasing) above the main sequence increase (decrease) sharply in dust temperature, providing the feedback described here.  Further, the characteristic star-forming main sequence dust temperature declines from ∼ 35 K to ∼ 25 K from z ∼ 2 to the present, in line with the evolution shown in Fig. <ref>.

It should be noted that these dust temperatures are not necessarily the same as gas temperatures, because the two may not have sufficient time or cross-section to reach equilibrium.  Further, the nature of photometric measurements is such that these are luminosity-averaged dust temperatures, which may not match those in the coolest, star-forming regions of these galaxies.  Nevertheless, if gas temperatures in star-forming regions exhibit similar behaviour to dust temperatures, this would match the qualitative predictions of the cosmic ray feedback model.  Other models with a similar evolution in SFR can also produce similar behaviour, although an empirical fit to the main sequence alone is likely to produce an exponentially increasing SFR instead <cit.>.

It is still perhaps tempting to compare the predictions of the model presented here more quantitatively with dust temperatures.  However, dust temperatures are currently not sufficiently certain to provide such a test.  The inferred dust temperature depends strongly on the chosen value of effective emissivity β_eff (cf. <cit.>).  Values ranging from 1 to 2.5 are common in the literature <cit.>.  As a result, estimates of dust temperature from these various groups can vary by ≳ 10 K for the same observations of the same galaxy due to modelling differences.  

We note that these different models do agree on general trends in dust temperature, which is why the trend in dust temperatures does provide a test of this model.  It should also be noted that because the cosmic ray-driven model does contain free parameters, most notably with respect to gas infall rates, it will be possible to tune the model (Fig. <ref>) to empirically match the gas temperatures produced by any specific choice of β_eff and other parameters.  For example, a closer match with the <cit.> measurements is shown in Fig. <ref>.  However, because this arises from additional tuning, it cannot provide an additional test of the model presented in this work.





 §.§ High-redshift Metallicity and Dust


A further feature of our model is that many massive, short-lived stars ought to be formed at high redshift and then almost immediately die.  This initially CMB-driven, and later cosmic ray-driven, phenomenon results in increased metal and dust production earlier than would otherwise be expected.  Because the total mass of the high-mass end of the IMF grows rapidly with temperature, and because at high redshift these temperatures are predicted to reach ∼ 40-75 K rather than the 20 K in the Milky Way, this can be a substantial effect, increasing dust and metallicity by a factor of up to 5-10.  

Recent observations have indeed found that the first galaxies are surprisingly dust- and metal-rich, to the point that it is difficult to produce dust <cit.> and possibly also metals <cit.> in such quantities without having had more star formation then currently believed.  The feedback described here therefore not only presents a possible solution to this problem but predicts that such a problem should exist.



 §.§ Chemical Abundances


The star formation history of a galaxy is a key factor in its chemical evolution.
As discussed by <cit.>, and illustrated in Fig. <ref> at higher redshifts we predict higher typical stellar masses.
As a result there ought to be fewer low-mass asymptotic giants (AGB stars) at early times.
This should enhance ^16O and other α–process isotopes relative to carbon <cit.>.
There may also be more ^14N relative to ^12C owing to more stars undergoing hot bottom burning <cit.>.
Both effects have been seen in the Milky Way <cit.>.



These effects ought to be strongest in galaxies which underwent a cosmic ray-dominated phase, and so should be stronger in more massive galaxies.
They should also be stronger in galaxies at in starburst galaxies, which may be in a runaway phase, as well as those at higher redshift.
Correspondingly they should be weaker in dwarf galaxies and in the local universe.
Hence while our model is too simple to predict the magnitude of these features, it does robustly predict the sign of their correlation with redshift, SFR, and galaxy mass.



§ DISCUSSION


 
This work develops a simplified model for galactic evolution based around feedback between gas temperature, star formation rate, and the stellar initial mass function.  For most of the history of a typical galaxy, the temperature is dominated by cosmic rays, which are perhaps the most important component currently not well modelled by numerical simulations.  It is therefore intriguing that considering cosmic rays nearly to the exclusion of all else produces, at least qualitatively, many of the features observed by large populations of star-forming galaxies.  

The most notable of these features is the prediction of two key features of the star-forming main sequence.  Cosmic ray feedback drives galaxies with similar gas mass but in different environments to a similar, attractor star formation rate.  Further, although this attractor solution is time-dependent, for proto-galaxies which produce enough metals to form main sequence stars sufficiently early (likely z≳15) that the initial temperature is dominated by the cosmic microwave background, the relevant time will be cosmic epoch rather than age since collapse.  Thus, galaxies do not evolve as a state machine, but rather with a dependence upon cosmic epoch, one of the most puzzling features of the star-forming main sequence.  

Thus, a cosmic ray-dominated model may be able to predict several of the behaviours that have been most difficult to produce in numerical simulations.  Further, the same model can describe a galaxy under very different conditions, as typical galactic history would include, in order, an initialisation phase, a star-forming main sequence phase, and a quiescent phase. However, a full treatment requires combining cosmic rays with the many other effects which are well-modelled by simulations, rather than the toy model presented here which excludes nearly everything else.  

Similarly, galaxies here are treated as monolithic, with one temperature, gas mass, cosmic ray density, and other properties.  This approximation is used both to make the problem tractable and because photometric template fitting typically produces only luminosity-averaged galaxy properties.  However, both observations and simulation demonstrate strongly that this is a poor model.  It is possible that the scatter in the main sequence, as well as conditions driving some galaxies to exhibit the runaway behaviour which is also a possible solution to these equations, is due to this effect.  However, a discussion of these solutions is reserved for a future companion paper.  

In short, it is intriguing that if one considers cosmic ray-driven evolution and excludes all of the other complex baryonic physics which simulations show to be important, the result is evolution which looks remarkably similar to observed features of evolving galaxies which have otherwise often been difficult to match.  However, we know that these other astrophysical effects are important, and we should not expect to produce a fully predictive model without including them.  It is hoped that presenting the mechanism in this format will motivate simulations which include these effects in combination with all of the others, despite the inherent difficulties.
 



 §.§ Impossibly Early Galaxies


Recent high-redshift surveys have revealed a surprisingly large population of very early, very luminous, and therefore presumably massive, galaxies <cit.>.  Indeed, if these galaxies as truly as massive as has been derived from template fitting (out to z∼ 6-8) or UV luminosities (z≳ 8), they are more massive than produced in simulations <cit.> and should reside in halos too massive to have formed in the available time <cit.>.  Several different solutions have been proposed for this problem, including an abundance matching-derived shift in the stellar baryon fraction <cit.> and evolution in dust properties <cit.>.

This problem was a key motivation for our initial consideration of the relationship between temperature and the stellar initial mass function <cit.>, which showed that the CMB contribution to temperature in star-forming regions would decrease mass-to-light ratios, leading to stellar mass overestimation when using a Milky Way-derived IMF.  However, this effect is insufficient to reduce stellar masses to a locally-measured halo mass to stellar mass relation.  

The model described here predicts that the stellar populations in these galaxies formed at even higher temperatures, because galaxies at those masses were cosmic ray-dominated rather than CMB dominated even out to z∼ 20.  Therefore, mass-to-light ratios would be even lower, potentially to the point that stellar masses could be brought back into agreement with theory.  The large parameter space of our toy model prevents revised stellar masses from being a direct prediction of the model, but a sufficiently large correction is plausible.  Spectroscopy with JWST/NIRSpec will likely be sufficient to determine whether this is the correct resolution of the "impossibly early" galaxy problem.



 §.§ Dark Matter Distributions and Interactions


Because massive galaxies become quiescent by the present epoch, most of their stars today are comprised of the low-mass end of the main sequences formed at higher redshifts.  The model developed here predicts that those main sequences formed at higher temperatures than present star-forming regions of the Milky Way.  Thus, again their stellar masses will be overestimated if a Milky Way-derived IMF is assumed.  For a typical evolutionary history, this results in a 20-40% overestimate of stellar mass <cit.>.

Because stellar mass estimate is inherently difficult, for most applications, this correction is negligible.
However, it may shed light on a critical puzzle relating to the nature of dark matter.  The dark matter profile is inferred for local galaxies and clusters by subtracting the baryonic mass, which is inferred from stellar mass, from the total mass inferred from peculiar velocities (cf. <cit.>).  Recent observations find two populations: some galaxies exhibit a “cuspy” dark matter profile (<cit.>; see also <cit.>), in which the central density increases sharply as would be expected from collisionless, cold dark matter <cit.> and others exhibit a “core” with more constant central dark matter density <cit.>, which would require an additional interaction, either a dark matter self-interaction <cit.> or a dark matter-baryon interaction <cit.>.  Therefore, core profiles provide our strongest current evidence for dark matter interactions beyond ΛCDM.

At most radii, the dark matter density is higher than the baryon density.  For cuspy profiles, the same is true even in the centre of the galaxy, so that a 20-40% overestimate in baryon mass is insignificant.  However, for core profiles, baryons become dominant in the central ∼ 1 kpc in galaxies and a larger radius for clusters, in some cases attaining a maximum density 100 times higher than the constant dark matter density.  For such a galaxy, 20% overestimate in the baryon density would mean that the dark matter density would not be constant, but rather would continue to increase, and will have been underestimated by as much as a factor of 20 in the centre.  This would be sufficient to turn a core profile into a cuspy profile, raising the possibility that with corrected stellar masses, no additional dark matter interactions would be required to match galaxy and cluster rotation curves.  

Whether the mass corrections would truly produce a cuspy profile as predicted for collisionless dark matter, or merely an increased central density, will depend upon model parameters.  It will also be sensitive to the additional effects described above which are beyond the scope of this work.  As with the other predictions of the toy model we develop here, these are only robust qualitatively, not quantitatively.  It should further be noted that in dwarf galaxies the stellar mass is always very low and, if the same is true of baryon densities, local dwarf galaxies with core profiles <cit.> will remain cores after this correction.



The simple feedback described here holds the possibility for producing many observed features of galaxy evolution which have heretofore been difficult to explain.
However its present incarnation is highly simplified and comes with a large permissible parameter space, which prevents us from making robust, quantitative predictions. 
Doing so will require fully integrating cosmic rays, in all of their complexity and difficulty, into numerical simulations which will be capable of making such predictions. 
Some progress on this front has already been made <cit.>, and despite the inherent difficulty, the sharp differences between the properties of cosmic ray-dominated feedback and other models, as well as the problems which might potentially be solved, now make including these effects essential in future simulations.



The authors thank Iryna Butsky, Daniel Ceverino, Kristian Finlator, Vasily Kokorev, Adrian Lopez, and Georgios Magdis for helpful discussions.  CLS is supported by ERC grant 648179 "ConTExt".  The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140.  ASJ thanks the UK Marshall commission and the Flatiron Institute of the Simons Foundation for financial support. Support for JL by the Rose Hills Foundation is appreciated. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958 and by the Gordon and Betty Moore Foundation through Grant GBMF7392. 








mnras















