Understanding the radio luminosity function of star-forming galaxies and its cosmological evolution

We explore the redshift evolution of the radio luminosity function (RLF) of star-forming galaxies using galform , a semi-analytic model of galaxy formation (SAMGF) and magnetizer , a dynamo model of the magnetic field evolving in a galaxy through its formation history. Assuming energy equipartition between the magnetic field and cosmic rays, we derive the synchrotron luminosity of each sample galaxy. In a model where the turbulent speed is correlated with the star formation rate, the RLF is in fair agreement with observations in the redshift range 0 ≤ z ≤ 2. At larger redshifts, the structure of galaxies, their interstellar matter and turbulence appear to be rather different from those at z ≲ 2, so that the turbulence and magnetic field models applicable at low redshifts become inadequate. The strong redshift evolution of the RLF at 0 ≤ z ≤ 2 can be attributed to an increased number, at high redshift, of galaxies with large disc volumes and strong magnetic fields. On the other hand, in models where the turbulent speed in galaxies is a constant or assumed to be an explicit function of the redshift, the redshift evolution of the RLF is poorly captured. The evolution of the interstellar turbulence and outflow parameters appear to be major (but not the only) drivers of the RLF changes. We find that both the small-and large-scale magnetic fields contribute significantly to the RLF but the small-scale field dominates at high redshifts. Polarisation observations will therefore be important to distinguish these two components and understand better the evolution of galaxies and their nonthermal constituents.


INTRODUCTION
The galaxy luminosity function (LF) -the comoving number density of galaxies as a function of their magnitude at a given wavelength -is one of the primary observables used to probe the physics of formation and evolution of galaxies.The observed LFs and their redshift evolution have been measured, to various degrees of detail and reliability, for a wide range of wavelengths from the UV to the radio.When assessed in the framework of a theoretical model of galaxy formation, the LFs provide key insights into the evolution of star formation and feedback processes, the nature of and conditions in the ISM (consisting of gas, dust, magnetic fields, ⋆ charles.jose@cusat.ac.in, lchamandy@niser.ac.in, anvar.shukurov@ncl.ac.uk, kandu@iucaa.in,luizfe-lippeSR@gmail.com, c.m.baugh@durham.ac.uk and cosmic rays), and dynamical processes such as galaxy mergers.Two major approaches used to model galaxy formation are: (i) semi-analytic models (SAMGFs) that account for the complex baryonic processes in evolving dark matter (DM) haloes derived from N -body simulations using a combination of analytic approximations and numerical prescriptions, (ii) hydrodynamic or magnetohydrodynamic (MHD) simulations based on fundamental dynamical equations, that include subgrid models to account for unresolved physical processes(often similar to those used in SAMGFs).There are several galaxy formation models (Trayford et al. 2015;Trčka et al. 2022;Croton et al. 2006;Samui et al. 2007;Jose et al. 2013;Lacey et al. 2016) that interpret the ultra-violet (UV) (Wyder et al. 2005;Page et al. 2021), optical (Blanton et al. 2003;Loveday et al. 2012) and infrared (IR) (Dunne et al. 2011;Gruppioni et al. 2013;Marchetti et al. 2016) LFs and shed light on the evolution of obscured star formation rate (SFR), stellar mass (M * ), ISM, dust and feedback processes in galaxy populations.
Complementary to the UV, optical and IR LFs, is the radio LF (RLF) of star-forming galaxies (SFGs) at the emission frequency 1.4 GHz obtained with the Very Large Array (VLA) and the Westerbork Synthesis Radio Telescope (WSRT) (Adams & van Leeuwen 2019;Condon et al. 1998).Recent estimates of RLFs of SFGs at 1.4 GHz include those of Sadler et al. (2002); Condon et al. (2002); Best et al. (2005); Mauch & Sadler (2007); Padovani et al. (2011) and Condon et al. (2019) for the local Universe and from Smolčić et al. (2009); Pracy et al. (2016); Novak et al. (2017); Bonato et al. (2021); Malefahlo et al. (2022); Enia et al. (2022) and van der Vlugt et al. (2022) for redshifts up to z ≃ 4.5.Moreover, several ongoing and upcoming deep surveys using present and next-generation radio telescopes, including the Square Kilometre Array (SKA), will improve our knowledge of the RLFs of SFGs to higher redshifts and fainter luminosities with unprecedented accuracy (Norris et al. 2013;Adams & van Leeuwen 2019).For example, deep surveys using SKA Phase-I will detect sub-µJy SFGs up to z ≃ 6 which will help to measure RLFs down to luminosities more than an order of magnitude lower than what is available now (Jarvis et al. 2015).
Synchrotron emission dominates the radio continuum at 1.4 GHz (Ginzburg & Syrovatskii 1965), so models of the galactic magnetic field and cosmic ray propagation should be the starting points in a prediction of the RLF.MHD simulations of evolving galaxies (e.g., Liu et al. 2022;Pfrommer et al. 2022) cannot provide the statistically significant galaxy samples needed to compute RLFs.Moreover, limited spatial resolution prevents such simulations from probing scales as small as 1-100 pc, which is crucial for capturing dynamo processes at the turbulent and global scales, and hence realistic magnetic fields and, consequently, any realistic distribution of cosmic rays controlled by them (see section 13.14 of Shukurov & Subramanian 2021, for a review).Turbulent magnetic fields, produced by the fluctuation dynamo, supernova shock fronts, and tangling of the large-scale magnetic field, require a spatial resolution of order 1 pc to be realistically reproduced (Gent et al. 2021(Gent et al. , 2023)), whereas the mean-field dynamo, that generates large-scale galactic magnetic fields, relies, apart from overall rotation, on the density stratification at scales of order 100 pc (chapter 11 of Shukurov & Subramanian 2021).MHD simulations at an adequate resolution are only available for local ISM regions of order a few kiloparsecs in size.(e.g.Gressel et al. 2008;Gent et al. 2013b,a;Hollins et al. 2017;Gent et al. 2024), and subgrid dynamo models are not yet available despite an effort to relate the statistical properties of magnetic fields in nearby galaxies to theories ( Van Eck et al. 2015;Chamandy et al. 2016; see Beck et al. 2019, for a review).
Therefore, at present, combining SAMGFs with a model for the evolution of magnetic fields and cosmic rays is arguably the only viable theoretical approach for modelling RLFs.
The earliest effort towards incorporating galactic dynamo theory into a hierarchical model for galaxy formation in cold DM cosmology to produce galactic magnetic fields in a cosmological volume-sized sample of galaxies was by Rodrigues et al. (2015).Rodrigues et al. (2019) (hereafter R19) extended this model into a computational framework called magnetizer (Rodrigues & Chamandy 2020) which simulates mag-netic fields in individual galaxies as they evolve.In particular, these authors coupled the magnetizer with the output of the galform SAMGF of Lacey et al. (2016) and Gonzalez-Perez et al. (2014) to obtain simulated magnetic fields for a large sample of galaxies in the redshift range 0 ≤ z ≤ 6 and compared the predictions of mean-field strengths and pitch angles with observations of local galaxies, finding reasonable agreement.
The present work extends the model of R19, to derive the RLF of SFGs at the rest-frame frequency of 1.4 GHz over cosmic history.We calculate the total radio luminosity due to the synchrotron emission based on our magnetic field and cosmic ray models.Any significant differences of this work from the model of R19 are mentioned in the text.
Key quantities that affect the amplification and sustenance of galactic magnetic fields are the gas density and root-meansquare (rms) turbulent speed, as well as the galactic rotation and stratification (the gas scale height).Some observational studies indicate that the gas velocity dispersion in galaxies (which includes a contribution of the turbulence) is correlated with the SFR surface density (SFRD) (Green et al. 2010;Lehnert et al. 2013;Moiseev et al. 2015) while other authors find a correlation with the global SFR (Genzel et al. 2011;Varidel et al. 2016Varidel et al. , 2020;;Zhou et al. 2017;Übler et al. 2019a).
We investigate how the predicted shape and redshift evolution of the RLF are affected by various phenomenological ISM turbulence models and compare the results with observations.By contrast, past works have probed the radio luminosity of isolated template galaxies and used the radio-FIR correlation to deduce the RLF (e.g.Werhahn et al. 2021;Pfrommer et al. 2022;Vollmer et al. 2022;Schober et al. 2023).
The RLFs obtained from observations contain contributions from both active galactic nuclei (AGN) and SFGs (Mauch & Sadler 2007;van der Vlugt et al. 2022).The physical nature and processes responsible for the radio emission in these objects are rather different and should be considered separately.These contributions can be separated, with varying degrees of confidence, using multi-wavelength observations to estimate the RLF of each.In this work, we model the RLF of star-forming disc galaxies only.
This paper is organized as follows.In Section 2, we discuss our models to compute the RLF of galaxies at 1.4 GHz.We present the predictions for the magnetic field in Section 3, and for the RLF in Section 4 along with a comparison with observational data and a discussion of the roles of various physical processes involved.Our results are put into a wider context in Section 5 and conclusions are formulated in Section 6.Throughout this work, we use cosmological parameters based on seven-year WMAP data (Komatsu et al. 2011), which are also used in the galform model of Lacey et al. (2016).We use the same parameter values in magnetizer as R19 unless otherwise stated.

THE MODEL
The RLFs presented in this work are the products of a threelevel model.In the first level, the output of the galform SAMGF of Lacey et al. (2016) is used to derive the properties of a representative galaxy population of about 7 × 10 5 disc galaxies in a comoving volume of 6 × 10 6 Mpc 3 at each redshift in the redshift range 0 ≤ z ≤ 6 (corresponding to 3.5% of the total volume of the N-body simulation used to provide halo merger histories), with stellar masses between about 10 8 M⊙ and 10 12 M⊙.In the second level, magnetizer is used to couple the galactic dynamo theory and the galaxy formation model to predict the evolution of both turbulent and global magnetic fields as a function of time and galactocentric radius in each galaxy.To avoid excessive and unnecessary complications, all galaxy properties are assumed to be axisymmetric; the dominance of axially symmetric large-scale magnetic fields in nearby galaxies is well established (e.g., Beck et al. 2019).In the final third stage, the synchrotron emission at 1.4 GHz is computed for each galaxy to derive the RLF over a wide redshift range.

Galaxy formation model
The galform version that we use (Lacey et al. 2016, see also Cole et al. 2000;Baugh et al. 2005;Bower et al. 2006) incorporates theoretically and empirically motivated models for complex physical processes (star formation, supernova, and AGN feedback, various dynamical processes, the evolution of the gas, stars, and dust, etc.) that baryons undergo in DM haloes.The assembly histories of haloes are described by halo merger trees extracted from N -body simulations (Guo et al. 2013).The model successfully reproduces the observed K-band, optical, near-infrared, and far-UV luminosity functions of SFGs, far-IR number counts and H i and stellar mass functions over a wide range of redshifts along with the Tully-Fisher, metallicity-luminosity and size-luminosity relations at z = 0. Nearly all radio-bright SFGs are spirals (Sadler et al. 2002;Condon et al. 2002).Therefore, as in R19, we only consider disc galaxies in the output of galform, which are defined as galaxies where the bulge mass is less than half the total stellar mass.
galform outputs several global properties of galaxies at every snapshot of the N -body simulation, like the SFR, stellar mass M * , the cold gas masses of the disc and the bulge, and the disc half-mass radius (r 1/2 ), which is assumed to be the same for stars and gas.Several ISM parameters for magnetizer are derived from these galform outputs in the same way as in R19, as briefly discussed below, but we refer the reader to that paper for further details.
R19 focused on the galactic mean fields and therefore selected galaxies with a large disc that can host the mean-field dynamo in a relatively large volume.However, as we discuss below, random magnetic fields generated independently of the mean-field dynamo action make a significant contribution to the total-intensity RLF, and central parts of galaxies can be very bright in the synchrotron.Therefore, we have changed the sample selection criterion and include all galaxies that have a gas disc, either large or small.The motivation is that a higher local gas density in the disc would lead to star formation, turbulence, and, hence, significant magnetic and cosmic ray energy densities.

Derived galactic quantities
To obtain the RLF, we need to deduce the distribution of the gas density and the magnetic field along the galactocentric distance r and distance Z from the mid-plane in each galaxy, as well as its rotation curve and the gas scale height.
The galaxy rotation curve, V (r), is computed by assuming that galaxies have a thin stellar disc with an exponential surface mass density profile, a bulge with the Hernquist (1990) profile and a DM halo with an adiabatically contracted Navarro-Frenk-White density profile (Navarro et al. 1997).The rotation curve is then used to determine the galaxy angular velocity, Ω(r), and the rotational shear rate, S(r) = r dΩ/dr.
As in R19, we adjust the half-mass radius r 1/2 of the galform galaxies to reproduce the observed stellar mass-r 1/2 relation of Lange et al. (2016) and r d = 2.7r 1/2 is adopted as the maximum gas disc radius for computing quantities like the volume-averaged magnetic field strength.
The surface mass densities of the stars and gas in the disc are both assumed to have an exponential radial profile with the scale length rs = r 1/2 /a, where a ≈ 1.68. 1 As in R19, the disc gas mass is separated into diffuse and 'molecular' phases.The ratio of the surface densities of the diffuse and molecular gas components is computed using the empirical relation from Blitz & Rosolowsky (2004, 2006).The density distributions of diffuse gas, molecular gas and stars perpendicular to the galactic disc are assumed to be exponential with different scale heights.Thus, for example, for the diffuse gas density we have where Z is the distance from the galaxy mid-plane, ρ d0 (r) is the mid-plane density and h d (r) is the density scale height, both functions of the galactocentric radius r, determined assuming that the ISM is in hydrostatic equilibrium with a total pressure that includes the thermal, turbulent, magnetic and cosmic-ray contributions.We use empirical relations for the scale heights of the molecular gas (0.032rs) and stars (0.1rs) (Kregel et al. 2002;Licquia & Newman 2016).

Interstellar turbulence
Galactic magnetic fields are amplified and sustained by a turbulent dynamo.Thus, magnetic field models depend on certain statistical properties of the turbulence, like the correlation length l, correlation time τ , and rms turbulent velocity v.As in R19, we adopt l(r) = min[100 pc, h d (r)] and τ = l/v (for estimates of interstellar turbulence parameters, see Chamandy & Shukurov 2020).Thus far, the magnetizer model described is the same as that of R19.However, R19 assumed that the rms turbulent speed is v = 10 km s −1 , close to the sound speed in the warm gas and independent of the redshift and SFR.However, the gas velocity dispersion (which includes an uncertain contribution from the turbulence) appears to depend on the intensity of star formation (see Fig. 1 and the text below).Given this uncertainty, we considered three alternative prescriptions for determining the turbulent speed.
1 To obtain a, we note that M (r) = 2πΣ 0 r 2 s r/rs 0 xe −x dx for the mass of either stars or gas within the cylindrical radius r, where Σ 0 is the surface density at r = 0.As noted in the text, the scale length rs is assumed to be the same for both stars and total gas.Since ∞ 0 xe −x dx = 1, we have M∞ ≡ limr→∞ M (r) = 2πΣ 0 r 2 s .Now, M (r 1/2 ) = M∞/2 by definition, which leads to a 0 xe −x dx = 1/2 and a follows.Interstellar turbulence has various drivers, including star formation and the associated supernova activity, gravitational instability (Krumholz & Burkhart 2016;Krumholz et al. 2018), cosmological gas accretion, and galaxy mergers (Jiménez et al. 2023;Ginzburg et al. 2022).Some models (e.g., Krumholz et al. 2018;Ginzburg et al. 2022) predict that the turbulent speed is independent of the SFR when SFR < SFR0 with a certain threshold SFR0 and increases with the SFR otherwise.The velocity dispersion data, consistent with more recent observations (Varidel et al. 2020;Yu et al. 2021), are shown in the upper panel of Fig. 1, and we use the form with certain constants v0 and c.As suggested by fig. 3 of Yu et al. (2021), we adopt SFR0 = 3 M⊙ yr −1 and estimate v0 = 25 km s −1 and c = 0.50 by fitting the median of the relation shown in Fig. 1.The fitted dependence is shown as a solid curve in Fig. 1, where the one-dimensional velocity dispersion σ has been converted to the three-dimensional speed v = σ √ 3 assuming the turbulence to be isotropic.Further, v is assumed to be independent of r for simplicity.
The data points in Fig. 1 have a large scatter.To account for this, we fit equation (2) to the 16 th and 84 th percentiles of v as a function of SFR, to obtain v0 = 17 km s −1 and c = 0.45 for the 16 th percentile and v0 = 40 km s −1 and c = 0.55 for the 84 th percentile, also shown in Fig. 1.These measures of the scatter are used to estimate the degree of uncertainty of the predicted luminosity functions.We note that the SFRv relation is assumed to be independent of redshift in this model, which we adopt as the fiducial model.
The turbulent speed of v0 = 25 km s −1 can be understood as the volume average of the rms turbulent speed in the multiphase ISM.As an illustration, if the turbulence is transonic in both warm and hot phases, and their respective fractional volumes and sound speeds are fw = 0.9, cw = 10 km s −1 and f h = 0.1, c h = 100 km s −1 , the volume average follows as (fwc 2 w + f h c 2 h ) 1/2 ≈ 33 km s −1 .We note, however, that the fractional volumes of the ISM phases are likely to vary between galaxies with different SFRs and between various locations within a given galaxy.
As an alternative, we consider a model with a constant turbulent speed v = 25 km s −1 suggested at lower SFRs by Fig. 1, independently of the SFR and redshift (whereas R19 used v = 10 km s −1 ).For the gas number density 0.1 cm −3 of the warm interstellar gas, the magnetic field strength at energy equipartition with the turbulence is then about 4 µG, close to the strength of the random magnetic field observed in nearby spiral galaxies (see Section 2.2.1).The increase in the velocity dispersion at large SFR in Fig. 1 can be interpreted to arise partly from chaotic galactic fountains and winds, which are especially vigorous in galaxies with high SFR.This model is called v25-R19 below.
The velocity dispersion is observed to be correlated with other parameters in addition to the SFR.These include the gas fraction, DM halo mass and the stellar mass (e.g., Krumholz et al. 2018;Ginzburg et al. 2022).This implies that the SFR-v correlation evolves with the redshift, which has been verified by the observations of Übler et al. (2019a) and also found by Jiménez et al. (2023) using the EAGLE hydrodynamic simulations (Schaye et al. 2015).Therefore, we also consider a third model, referred to as v(z)-J23, where v is an explicit function of z, as given by the REF-L100 model of Jiménez et al. (2023) and plotted in the bottom panel of Fig. 1.In this model, v does not depend on the SFR explicitly.

Correction to the SFR from used in galform
In the version of galform described in Lacey et al. (2016), stars form out of cold gas in the disc and spheroid, and the corresponding SFR is proportional to the mass of the cold, molecular gas.The ratio of molecular to atomic gas surface densities is assumed to depend on the pressure in the midplane of the disk, using the empirical prescription of Blitz & Rosolowsky (2006) (see section 3.4 of Lacey et al. (2016) for details).The comoving SFR density (SFRD), the SFR per unit comoving cosmological volume (not to be confused with the star formation density within an individual galaxy), from Lacey et al. (2016) is shown as a function of redshift in the top panel of Fig. 2. Also shown is the best fit to the observed redshift evolution of the SFRD, derived by Hopkins & Beacom (2006) from multi-wavelength photometric and spectroscopic observations across a wide range of wavelengths from X-ray to radio.Both the theoretical and observed SFRDs are computed after assuming the IMF of Kennicutt (1983) (see section 6.3 of Lacey et al. 2016, for details).It is clear from Fig. 2 that, while observations and predictions agree quite well at z = 0, the predicted SFRD lies below the fit to the observations at all redshifts, with their ratio shown in the lower panel.To account for this discrepancy and reconcile the SFR with that observed, we multiply the SFR of galaxies from galform at a given redshift with this ratio, taken at the same redshift.

Magnetic field
Large-scale (mean) B and small-scale (random) b magnetic fields contribute comparably to the galactic synchrotron emission (here and elsewhere, a bar above a variable denotes ensemble or volume average).In nearby galaxies, the ratio of their mean energy densities, b 2 /B 2 , often significantly exceeds unity (Beck 2015a;Beck et al. 2019), but this is subject to a range of assumptions involved in the interpretation of Faraday rotation and synchrotron observations, including the assumption of energy equipartition between cosmic rays and magnetic fields.Beyond semi-quantitative observational and theoretical estimates, the absolute and relative strengths of the two parts of an interstellar magnetic field remain somewhat uncertain.
The results of the mean-field dynamo theory used in the text are introduced in Appendix A; further details and references can be found in Shukurov & Subramanian (2021), among a number of other books and reviews (e.g., Brandenburg & Subramanian 2005).

The random (small-scale) field
Random (turbulent) interstellar magnetic fields are produced by the fluctuation (or small-scale) dynamo, the tangling of the large-scale (or mean) magnetic field by turbulence and compression at random shock fronts.
The fluctuation dynamo amplifies any seed magnetic field on time scales comparable to the turbulent eddy turnover time.Order-of-magnitude estimates suggest that the random field reaches energy equipartition with turbulence after a time of the order of 10 Myr (e.g., Beck et al. 1994, andsections 6.7 and13.3 of Shukurov &Subramanian 2021).Since this time is short compared to the time scales on which galaxies evolve, we assume that at any given time and position within the galaxy the random field has already saturated, and its rms strength is proportional to Beq, the field strength corresponding to equipartition with the local turbulent kinetic energy density, where ρ is the gas density (a function of r and Z as well as of the redshift) and f b is a constant.While f b may depend on various parameters of the interstellar turbulence, there is no reliable model at present, so we set it to be constant, chosen based on the existing evidence.
The fluctuation dynamo produces intense magnetic ropes and filaments that do not fill the volume and the rms field strength depends on both the field strength within the filaments and their fractional volume.Most existing simulations of the fluctuation dynamo suggest f b ≃ 0.5 due to this mechanism (Brandenburg & Subramanian 2005, see also section 13.3 of Shukurov & Subramanian 2021, and references therein) in incompressible turbulence and even lower in supersonic flows (e.g., Federrath et al. 2014).However, the simulations are severely limited to modest values of the kinetic and magnetic Reynolds numbers, while f b can depend on them and is likely to do so.Much stronger, superequipartition magnetic fields with f b > 1 cannot be excluded.For example, this would be the case as long as the fluctuation dynamo produces locally helical magnetic fields, thus diminishing the associated Lorentz force and their back-reaction on the velocity field (see section 4.2 of Belyanin et al. 1994, for a heuristic model).The ISM multi-phase structure further complicates the action of the fluctuation dynamo in galaxies (Gent et al. 2023).
The tangling of the mean (large-scale) magnetic field by the turbulence (an integral part of the mean-field dynamo action) is thought to produce a volume-filling random magnetic field of a strength comparable to that of the mean-field (section 13.3 of Shukurov & Subramanian 2021).Since the mean magnetic field can be stronger than Beq if the mean-field dynamo is strong, e.g., in the central parts of galaxies where the differential rotation is especially strong (section 13.7 of Shukurov & Subramanian 2021, and Appendix A), this may lead to f b > 1.
Compression of magnetic fields at random shock fronts driven by supernovae can produce very strong local random magnetic fields, but their fractional volume is likely to be small; nevertheless, the resulting rms field strength can correspond to f b ≃ 1 (Bykov & Toptygin 1987).
We adopt f b = 1, whereas R19 used 0.5.As the random field is taken to be isotropic, b 2 r = b 2 ϕ = b 2 Z = b 2 rms /3 for the cylindrical field components.

The mean (large-scale) field
The mean-field dynamo, by contrast, can amplify a largescale magnetic field B on length scales of the order of a kiloparsec or larger, on time scales that depend on the galactic rotation period, velocity shear due to the differential rotation and the time scale of the turbulent magnetic diffusion across the gas layer.The relevant parameters can be combined into the dynamo number of equation (A2).The amplification time of the mean field ranges from a fraction of a gigayear in the inner parts of spiral galaxies to a few gigayears in the outer parts.We simulate the evolution of the mean component of the magnetic field in galactic discs by solving numerically the mean induction equation supplemented by a dynamical equation that models the non-linear back-reaction from the Lorentz force on the velocity field that quenches the dynamo to establish a stationary magnetic field.The resulting strength of the large-scale magnetic field is proportional to Beq and also depends on the dynamo number and other parameters (Shukurov & Subramanian 2021).The equations solved and the numerical implementation can be found in R19, and they rely on the thin-disc approximation in the mean-field dynamo equations applicable for h/r ≲ 0.1 (Baryshnikova et al. 1987).The upper left panel of Fig. 3 shows that this ratio is less than 0.4 at the half-mass radius for the vast majority of galaxies, and even smaller for the galaxies in the larger stellar mass range and at a high redshift.To ensure that the solutions of the dynamo equations are computed as fast as required with a large sample of galaxies, we employ well-tested (chapters 11 and 13 of Shukurov & Subramanian 2021, see also Appendix A) approximations applicable to a thin disc (R19 and references therein), including the no-Z approximation, which replaces Z (the distance from the mid-plane) derivatives with divisions by the appropriately scaled gas scale height, allowing the equations to depend on only one spatial dimension (r).The model parameters are the same as in R19, with two exceptions.We set Rκ (Chamandy et al. 2014) to 1.5 (Gopalakrishnan & Subramanian 2023), whereas R19 used Rκ = 1.The strength of the mean magnetic field is approximately proportional to Rκ (see Appendix A).The total pressure (the sum of turbulent, magnetic, and cosmic ray pressures) relative to the turbulent pressure is also estimated to be slightly larger than that adopted by R19.
The initial (seed) magnetic field required to launch the mean-field dynamo is the large-scale component of the random magnetic field produced by the fluctuation dynamo: it does not vanish when the random field is averaged over a finite volume of the galactic disc.This provides an initial kiloparsec-scale random magnetic field of the order of 10 −9 G in strength, stronger than any plausible primordial magnetic field (Ruzmaikin et al. 1988;Poezd et al. 1993;Bhat et al. 2016;Gent et al. 2024; see chapter 13 of Shukurov & Subramanian 2021, for a review).We use the same prescriptions as in R19 for initializing the magnetic field at t = 0 and for preventing the mean field from decaying to very low values when the dynamo is subcritical.
When an evolving galaxy experiences a major merger (where the ratio of the combined stellar and cold gas masses involved exceeds 0.3), the mean magnetic field is reset to the seed field assuming conservatively that such a merger destroys the large-scale field or distorts it to such an extent that it is so different from an eigenfunction of the dynamo equations that most of it decay rapidly because of the turbulent diffusion and the dynamo starts again with a weak seed magnetic field described above.

Field distribution across the disc
The procedures described above specify the variation of the random and mean magnetic fields along the galactocentric distance in each galaxy of the sample, but the field distribution across the gas layer (in Z) also has to be specified.In our fiducial model, the scale heights, hB, of both the random and the mean magnetic field are assumed to be the same as that of the diffuse gas density h d .Thus, for example, the distribution in Z of the mean field is adopted as where hB = h d and BM(r) is the solution of the mean-field dynamo equations produced using magnetizer.
To provide a simple measure of the magnetic field strength, we use its strength B0 corresponding to the volume-averaged total magnetic energy density, for the total magnetic field B (where B ≡ |B|) and similarly B0 for its mean part and b0 for the random part.We also introduce a similar quantity for the equipartition magnetic field of equation ( 4),

The synchrotron emission
The mid-radio (1-10 GHz) continuum emission of an SFG consists of the non-thermal (synchrotron) emission due to relativistic (cosmic ray) electrons and the free-free emission from thermal electrons.The radio emission at 1.4 GHz is dominated by the synchrotron, with the thermal emission contributing about 10 per cent (Tabatabaei et al. 2017), and we neglect the thermal contribution to the RLF (see Schober et al. 2023 for a discussion of the thermal contribution at other frequencies).
We assume that the relativistic electrons have an isotropic energy spectrum with a factor KE, where N (E) dE is the number density of electrons with energy between E and E + dE and s is the spectral index.The synchrotron emissivity (the energy emitted by the unit volume of the source per unit time and unit frequency interval within the unit solid angle) in a homogeneous magnetic field B is then given by (e.g., Ginzburg & Syrovatskii 1964;Rybicki & Lightman 1979) where we use the standard notation for physical constants, ν is the emission frequency (= 1.4 GHz in our case), B ⊥ is the and random (lower left) magnetic field strengths, and the ratio of B 0 to B eq,0 (lower right).The magnetic field parameters presented are defined in equations ( 6) and ( 7).At each redshift, the galaxies are separated into two stellar mass bins defined in the legend.The area under each curve is proportional to the number of spiral galaxies in the given mass bin at the given redshift.
Figure 4.The galaxy number density (in comoving coordinates) as a function of various volume-averaged magnetic field properties defined in equations ( 6) and ( 7) for the mean (solid) and random (dotted) fields (upper panel) and the ratio of B 0 to B eq,0 (lower panel) for various redshifts.
total magnetic field strength in the sky plane perpendicular to the line of sight, and where Γ is the gamma-function (see chapter 3 of Shukurov & Subramanian 2021 for details).This expression is applicable to a partially ordered magnetic field, B = B + b in our case, unless the field is inhomogeneous on very small scales: the synchrotron emission pulse from a relativistic electron is formed over a distance of order rB/γ ≃ 10 9 cm(B/1 µG) −1 (where rB is the electron Larmor radius and γ is the particle Lorentz factor), and magnetic fluctuations at larger scales do not affect the applicability of equation ( 9).We adopt s = 3 for the electron spectrum, close to the observed value (Bisschoff & Potgieter 2014; Bisschoff et al. 2019).For s = 3, the statistically averaged value of ε involves the average of b 2 ⊥ , a quantity predictable theoretically as discussed in Section 2.2.1.
To derive the synchrotron luminosity, the coefficient KE in equation ( 8) has to be determined.An often-adopted assumption regarding the energy density of cosmic rays (dominated by cosmic ray protons) is their energy equipartition with the magnetic field, where κcr is the ratio of the energy densities of the relativistic protons and electrons.For E2 → ∞, this yields KE = (s − 2)B 2 E s−2 1 /(8πκcr) and, for s = 3 and κcr = 100 (e.g., Beck & Krause 2005), we have For E1, we adopt the energy above which the electron spectrum measured in the Solar vicinity has the spectral index s ≈ 3, E1 ≈ 8 GeV (see fig. 10.10 of Shukurov & Subramanian 2021, and references therein).Since cosmic rays have a very high diffusivity (e.g., Longair 1994), we assume that their energy density depends on the spatially averaged magnetic field and use B 2 = B 2 + b 2 rms in equation ( 12).The assumption of the local energy equipartition between cosmic rays and magnetic fields lacks any general and convincing evidence.Both the observed synchrotron intensity fluctuations in galaxies (Stepanov et al. 2014) and simulations of cosmic ray propagation rather suggest a weak anticorrelation between their energy densities (Tharakkal et al. 2023c,b,a).However, there are no readily available, simple alternative approaches to estimate the energy density of cosmic rays in galaxies.For the exploratory study of this paper, we adopt the equipartition assumption but a refinement of our model would include an estimate based on the SFR and a cosmic ray propagation model.
Using the large-scale and small-scale magnetic fields computed for each galaxy, we derive B ⊥ and obtain the synchrotron luminosity assuming that the galaxies in the sample have random orientations with respect to the line of sight, (we use the uniform distribution of cos i, where i is the inclination angle of a galaxy to the line of sight), where Ω is the solid angle.The method for computing the integral in equation ( 13) is described in Appendix B.

MAGNETIC FIELDS IN EVOLVING GALAXIES
The top panel of Fig. 4 shows the number density of galaxies in the fiducial model as a function of the volume-averaged field strengths of the mean and random magnetic fields, B0 and b0 defined in equation ( 6) in each galaxy at various redshifts.At z = 0, galaxies with 1 µG ≲ b0 ≲ 5 µG are more common than galaxies with B0 in the same range, while galaxies with B0 < 0.4 µG or B0 > 10 µG are far more common than those with b0 in that range.With redshift, these limits shift to higher values of magnetic field strength.At lower redshifts, the number density of galaxies as a function of b0 has two distinct maxima.The distribution of B0 is wider than that of b0 because brms is directly related to Beq through equation ( 4), whereas B depends on other galactic parameters as well.The lower panel of Fig. 4 shows the number density of galaxies as a function of B0/Beq,0 and since brms = Beq for f b = 1 (Section 2.2.1) and hence b0 = Beq,0, it is clear that most galaxies in our model have a stronger random magnetic field compared to the mean field.As a larger volume-averaged field strength typically results in a higher synchrotron luminosity, magnetic fields in radio-bright galaxies are predominantly random, especially at high redshift (see also Section 4).
Figure 3 shows the evolution of the probability density of the mean and random magnetic field components in the redshift range 0 < z < 3 with galaxies grouped into two stellar mass bins, 8.5 ≤ log(M * /M⊙) < 10.5 and log(M * /M⊙) ≥ 10.5.The probability density is normalized such that the total probability over all redshifts is unity in each mass bin, separately.This allows us to visualize the relative number of galaxies with a given magnetic field both at each redshift and between different redshifts.The highest mass bin contains about 3% of the total number of galaxies in our sample, and the number of galaxies in a given mass bin changes with z owing to star formation, mergers, etc.The probability distribution of B0 has a single peak, with most galaxies having B0 > 0.1 µG at all redshifts.This is different from the results of R19, their fig.6, where a significant fraction of galaxies have very weak magnetic fields of the order of 10 −4 µG (the strength of the seed field) because the mean-field dynamo is subcritical in some galaxies, particularly at higher redshifts.The mean-field dynamo is significantly stronger in our fiducial model because of the higher turbulent speed v. Stronger turbulent flows promote the mean-field amplification despite the increased turbulent diffusion because the dynamo number in the kinematic regime scales as h 2 d /v 2 (see Appendix A) and h d increases with v faster than linearly.Indeed, under hydrostatic equilibrium, the disc scale height is proportional to the gas pressure, which includes the turbulent contribution ρv 2 , so h d increases with v approximately as v 2 , and therefore, the dynamo number increases with v. Our sample includes massive galaxies with a volume-averaged mean field strength as high as 100 µG, even at high redshifts.Since B0 is a volumeaveraged quantity, the local mean magnetic field strength can be significantly higher.This is consistent with the observational results of Geach et al. (2023) where a density-weighted ordered field as strong as 500 µG is estimated for the dense, dust-rich regions of a galaxy at z = 2.6.
Figures 4 and 3 also show that the mean magnetic field across a sample of galaxies decreases with time because of the depletion of the interstellar gas (although individual galaxies may have monotonically increasing field strength depending on their formation history), which is consistent with the results of R19 and Rodrigues et al. (2015).We also find, as R19 do, that the mean field strength shows a stronger scatter for galaxies with higher stellar masses.However, a new feature of the present model is that galaxies with larger stellar masses have stronger mean magnetic fields, on average, whereas in R19 there was less of a distinction.Galaxies with higher stellar mass generally have higher SFR, as seen in the right-hand panels of Fig. 5, where we show the scatter plot of M * versus SFR, with data points coloured according to the volume-averaged total magnetic field strength.The correlation between the SFR and M * (the so-called galaxy main sequence) is consistent with several observational studies and galaxy formation models (Davé 2008;Davies et al. 2017;Bauer et al. 2011;Whitaker et al. 2014;Sparre et al. 2015;Katsianis et al. 2016;Cowley et al. 2017).In our fiducial model, v increases with SFR.The strength of the saturated mean field increases with v because of the increased dynamo number and the enhancement of Beq of equation ( 4).This explains the increase of B0 with the stellar mass seen in Fig. 3.The increase of the total magnetic field strength B0 with the stellar mass, evident in Fig. 5, is explained by the increase of both B0 and b0 (≡ Beq,0) with M * .
Figure 3 also shows the probability density of the volumeaveraged random magnetic field strength b0.Similarly to B0, galaxies of a higher stellar mass have larger b0.This can be explained by the higher SFR of such galaxies, hence higher v according to equation (2) and larger Beq according to equation (4) (we have confirmed that the correlation of ρ with M * is very weak).The distribution of b0 is bimodal for low-mass galaxies, but not for galaxies of a high mass, for reasons not yet quite clear to us.
The lower panel of Fig. 3 shows the distribution of the ratio of the characteristic field strengths B0 and Beq,0 introduced in equations ( 6) and ( 7).The ratio B0/Beq,0 tells us how the volume-averaged magnetic energy density of a galaxy B 2 0 /8π compares to its averaged turbulent energy density B 2 eq,0 /8π and how close are the dynamos to the saturated state.The random magnetic field saturates at Beq in our model while the mean field can be amplified to significantly larger strengths.
The ratio B0/Beq,0 is sensitive to the distribution of the magnetic field within the disc.The energy density of the turbulence, given by B 2 eq /8π from equation (4), decreases with the galactocentric radius r because of the reduction in the gas density ρ alone as we assume that the turbulent speed v is independent of r.The strength of the mean magnetic field B, while generally decreasing with r, can have a different radial distribution as it depends, in particular, on the angular velocity and shear rate of the galactic rotation.As an illustration, consider the radial dependencies B = B exp(−r/rB), Beq = Beq exp(−r/req) and h d = hd e r/r h with certain length scales rB, req and r h for 0 ≤ r < ∞ (and the variables with the tilde denoting the central values, those at r = 0), so that At all redshifts, the probability distributions of B0/Beq,0 have extended tails, especially in the case of massive galaxies, where this ratio exceeds unity (up to three) suggesting that most galaxies have B/ Beq ≳ 1 and/or rB ≳ req.The mean magnetic field strength likely exceeds Beq in the central parts of galaxies where the dynamo is the strongest.A very slow decrease in the strength of the mean magnetic field with r, rB ≫ req, has been observed in nearby galaxies (section 4.2 of Beck 2015b, and references therein).A large radial extent of magnetised region in spiral galaxies is explained by the slow decrease of the dynamo number with r because of the disc flaring which compensates for the relatively slow decrease of the rotation rate in a flat rotation curve, as can be seen from equation (A3).The radial extent of the magnetised region is controlled by the decrease of the gas density and Beq with r.

THE RADIO LUMINOSITY FUNCTION
Figure 6 shows the variation with the redshift of the RLF Φ(L1.4) of SFGs at the rest-frame frequency 1.4 GHz, where Φ(L1.4) is the number of galaxies per unit comoving volume per decade in luminosity, at luminosity L1.4.For comparison, we note that the Milky Way luminosity at 1.4 GHz is L0 = 3 × 10 21 W Hz −1 (Berkhuijsen 1984).The predictions are computed at the redshifts (shown on each panel) closest to the median of the redshift bins used for RLF measurements by Condon et al. (2002) and van der Vlugt et al. (2022) where simulation snapshots are available.We show the RLFs obtained from two models described in Section 2.1.2:the fiducial model (where the turbulent speed v depends on the SFR) and the model v25-R19 where v = 25 km s −1 in all galaxies.The observed luminosity functions are from Condon et al. (2002Condon et al. ( , 2019) ) and Mauch & Sadler (2007) for the local Universe and Novak et al. (2017) and van der Vlugt et al. ( 2022) for the higher redshifts.
The agreement of the fiducial model with the data is quite satisfactory, if not perfect, over a wide redshift range.In particular, for z ≲ 1.2, the high-luminosity end of the RLF is reproduced remarkably well.Although the median luminosity is under-predicted at higher redshifts for all values of L1.4, the model and observations still marginally agree for z ≲ 2, including the changes in the slope around L1.4 = 10 23 W Hz −1 for z ≲ 0.5.Moreover, when extrapolated to higher luminosities, the median outcome of the fiducial model is in agreement with the data point at the highest luminosity at small z (our galaxy sample is not large enough to probe such bright, rare galaxies at certain redshifts, including z ≈ 0).However, the number density of galaxies in the luminosity range 20.5 ≤ log(L1.4/W Hz −1 ) ≤ 21.5 in the local Universe is significantly higher than what is observed.This local maximum at lower luminosities is sensitive to the poorly constrained ratio of the mean and random magnetic fields (Section 4.2) and might be adjusted by the fine-tuning of our model, which ) of 10 3 randomly selected galaxies with L 1.4 ≥ 10 21 W Hz −1 at z = 0.0, 0.5, 1.0 and 1.5 versus the fourth power of the volume-averaged magnetic field strength and the volume of the emission region (disc volume).Right: the colour-coded volume-averaged magnetic field strength of galaxies shown in the left-hand panels versus their star formation rate and stellar mass.
we avoid.Also, for z > 2, the model number densities are systematically and significantly smaller than those observed, and the discrepancy between the predictions and the data increases with redshift.

The role of the turbulent speed
The turbulent speed v is among the model parameters that affect directly the RLF because the strength of the magnetic fields on both large and small scales increases in proportion to the turbulent kinetic energy via their dependence on Beq through equation ( 3), and the large-scale field strength has additional dependence on v via the α effect and turbulent diffusion.The models v25-R19 and v(z)-J23 are designed to explore the role of this parameter.
Apart from the fiducial model results, Fig. 6 also shows the prediction of the v25-R19 model (dotted black) which has the constant turbulent speed v = 25 km s −1 in all galaxies, in contrast to the fiducial model where v is larger in galaxies with SFR ≥ 3 M⊙ yr −1 (Section 2.1.2).The two models differ little in the predicted RLFs at low luminosities, especially at low redshifts where SFR < 3 M⊙ yr −1 and hence v = 25 km s −1 in most galaxies.However, the increase in the SFR up to z ≃ 2 makes the v25-R19 model predictions differ from those of the fiducial model, and they become less and less adequate with increasing z, severely under-predicting the number density of galaxies at the high luminosity end.Even if the overall luminosity were increased by a factor by adjusting model parameters (as discussed in Section 5), the shape of the v25-R19 RLF does not match the data as well as the fiducial model.overall shape nor the magnitude of the resulting RLFs agree with the observations: this model is inferior to the fiducial one.
In fact, neither the fiducial nor the v(z)-J23 model is fully justified because model v(z)-J23 does not capture the SFRdependence of turbulence at constant z, while the fiducial model does not allow for such effects as the gravitational instability (Krumholz et al. 2018) and gas accretion and mergers (Ginzburg et al. 2022) which are likely to result in an additional dependence of the turbulent intensity on redshift independently of the SFR.A more realistic model for turbulence will be among the refinements of our model to follow.

Relative importance of the random and mean magnetic fields
Under the assumption of energy equipartition between cosmic rays and magnetic field adopted here, the contributions of the mean and random magnetic fields to the RLF are not additive because the radio luminosity includes the integrals of Therefore, to assess the significance of the mean and random magnetic fields, we consider results obtained when one of them is ignored.In Fig. 7, we show the predictions with only the random field (dash-dotted cyan line) or the mean field (dashed blue line) retained.Both random and mean fields are clearly important.However, the RLF that relies solely on the mean magnetic field decreases with the luminosity too rapidly and by z ≈ 1 the contribution from the random field dominates that of the mean-field at least at those high luminosities where observational data are available.On the other hand, at z ≈ 0, we obtain a significantly better fit to the data by neglecting the random field and including only the mean field, particularly at low luminosities where the fiducial model over-predicts the RLF.
We stress that the relation between the contributions of the random and mean magnetic field components to the RLF depends on poorly constrained parameters like fB for the random field and the magnitude of the mean helicity of interstellar turbulence for the mean field.Since the radio luminosity scales as the fourth power of the magnetic field (Section 4), the relative influence of the random and mean magnetic fields is rather sensitive to such parameters.Clearly, the same magnitude of the luminosity function could be obtained by changing the parameter values so as to boost one field component and reduce the other.Therefore, the relative importance of mean and random field is not well-constrained by the model and for realistic parameter values, neither can be neglected.
The current interpretations of the observations of the synchrotron emission in nearby galaxies suggest that the energy density of the random magnetic field exceeds that of the large-scale field by a factor of three or even larger (Beck et al. 2019).However, the interpretations rely on the assumption of the local energy equipartition between cosmic rays and magnetic fields, and so are model-dependent.
To conclude, there are several plausible reasons for the apparent complexity in the relative roles of the mean and random magnetic fields at different redshifts.It is likely that a better model for cosmic rays is necessary (e.g. because their number density might be more sensitive to the mean magnetic field than to its random part).Another reason could be our limited understanding of the interaction of the fluctuation and mean-field dynamos (chapter 8 of Shukurov & Subramanian 2021).Both types of dynamo action are sensitive to the multi-phase ISM structure but in a poorly understood manner while it is clear that the ISM structure depends on both the SFR and redshift.Remarkably, the RLF can shed light on rather subtle aspects of both the ISM structure and galactic dynamos.The relative importance of mean and random fields thus deserves careful analysis beyond the scope of this paper.

The redshift evolution of the RLF
The redshift evolution of the observed RLF for bright starforming galaxies (L1.4 ≳ 10 21 W Hz −1 ) is discussed by Novak et al. (2017) and van der Vlugt et al. (2022van der Vlugt et al. ( , 2023)); their observational data are shown in Figs 6 and 7 together with our results.Both the fiducial model (where the turbulent speed increases with the SFR) and the observations agree that the range of the radio luminosities extends to larger L1.4 as z increases.The galaxy luminosity functions observed in far infrared (Gruppioni et al. 2013;Koprowski et al. 2017;Lim et al. 2020) and K-band (Cirasuolo et al. 2010;Mortlock et al. 2017) evolve similarly with the redshift.However, the RLF shows a slightly stronger evolution compared to that in the infrared as the infrared-to-radio luminosity ratio is observed to decrease mildly with increasing redshift (Delhaize et al. 2017;Calistro Rivera et al. 2017;Magnelli et al. 2015;Ivison et al. 2010).
Both the size of the emission region and the strength of the magnetic field contribute to the increase of the synchrotron luminosity with z.This is illustrated in the left-hand panels of Fig. 5, which presents the colour-gradient scatter plot of the radio luminosity of a thousand randomly selected galaxies with L1.4 ≥ 10 21 WHz −1 from our sample at different redshifts as a function of their disc volume and B 4 0 (as L1.4 ∝ B 4 under the assumption of energy equipartition between cosmic rays and magnetic fields).It is clear that the brightest galaxies typically have both large emitting volume and volumeaveraged magnetic field strength.For example, the volumeaveraged magnetic field strength of the brightest galaxies at z = 1.5 is about a factor of two larger than at z = 0 Certain factors that contribute to stronger magnetic fields at high redshifts are clarified in the right-hand panels of Fig. 5. Galaxies with larger M * tend to contain stronger magnetic fields, for reasons explained in Section 3. The SFR density has a maximum at z ≃ 1.0 (Fig. 2), indicating high SFR in galaxies at high redshifts.This results in a higher turbulent speed in galaxies via equation 2 leading to a systematic increase of the magnetic field strength with z at 0 ≲ z ≲ 1.5.
Thus, our fiducial model finds a shift to higher luminosities with increasing z, in agreement with the observations up to z ≈ 1.8.However, as mentioned at the beginning of Section 4, the z-dependent shift in the model RLF is weaker than that inferred from observations.The median predictions of the model seem to underestimate the observed RLF at z ≳ 1, and the difference between them increases with increasing z, and at z ≳ 2 they clearly disagree for reasons discussed in the next section.

DISCUSSION
Our goal in this paper is, firstly, to explore the applicability of the theory and models of galactic magnetic fields to young and evolving galaxies and, secondly, to identify the most important galactic properties and parameters that can be deduced and understood with the help of their radio luminosity function.The synthetic RLFs of star-forming galaxies obtained from their modelled synchrotron emission at 1.4 GHz are in broad agreement with the available observations, reproducing the number density of galaxies in the luminosity range 10 19 ≲ L1.4 ≲ 10 25 W Hz −1 (3 × 10 −3 to 3 × 10 3 in terms of the Milky Way luminosity at 1.4 GHz -Berkhuijsen 1984) in the redshift range 0 ≲ z ≲ 2, within the uncertainties in the turbulent speed discussed in Section 2.1.2.Moreover, the form (e.g., changes in the slope) of the RLF agrees well with what is observed.We do not consider z ≳ 2, because the synchrotron luminosity predicted by the model becomes significantly lower than what is observed, continuing the trend with z already visible in Fig. 6.
We have made no persistent attempt to achieve any better agreement with the observational results even though this would be possible given that many properties of the ISM in evolving galaxies leave much freedom for adjustment.We avoid making heuristic assumptions that do not have explicit observational or theoretical justification.
We have identified the turbulent speed v as one of the predominant factors affecting the RLF.Our results firmly indicate that v must be an increasing function of the star formation rate for the SFR exceeding a certain threshold (here adopted as SFR0 = 3 M⊙ yr −1 ), and the model with v = const (v25-R19) cannot explain the increases in the galactic luminosity with z even though its results are close to those with the SFR-dependent turbulent speed (the fiducial model) at z ≈ 0 (Fig. 6).
We also find that the shape and redshift-dependence of the RLF cannot be reproduced by making v depend explicitly on z instead of on the SFR.This is elucidated by the model v(z)-J23, which predicts a decrease of the RLF with luminosity that is too steep, failing badly to reproduce observations.We appreciate, however, that the variation of the turbulent speed with the SFR adopted in the fiducial model may just serve as a proxy for other effects discussed below.

Interstellar magnetic fields
The energy density of the random magnetic field is proportional to the energy density of the turbulence with the proportionality factor f 2 b of equation ( 3) only known to be of order unity.Adopting f b = 1, we find a reasonable agreement between the predicted RLF and observations at the redshifts z ≲ 1.2 (see Fig. 6).Our knowledge of the strength and structure of turbulent magnetic fields largely relies on numerical simulations severely limited to modest Reynolds numbers.For example, adopting f b = 1.5 improves significantly the agreement at higher redshifts.Values of f b exceeding unity can be due to a stronger magnetic helicity since partially helical fields are close to the force-free state, so their back-reaction on the flow is reduced, which allows them to have a higher energy density.Thus, our model might fit the observations much better, had we assumed that f b increases slightly with z.Moreover, f b can well depend on various galactic parameters and be sensitive to the multi-phase ISM structure rather than be a constant or a simple function of z.Our results highlight the need to understand better how random magnetic fields are maintained in galaxies.The theory of the mean magnetic fields is developed better although many aspects of the action of the mean-field dynamo in the multi-phase, evolving ISM remain speculative.Since the energy density of both the turbulent and large-scale magnetic fields depends on the turbulent energy density, the gas density and turbulent speed are identified as major factors affecting the RLF.However, the mean magnetic field has a more complicated dependence on the galactic parameters than its random part, including the differential rotation and the gaseous disc thickness in addition to the turbulent energy density.The energy density of the mean magnetic field (unlike the random field in our model) can exceed the turbulent energy density.As a result, the total magnetic energy density in a significant fraction of galaxies (especially massive ones) exceeds the turbulent energy density (see also fig.6 of Rodrigues et al. 2019).The RLF in polarised emission would help to isolate the contribution of the mean magnetic field to the synchrotron luminosity and thus assess rather subtle properties of evolving galaxies.Our predictions for the polarised RLF will be published elsewhere.

Galactic structure
The predictions of our fiducial model, which assumes that the structure of a galaxy is broadly similar to a spiral galaxy in the local Universe, are reasonable until z ≈ 2, the redshift at which galaxies appear to develop persistent thin discs dominated by rotation similar to those at z ≈ 0 (Ginzburg et al. 2022;Jiménez et al. 2023, and references therein).Our models implicitly assume that the disc is formed instantaneously whereas the history of the disc formation can affect the outcome of the evolution in the strongly nonlinear galactic systems.Certain effects of enhanced star formation on the galactic structure and environment, such as outflows (fountains and winds) and the emergence of radio haloes, are not captured by our model which only includes magnetic fields in galactic discs.Such effects can affect high-redshift galaxies more profoundly and may help to explain why our model does worse at matching observations as z increases.For example, galaxies with a high SFR may have radio haloes with vertical extent 5-10 kpc and magnetic field strength comparable to that in the disc.Thus, future improvements to the magnetic field model should include the radio halo and its contribution to the synchrotron luminosity.

Cosmic rays
To obtain the number density of cosmic ray electrons, we have adopted a widely used (albeit poorly justified) assumption of energy equipartition between cosmic rays and magnetic fields.However, the energy density of cosmic rays (including the number density of relativistic electrons) is likely to be an increasing function of the supernova rate directly related to the SFR.We shall explore elsewhere an alternative model where the abundance of cosmic rays is controlled by their sources and propagation in the evolving magnetic field.
The predicted RLF is also sensitive to the slope s and amplitude KE of the energy spectrum of cosmic ray electrons given by equation ( 8).Under the assumption of energy equipartition between cosmic rays and magnetic field, changes to the value of s also affect KE (in proportion to s − 2) and hence the predicted luminosity of all galaxies, shifting the RLF along the horizontal (log L) axis.The shape of the RLF also can be mildly modified as the synchrotron emissivity is proportional to B (s+1)/2 ⊥ (see equation 9).Since the emissivity changes as ν −(s−1)/2 , multi-wavelength radio observations can constrain s and resolve the degeneracy between KE and s.
The estimate of the synchrotron emissivity of equation (9) needs to be refined at high redshifts because the energy losses of relativistic electrons to the inverse Compton scattering off the CMB photons increase as (1 + z) 4 .The ratio of the rates at which a relativistic electron loses energy to the inverse Compton scattering and synchrotron is given by 8π wCMB/B 2 , where wCMB is the CMB energy density (section 3.1.4of Shukurov & Subramanian 2021).For electrons with the Lorentz factor γ, the frequency of the CMB photons is boosted by the factor γ 2 , so 1 GeV electrons emit in the X-ray range.Schleicher & Beck (2013a) suggest that the inverse Compton losses dominate over the synchrotron emission at z ≳ 2 if the typical strength of the galactic magnetic fields B does not increase with z but at a larger z if B increases with the redshift.The effect of the inverse Compton losses on our results at z ≤ 2 is likely to be only modest, especially for the radio-bright galaxies.

Interstellar turbulence and galactic outflows
It appears that the mechanisms of magnetic field generation, and the ISM structure on which they rely, experience a significant change at z ≃ 2. This is suggested, in particular, by the increase in the gas velocity dispersion with the SFR (Fig. 1) above the sound speeds in the diffuse ISM phases, cs ≃ 10 km s −1 in the warm gas and cs ≃ 130 km s −1 in the hot phase.Meanwhile, simulations of the supernovadriven multi-phase ISM show that the fractional volume of the hot gas at the galactic mid-plane increases only slightly (as SFR 0.36 assuming that the supernova rate is proportional to the SFR) from about 0.2 to 0.28 when the supernova rate increases by a factor of four in comparison with the Milky Way rate and does not exceed 0.5 when the supernova rate is 16 times the local rate, while the abundance of the cooler phases reduces as the supernova rate increases (de Avillez & Breitschwerdt 2004; Breitschwerdt & de Avillez 2021, and references therein).The main consequence of the increase in the SFR is a more vigorous outflow of the hot gas in the form of the galactic fountain or wind.Although the off-planar gas is hot and highly ionised, significant amounts of neutral hydrogen are present, entrained from the disc, or cooled down in situ).Its H i velocity dispersion is about 60 km s −1 in the Milky Way (table 1 of Kalberla 2003) while the velocity dispersion of the hot ionised gas is 75 km s −1 (section 8.3 of Savage et al. 1997).A review of relevant observations and further references can be found in section 10.2.2 of Shukurov & Subramanian (2021).Meanwhile, the turbulent velocity in the disc hardly changes as the supernova rate varies by a factor 512, remaining close to the sound speed in the warm gas while the H i line width is larger than the (one-dimensional) turbulent velocity dispersion (Joung et al. 2009, see also Dib et al. 2006).Therefore, an increase of the galactic outflow intensity and an enhanced contribution of turbulence in the off-planar gas appear to be the dominant effects of the increased SFR on a spiral galaxy.In this context, the increase of the turbulent speed with the SFR adopted in the fiducial model can be interpreted as a reflection of the increasing importance of the off-planar gas with its more intense turbulence.
Similarly to R19, we assume that the rms turbulent speed is independent of the galactocentric distance r.Some observations of the H i velocity dispersion in spiral galaxies show that it decreases with the galactocentric radius r from often supersonic values in the inner parts of galaxies (Boulanger & Viallefond 1992;Petric & Rupen 2007;Tamburro et al. 2009;Mogotsi et al. 2016).On the other hand, there is significant observational evidence for a very weak variation of the H i velocity dispersion with r in the discs of spiral galaxies (Dickey et al. 1990;Blitz & Spergel 1991, Section 12.2.3 of Kamphuis 1993).It is plausible that the turbulent speed is indeed independent of the galactocentric radius and the spatial variation of the gas velocity dispersion is due to variations in the outflow intensity.This aspect of the model requires further analysis.
The models and observational estimates of the H i and Hα velocity dispersions discussed in Section 2.1.2are often interpreted as suggesting strongly supersonic turbulence, especially at high SFR.However, the turbulence in the diffuse (warm and hot) interstellar gas is likely to be transonic (or subsonic, depending on the kinetic energy injection rate) because of the self-regulation of the supersonic turbulence (Section 2.10.2 of Shukurov & Subramanian 2021).The kinetic energy of supersonic turbulence efficiently dissipates at shock fronts to heat the gas until an equilibrium state is reached in which the turbulent speed is comparable to the sound speed.Thus, turbulent flows with sufficiently strong energy injection rates are likely to be transonic.Such a self-regulation of a turbulent system can be affected by radiative cooling which can prevent the gas heating by removing the dissipated turbulent energy via radiation (Enrique Vázquez-Semadeni, private communication).The cooling is especially efficient in dense regions (like molecular clouds which occupy a negligible fraction of the disc volume), where supersonic turbulence can be maintained, but perhaps not in the warm and hot ISM phases which dominate the radio luminosity.The self-regulation of supersonic turbulence is not included in the turbulence models of Krumholz et al. (2018) and Ginzburg et al. (2022), who assume that the disc is marginally stable concerning the gravitational instability at all times and consider only the energy balance of turbulence driven by supernovae and accretion flows onto and through the disc.The sample of SFGs where Varidel et al. (2020) measured the vertical velocity dispersion in the Hα spectral line deliberately includes only nearly face-on galaxies (0 < i < 60 • for the inclination angle), so not only the turbulence in the galactic discs but also outflows and turbulence of the off-planar gas are likely to contribute to the velocity dispersion obtained.
The dependence of v on the SFR, equation (2), is rather poorly constrained by the data available, and other fits can be equally acceptable.For example, adopting SFR0 = 1 M⊙ yr −1 and c = 1/3 in our fiducial model produces an RLF that is in reasonable agreement with observations, indicating a degeneracy between these parameters (lower values of the scaling exponent c would require lower values of SFR0).Furthermore, there are alternative prescriptions to derive the velocity dispersion from the star formation surface density (Niklas & Beck 1997;Chyży et al. 2011;Schleicher & Beck 2013b;Chamandy et al. 2024).Finally, the velocity dispersion may increase faster with redshift than expected from the SFRredshift dependence alone (Wisnioski et al. 2015;Übler et al. 2019b), so a better understanding of the interstellar turbulence and the role of the evolving multi-phase ISM structure appears to be essential for further progress in modelling the RLF, which we will address in the future.

CONCLUSIONS
The magnetic field and dynamo models presented here lead to a satisfactory agreement, within uncertainties, of the predicted RLF of star-forming galaxies with observations up to z ≃ 2. At higher redshifts, although the form of the theoretical luminosity function is similar to what is observed, the theory predicts a smaller number of galaxies of high luminos-ity, and the discrepancy increases with the redshift.While the observational data at higher redshifts still lack galaxies of faint and moderate luminosity, our findings enable us to identify a range of effects that need to be understood and included to improve the agreement with the data.
We have identified the strength of interstellar turbulence as one of the major factors affecting the RLF.A better understanding of the interstellar turbulence, primarily the turbulent speed (which makes only a limited and uncertain contribution to the observed spectral line widths) at high redshifts and the effects of the evolving multi-phase structure of the interstellar medium on magnetic fields and cosmic rays are required to improve our models.
Models presented here are based on a well-tested galaxy formation model; however, we rely on an implicit assumption that the overall structure of the interstellar medium does not change as the galaxies evolve.This is an oversimplification but little is known about some aspects of galactic evolution which are relevant to magnetic fields and cosmic rays (the turbulent parameters, gas disc thickness, the multi-phase ISM structure, etc.).It is then not surprising that our model becomes inapplicable beyond the redshift z ≃ 2 at which galaxies undergo a significant structural change, the development of pronounced gas discs in particular.Among distinct structural features of young galaxies, important for their nonthermal constituents, might be widespread radio haloes associated with stronger star formation and correspondingly vigorous galactic fountain flows.Our limited knowledge of the properties of turbulent magnetic fields, especially in the evolving multi-phase ISM, is another question that requires a better answer.
We have shown that the galactic luminosity function in the total radio intensity is a sensitive indicator of the state of the interstellar gas, especially its turbulence, density, and galactic fountains and winds.Our results are consistent with the notion that the structure of star-forming galaxies, and their magnetic fields, are different at z ≲ 2 and the higher redshifts.The luminosity function in polarised radio emission would provide even richer information on galactic rotation, gas stratification, and anisotropy of interstellar turbulence.

Figure 1 .
Figure 1.Top panel: the observed three-dimensional velocity dispersion as a function of SFR (Krumholz et al. 2018 and references therein, circles) and the best-fitting curves of the form of equation (2) to the median (solid/green) and the 16th (dash-dot/red) and 84th (dashed black) percentiles.The dotted blue line corresponds to v = 25 km s −1 .Bottom panel: The turbulent speed of the cold interstellar gas (T ≤ 10 K) as a function of redshift from the EA-GLE cosmological simulations for z ≤ 4 (solid line: Jiménez et al. 2023), extrapolated to z > 4 as a constant (dotted).

Figure 2 .
Figure 2. Top panel: the best-fitting dependence of the observed SFRD on the redshift from Hopkins & Beacom (2006) (dashed/blue) and the SFRD from Lacey et al. (2016) (solid/red).Bottom panel: the ratio of the SFRDs shown in the top panel.

Figure 3 .
Figure3.The probability distributions of the relative gas disc thickness h d (r 1/2 )/r 1/2 (upper left), volume-averaged mean (upper right) and random (lower left) magnetic field strengths, and the ratio of B 0 to B eq,0 (lower right).The magnetic field parameters presented are defined in equations (6) and (7).At each redshift, the galaxies are separated into two stellar mass bins defined in the legend.The area under each curve is proportional to the number of spiral galaxies in the given mass bin at the given redshift.

Figure 5 .
Figure5.Left: the scatter plots of the colour-coded radio luminosity log(L 1.4 ) of 10 3 randomly selected galaxies with L 1.4 ≥ 10 21 W Hz −1 at z = 0.0, 0.5, 1.0 and 1.5 versus the fourth power of the volume-averaged magnetic field strength and the volume of the emission region (disc volume).Right: the colour-coded volume-averaged magnetic field strength of galaxies shown in the left-hand panels versus their star formation rate and stellar mass.

Figure 6 .
Figure 6.The 1.4-GHz radio luminosity function of star-forming galaxies at various redshifts z specified within each frame.The median outcome of the fiducial model is shown as a red solid line with the shaded region representing the 16-84 percentile range of the turbulent speed shown in Fig. 1.The v25-R19 model (with the constant turbulent speed v = 25 km s −1 ) is represented with the dotted black line.The observational data are taken from Condon et al. (2002) (black crosses), Condon et al. (2019) (blue crosses), Mauch & Sadler (2007) (red circles), Novak et al. (2017) (black circles) and van der Vlugt et al. (2022) (green triangles).

Figure 7
Figure 7 compares the fiducial model with the model v(z)-J23 (dotted black), where v is an explicit function of the redshift taken from the bottom panel of Fig. 1.Neither the

Figure 7 .
Figure 7.As Fig. 6, but the fiducial model compared with the v(z)-J23 model (black dotted) where the turbulent speed is an explicit function of the redshift.The dashed blue and dash-dotted cyan lines show the contributions to the synchrotron luminosity due solely to the mean and random magnetic fields, respectively.