The LOFAR Two-metre Sky Survey: the radio view of the cosmic star formation history

We present a detailed study of the cosmic star formation history over $90$ per cent of cosmic time ($0\lesssim z\lesssim4$), using deep, radio continuum observations that probe star formation activity independent of dust. The Low Frequency Array Two Metre Sky Survey has imaged three well-studied extragalactic fields, Elais-N1, Bo\"otes and the Lockman Hole, reaching $\sim20\,\mu\rm{Jy/beam}$ rms sensitivity at $150\,\rm{MHz}$. The availability of high-quality ancillary data from ultraviolet to far-infrared wavelengths has enabled accurate photometric redshifts and the robust separation of radio-bright AGN from their star-forming counterparts. We capitalise on this unique combination of deep, wide fields and robustly-selected star-forming galaxies to construct radio luminosity functions and derive the cosmic star formation rate density. We carefully constrain and correct for scatter in the $L_{150\,\rm{MHz}}-\rm{SFR}$ relation, which we find to be $\sim0.3\,\rm{dex}$. Our derived star formation rate density lies between previous measurements at all redshifts studied. We derive higher star formation rate densities between $z\sim0$ and $z\sim3$ than are typically inferred from short wavelength emission; at earlier times, this discrepancy is reduced. Our measurements are generally in good agreement with far-infrared and radio-based studies, with small offsets resulting from differing star formation rate calibrations.


INTRODUCTION
Characterising the history of cosmic star formation and understanding the drivers of star formation in galaxies from high to low redshift have been fundamental goals of extragalactic research for several decades (see Madau & Dickinson 2014, for a thorough review).Early high-redshift surveys presented a broad view of star formation over cosmic time, with the volume-averaged star formation rate density (SFRD) increasing from high redshift to peak somewhere in the range 1 ≲  ≲ 2.5 and then declining towards the present day (Lilly et al. 1995(Lilly et al. , 1996;;Madau et al. 1996;Connolly et al. 1997;Pascarelle et al. 1998).By the mid-2000s, the SFRD had been constrained fairly tightly back to  ∼ 1, using a range of tracers (Hopkins & Beacom 2006;Wilkins et al. 2008).However, its form at higher redshift and the exact position of the peak remained less well determined.
Since then, numerous studies have attempted to constrain the SFRD more tightly, particularly at higher redshifts, and ultraviolet ★ E-mail: rcochrane@flatironinstitute.org (UV) studies have probed unobscured star formation back to  ∼ 10 (e.g.Bowler et al. 2015Bowler et al. , 2017;;Bouwens et al. 2015;Oesch et al. 2018;Bouwens et al. 2019) and beyond (e.g.Donnan et al. 2023).This is motivated, in part, by the need for better constraints on the physics of reionization; at  > 5, the SFRD determines the contribution of star-forming galaxies (SFGs) to the budget of ionizing photons (see the review by Stark 2016).However, galaxy selections based on unobscured emission (e.g. using the Lyman break or Lyman alpha emission line) are biased towards galaxies that are both young and fairly dust poor (Shibuya 2020), and substantial corrections are required to scale the UV-derived SFRD and bring it into line with infrared (IR)-derived values where the two overlap (at  < 3; Madau & Dickinson 2014).Such corrections are subject to considerable uncertainties on the degree of dust obscuration in galaxies.This is particularly unconstrained in the early Universe (see Ma et al. 2019, and references therein), and possibly underestimated: recent ALMA observations of the dust continuum emission from Lyman Break Galaxies suggests that individual galaxies as distant as  ∼ 8 can harbour significant amounts of dust (Watson et al. 2015;Laporte et al. 2017;Bowler et al. 2018).There are methods of estimating dust corrections from UV observations alone, primarily via the empirical IRX −  relation between the ratio of the infrared luminosity to UV luminosity, IRX, and the UV spectral slope,  (Meurer et al. 1999).However, the considerable scatter in the relation, due to complex geometries, older stellar populations and different intrinsic extinction curves (see Popping et al. 2017;Narayanan et al. 2018), add uncertainty to such corrections.Cosmic variance is also a concern for rest-frame UV studies of the highest redshift galaxies (Trenti & Stiavelli 2008;Driver & Robotham 2010;Ventou et al. 2017), as the depth needed to find these faint objects often comes at the expense of area.
An alternative tracer of star formation is rest-frame far-infrared (FIR) emission; since this is driven by the thermal output of dust heated by young stars, it is used to quantify star formation that is missed by the UV.This is particularly critical around the peak of cosmic star formation ( ∼ 1 − 3), where ∼ 85 per cent of the total star formation is obscured (Dunlop et al. 2017).While the large area surveys enabled by Herschel and the South Pole Telescope have discovered some extreme sources as distant as  ∼ 7 (e.g.Weiss et al. 2013;Strandet et al. 2017;Marrone et al. 2018;Casey et al. 2019), deeper surveys that have the sensitivity to detect more typical sources tend to be limited to small areas of sky (e.g.Aravena et al. 2016;Hatsukade et al. 2016;Dunlop et al. 2017;Franco et al. 2018).Since the redshift distribution of SMGs observed at ∼ 1 mm peaks around  = 2.0 − 2.5 (Chapman et al. 2005;Koprowski et al. 2014;Simpson et al. 2014;Danielson et al. 2017;Stach et al. 2019), constraints on the abundance of dusty sources and their contribution to the SFRD at  > 4 are few, due to the lack of detected sources (Casey et al. 2018).Recent work has shown the promise of untargeted, longer wavelength surveys in isolating higher redshift sources (Zavala et al. 2021;Casey et al. 2021;Cooper et al. 2021;Manning et al. 2022).Zavala et al. (2021) present a large (for mm), 184 arcmin 2 , 2 mm survey, from which they identify 13 sources.Combining their new data with an empirically-based model, they show that dust-obscured star formation dominates the cosmic star formation rate budget to  = 4, dropping to a 35 per cent contribution at  = 5, and 20 − 25 per cent at  = 6 − 7 (broadly in line with Dunlop et al. 2017 andBouwens et al. 2020).
The small field of view imaged in a single ALMA pointing makes extending untargeted sub-millimeter surveys to degree scales technically challenging (Chen et al. 2023).However, building samples of robustly-characterised star-forming galaxies is critical in order to answer key questions in galaxy evolution.These include understanding the timing and drivers of cessation of star-formation in different types of galaxies, which is fundamentally linked to a robust measurement of the cosmic star formation rate density as a function of redshift.We would also like to be able to constrain better the amount of unobscured versus obscured star-formation at different epochs, and how this varies across the galaxy population.This requires large statistical samples that span a broad range of cosmic time (i.e. from  ∼ 0 to well beyond the peak of cosmic star formation), as well as wide fields to overcome cosmic variance.
The new generation of radio interferometers offer a unique opportunity to provide such samples and resolve these issues.Unlike at UV and optical wavelengths, light at radio wavelengths is unaffected by dust obscuration.Non-thermal emission from supernovae at centimetre wavelengths has been shown directly to be a delayed, indirect tracer of star formation (Condon 1992;Cram 1998).Relativistic electrons spiralling in weak magnetic fields emit synchotron radiation, characterised by a smooth spectrum (   ∝   ,where  ∼ −0.7) over a large wavelength range.The broad utility of synchotron emission as a probe of star formation is also supported by the tight far-infrared to radio correlation (FIRC), which has been shown by many to hold over several orders of magnitude in radio luminosity (e.g.van der Kruit 1971;Ivison et al. 2010;Sargent et al. 2010;Bourne et al. 2011;Delhaize et al. 2017;Read et al. 2018;McCheyne et al. 2021).The sensitive, dust-independent nature of radio continuum emission as a star-formation rate tracer has been capitalised on by previous work estimating the star-formation history (e.g.Haarsma et al. 2000;Seymour et al. 2008), including the detailed studies made possible by VLA observations of the COSMOS field (Schinnerer et al. 2007;Sargent et al. 2010;Schinnerer et al. 2010;Delhaize et al. 2017;Novak et al. 2017;Smolčić et al. 2017;Leslie et al. 2020; Van der Vlugt et al. 2021).One key advantage of the radio is its ability to probe star formation and AGN activity in both local and very distant galaxies (radio-loud sources have been discovered at redshifts as high as  = 6.82;McGreer et al. 2006;Saxena et al. 2018;Belladitta et al. 2020;Banados et al. 2021;Ighina et al. 2021).However, to probe the bulk of the star-forming population, we need highly sensitive observations, and the deepest radio surveys typically cover limited areas of sky (2 deg 2 in the case of VLA-COSMOS).
Census studies are now becoming feasible with the current generation of radio telescopes, which can map large sky areas with high sensitivity and good angular resolution in an efficient manner.The International Low Frequency Array (LOFAR; Van Haarlem 2013), is a large array of radio antennas, centred in the Netherlands but with antenna stations around Europe.The large primary beam (full-width at half-maximum 3.8 deg 2 for stations in the Netherlands) enables 10 deg 2 regions of sky to be mapped in a single pointing.Making use of this, the LOFAR Two Metre Sky Survey (LoTSS; Shimwell et al. 2017Shimwell et al. , 2019Shimwell et al. , 2022) ) project is adopting a multi-pronged approach to surveying the Northern sky at radio wavelengths.One strand of this is a wide-field survey of the whole Northern sky at 120 − 168 MHz and ∼ 6 ′′ angular resolution (see Shimwell et al. 2017Shimwell et al. , 2019Shimwell et al. , 2022;;Duncan et al. 2019;Williams et al. 2019).The second strand is a series of deep-field pointings known as the LoTSS Deep Fields.
The LoTSS Deep Fields currently comprise deep observations of three well-studied Northern extragalactic fields: the European Large-Area ISO Survey-North 1 (Elais-N1; Kessler et al. 1996;Oliver et al. 2000), Boötes (Jannuzi & Dey 1999) and the Lockman Hole (Lockman et al. 1986), which are expected to reach eventual depths of ∼ 10Jy/beam rms (Best et al. 2023).The first data release of these Deep Fields data reach ∼ 20 Jy/beam at 150 MHz1 (Duncan et al. 2021;Kondapally et al. 2021;Tasse et al. 2021;Sabater et al. 2021).The radio imaging has been accompanied by a detailed programme of source association and cross-identification (Kondapally et al. 2021), photometric redshift estimation (Duncan et al. 2021) and host galaxy characterisation (Best et al. 2023).In this paper, we perform a detailed study of the radio view of cosmic star formation using the three LOFAR Deep Fields, Elais-N1, Boötes and the Lockman Hole.The sensitivity of our observations is comparable to that of the 3 GHz COSMOS-VLA survey (this reached 2.3 Jy/beam, which is equivalent to ∼ 19 Jy/beam at 150 MHz, assuming a radio spectral index of −0.7).Our data cover a substantially larger area (∼ 26 deg 2 of overlap with ancillary data, across the three fields), providing > 80, 000 radio-identified galaxies with optical counterparts.Together, the multi-wavelength catalogues we have constructed identify diverse populations of galaxies, out to  ∼ 6 (Best et al. 2023).A complementary paper, Kondapally et al. (2022), presents the cosmic history of low-excitation radio galaxies.
The structure of this paper is as follows.In Section 2, we first describe the multi-wavelength coverage of the three fields, and introduce the methods used to match radio sources with multi-wavelength counterparts.We briefly present an overview of the derivation of photometric redshifts and physical properties of the radio sources, as well as the separation of star-forming galaxies from AGN.We describe the methods used to construct luminosity functions and present the evolution of the 150 MHz luminosity functions of star-forming galaxies in Section 3. In Section 4, we construct star formation rate functions (SFRFs) from these luminosity functions, and compare these to SFRFs derived using SED-estimated star SFRs.We also estimate the scatter on the relation between 150 MHz luminosity and star formation rate.In Section 5, we construct the cosmic star formation history, from  ∼ 0 to  ∼ 4. We draw conclusions in Section 6.

THE DATA: PANCHROMATIC OBSERVATIONS OF ELAIS-N1, BOÖTES AND THE LOCKMAN HOLE
In this section, we present an overview of the LOFAR 150 MHz observations as well as the cross-matched UV-FIR photometric catalogues used in this work.

Radio observations with LOFAR
Elais-N1, Boötes and the Lockman Hole were observed by LOFAR at 150 MHz frequency with the HBA (high-band antenna) array, in a series of 8-hour pointings.Observations of Elais-N1 (phase center 16h11m00s +55°00'00"; J2000) took place as part of cycles 0, 2 & 4, with 164 hr of integration time in total.Observations of Boötes (phase center 14h32m00s +34°30'00") were taken in cycles 3 & 8, with 80 hr of integration time in total.The Lockman Hole (phase center 10h47m00s +58°05'00") was observed in cycles 3 & 10, with integration time summing to 112 hr.All three fields were calibrated and imaged using Netherlands-only baselines, which gives rise to an angular resolution of 6 ′′ (note that imaging using international stations is also possible; see Jackson et al. 2022;Morabito et al. 2022b;Sweĳen et al. 2022).The rms sensitivity reached at the pointing centre was 20 Jy/beam for Elais-N1, 32 Jy/beam for Boötes and 22 Jy/beam for the Lockman Hole.Sensitivity decreases further from the pointing center; for each field, the area enclosed by the 30 per cent power point is ∼ 25 deg 2 .A full description of these observations and the radio data reduction process is presented in Tasse et al. (2021) and Sabater et al. (2021).
Radio source extraction was performed on the Stokes I radio image using the Python Blob Detector and Source Finder (PyBDSF; Mohan & Rafferty 2015).The final radio catalogue comprises 84, 862 sources in Elais-N1, 36, 767 sources in Boötes, and 50, 112 sources in the Lockman Hole (Sabater et al. 2021;Tasse et al. 2021).As discussed in the following section, we use a subset of these sources in the following analysis, limiting the sample to radio-identified sources that lie in regions of overlap with key optical and infrared surveys.

Multi-wavelength data
The three fields have different photometric coverage.Here, we review the ancillary catalogues generated and described fully by Kondapally et al. (2021).
In Elais-N1 and the Lockman Hole, -band data are drawn from the Spitzer Adaptation of the Red-sequence Cluster Survey (SpARCS; Muzzin et al. 2009), which used the Canada-France-Hawaii Telescope (CFHT).In Lockman Hole, the SpARCS data also adds images in the , , , and  bands, and there are additional observations from the Red Cluster Sequence Lensing Survey (RCSLenS; Hildebrandt et al. 2016) in the , , , and  bands.In Elais-N1, optical (, , ,  & ) broad-band imaging is provided by the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) 1 survey (PS1; Chambers et al. 2016).Further optical imaging from the Hyper-Suprime-Cam (HSC; Aihara et al. 2018) et al. 2007; see below) with the broad-band filters , , ,  &  , as well as the narrow-band filter NB921.In Boötes, deep optical photometry in the   ,  and  bands is drawn from the NOAO Deep Wide Field Survey (NDWFS; Jannuzi & Dey 1999), with additional -band data from the zBoötes survey (Cool 2007) and additional  spec and  -band imaging from the Large Binocular Telescope (Bian et al. 2013).

Infrared data
At near-infrared (NIR) wavelengths, the UKIDSS-DXS DR10, which used the UK Infrared Telescope (UKIRT), provides  and  band coverage for Elais-N1 and the Lockman Hole.In Boötes, ,  and   data are drawn from the Infrared Boötes Imaging Survey, conducted with NEWFIRM on the Kitt Peak National Observatory Mayall 4-m telescope (Gonzalez 2010).
In the mid-infrared (MIR), Spitzer-IRAC observations at 3.6 m, 4.5 m, 5.8 m and 8.0 m are drawn from the Spitzer Wide-area Infra-Red Extragalactic (SWIRE; Lonsdale et al. 2003), which covers ∼ 8 deg 2 in Elais-N1 and ∼ 11 deg 2 in the Lockman Hole.In these two fields, we also draw data from the Spitzer Extragalactic Representative Volume Survey (SERVS) project (Mauduit et al. 2012), which covers 2.4 deg 2 of Elais-N1 and 5.6 deg 2 in the Lockman Hole, with the 3.6 m and 4.5 m channels, reaching ∼ 1 mag deeper than SWIRE.In Boötes, data at 3.6 m, 4.5 m, 5.8 m and 8.0 m are primarily drawn from the Spitzer Deep, Wide-Field Survey (SD-WFS; Ashby et al. 2009), with additional data from the Decadal IRAC Boötes Survey (M.L.N.Ashby PI, PID 10088).
24 m data are provided by Spitzer-MIPS (Rieke et al. 2004); for Elais-N1 and the Lockman Hole, the data were drawn from the SWIRE survey.At longer FIR wavelengths, data for all three fields were drawn from the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012), which used the Herschel Space Observatory (Pilbratt et al. 2010).Herschel imaging at 100 m and 160 m comes from Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010), and at 250 m, 350 m and 500 m from the Spectral and Photometric Imaging Receiver (SPIRE;Griffin et al. 2010).

Multi-wavelength catalogues
The ultraviolet to mid-IR flux densities described above are compiled in the forced-and matched-aperture, aperture-corrected, multiwavelength catalogues presented by Kondapally et al. (2021).Sources were identified on combined chi-squared stack images; those with signal-to-noise ratio (SNR) less than 3 in all filters were removed from the catalogues.
Full details of the FIR catalogues are provided in McCheyne et al. (2021), so we provide only a summary here.For all fields, FIR flux densities derived by the Herschel Legacy Project (HELP; Shirley et al. 2021) provided the basis of our measurements.As part of HELP, source deblending was performed using the Bayesian tool XID+ (Hurley et al. 2017).Where a LOFAR source could be crossmatched with a HELP catalogue entry within 0.5 ′′ , the HELP fluxes were adopted.Where a LOFAR source had no HELP match, XID+ was re-run with the radio host galaxy position added to the prior list.

A cross-matched radio and photometric catalogue
Full details of the cross-matching of radio-identified sources to the photometric catalogues summarised in Section 2.2 are provided in Kondapally et al. (2021).At the depths reached by the LOFAR imaging in the Deep Fields, radio emission from multiple sources can be incorrectly linked by the PyBDSF source extraction procedure, leading to 'blended sources'.The opposite scenario, extended emission from individual sources being split up into many components, can also occur.In this section, we summarise the approach employed to associate radio and multi-wavelength sources within the LOFAR Deep Fields.This was only performed over the area of sky with the best multi-wavelength data: ∼ 6.74 deg 2 in Elais-N1, ∼ 8.63 deg 2 in Boötes, and ∼ 10.28 deg 2 in the Lockman Hole.
The likelihood ratio (LR) method (de Ruiter et al. 1977;Sutherland & Saunders 1992) uses magnitude and colour information to match radio sources with their optical and IR counterparts, as described in Kondapally et al. (2021).This technique yields robust associations for the majority of radio sources (80 − 85 per cent of sources in each field).For sources with extended or complex radio emission, sources were visually inspected using LOFAR Galaxy Zoo (LGZ; see also Williams et al. 2019 for an earlier use of the same framework), a private Zooniverse project accessible to members of the LOFAR consortium.Other sources, primarily those that were potential radio blends, were sent for identification by an 'expert user', who de-blended the PyBDSF source using the PyBDSF Gaussian components.
Within the overlap region of PanSTARRS, UKIDSS and SWIRE in Elais-N1, there are 1, 470, 968 optically-detected sources and 31, 610 radio-detected sources.After cross-matching to optical/NIR counterparts and the flagging of spurious sources, there were 30, 839 sources.Within the overlap region of the NDWFS and the SDWFS in Boötes, there are 1, 911, 929 optically-detected sources, 19, 179 radio-detected sources, and 18, 579 sources in the final cross-matched catalogue.Within the overlap region of the SpARCS r-band and SWIRE in the Lockman Hole, there are 1, 906, 317 optically-detected sources, 31, 162 radio-detected sources, and 30, 402 sources in the cross-matched catalogue.As discussed in some detail by Kondapally et al. (2021), 2−3 per cent of sources (771 in Elais-N1, 600 in Boötes, and 760 in the Lockman Hole) have no robust optical/NIR ID.Some of these sources have FIR counterparts, and are likely dusty AGN and SFGs.Note that Novak et al. (2017) reported that ∼ 4 per cent of VLA-COSMOS sources are optically faint shortward of -band for a similar radio-selected population.NIR/optically-dark sources have also been reported as common in sub-millimeter selected samples (e.g. Simpson et al. 2014Simpson et al. , 2017;;Franco et al. 2018;Wang et al. 2019;Dudzeviciute et al. 2020;Smail et al. 2021).Kondapally et al. (2021) investigated the LOFAR-detected radio sources without IDs, and found that these are likely dominated by  ∼ 2 − 4 AGN.Since we focus on star forming galaxies in this paper, they are unlikely to have much effect on our analysis.

Deriving photometric redshifts of LOFAR-identified sources
The process of deriving photometric redshifts for sources in all three fields is described fully in Duncan et al. (2021) and we review only the most important details here.Both template fitting and machine learning techniques were used to derive photometric redshifts for both the radio-selected LOFAR sources, and the full optical catalogue (∼ 5 million sources in total).This 'hybrid' approach was shown by Duncan et al. (2018b) to improve upon traditional template fitting, particularly for intermediate redshift AGN, which had proved a challenging population to obtain redshifts for (Duncan et al. 2018a).Duncan et al. (2021) showed that outlier fractions, defined as | phot −  spec |/(1 +  spec ) > 0.15, range from 1.5 to 1.8 per cent for galaxies, and from 18 to 22 per cent for optical, IR and X-ray selected AGN.In this paper, we make use of the photometric redshift posteriors for sources without spectroscopic redshifts (this is the majority; only ∼ 8.6 per cent of all LoTSS Deep Fields radio sources have spectroscopic redshifts).Since we exclude the majority of AGN from our analysis, these redshifts are reasonably robust (≲ 2 per cent outlier fractions).

Classification of star-forming galaxies and AGN
Reliable classifications of star-forming sources and AGN are essential for this study.Best et al. (2023) present the multiple methods used to identify AGN from emission in different wavebands, which we summarise here.A combination of spectral energy distribution (SED) fitting codes was used to identify the majority of radiativemode AGN (i.e.those identifiable via their optical/IR/X-ray emission).Each radio-identified source was fitted with four different SED fitting codes: MAGPHYS (da Cunha et al. 2008), BAGPIPES (Carnall et al. 2018(Carnall et al. , 2019))  Million Quasar catalogue, which is mainly based on the Sloan Digital Sky Survey (Alam et al. 2015).X-ray AGN (defined by X-rayto-optical flux ratios or hardness ratios) were identified in Boötes thanks to the deep X-ray data provided by Chandra, as part of the X-Boötes survey of NDWFS (Kenter et al. 2005).In Elais-N1 and the Lockman Hole, the Second ROSAT All-Sky Survey (Boller et al. 2016) and the XMM-Newton Survey provide shallower data, enabling the identification of brighter X-ray sources.
Radio-selected AGN were identified based on their excess emission at 150 MHz.For each source, the radio emission expected for a purely star-forming source with a given SED-derived SFR (see Section 2.6) was determined using the 'ridgeline' relation between SFR and  150 MHz derived by Best et al. (2023); note that the 'ridgeline' agrees well with the relation derived by Smith et al. (2021).Then, excess radio emission relative to this expected value was calculated.Sources that were either offset from the relation by more than 0.7 dex, or had extended (> 80 kpc) radio emission, are likely not powered solely by star formation, and were hence classified as radio-selected AGN.
In our final classifications, star-forming galaxies were defined as those without AGN signatures in optical/IR/X-ray, and without excess radio emission.Sources with excess radio emission form the radio-loud AGN classes.In this paper, we focus on galaxies without radio excess, i.e. galaxies classified as either star-forming or radioquiet AGN.These comprise 77% of the total LoTSS sample (81% of sources with optical/NIR counterparts and classifications; see Best et al. 2023).For the radio-quiet AGN, star-formation is expected to drive the majority of the radio continuum emission detected by LOFAR (though there will be a minority of cases where a weak jet dominates; see Macfarlane et al. 2021 for a model of quasar radio luminosity that comprises contribution from star formation and an AGN jet.They concluded that the jet-launching mechanism operates in all quasars, but with different efficiency.See also Gürkan et al. 2019).
Our final sample comprises 55, 581 radio sources within the redshift range 0.1 <  < 5.7 for which radio emission is dominated by star formation.21, 638 of these are in Elais-N1, 12, 787 are in Boötes, and 21, 156 are in the Lockman Hole.

Stellar masses and star formation rates for star-forming galaxies
As noted in Section 2.5, four different SED fitting codes are used to fit all radio-selected sources.SED fitting provides an alternative way to estimate galaxy SFRs, compared with single wavelength flux calibrations.'Energy balance' based SED fitting ensures that the unobscured and obscured components of galaxy emission are fitted simultaneously and self-consistently with a combination of simple stellar populations, a dust model, and models of star formation history and galactic chemical evolution.Physical parameters may then be derived from the best-fitting combination of models.These codes have been tested on simulated galaxies (e.g.Hayward & Smith 2014;Dudzeviciute et al. 2020).Each of the four SED fitting codes provides estimates for the physical properties of the sources, including stellar mass and star-formation rates.Best et al. (2023) describe the process by which these SED fitting outputs are combined to generate 'consensus' estimates.For all sources that were not classified as radiative-mode AGN, stellar mass and SFR were both generally calculated using the mean of the logarithm of the values derived using MAGPHYS and BAGPIPES, providing both codes yielded a statistically acceptable fit (note that this was the case for ∼ 90 per cent of these sources).Where one fit was unacceptable, the estimate from the acceptable fit was adopted.For the radio-quiet AGN in our sample, stellar masses and SFRs were derived using more appropriate SED fitting techniques, which included a component of emission from the AGN.For the vast majority of these (> 94 per cent), stellar masses and SFRs were derived using CIGALE (see Best et al. 2023 for full details of the small number of cases where AGNfitter was preferred).
The distributions of 150 MHz luminosity, star formation rate and stellar mass, for our selected population of star-forming galaxies and radio-quiet AGN, are shown in Figure 1.

Constructing the 150 MHz luminosity function
With our sample of star-forming galaxies and radio-quiet AGN inhand, we construct 150 MHz luminosity functions at a range of redshifts.Approximately 90 per cent of the radio sources in our sample are classified as 'pure' SFGs, but we include radio-quiet AGN to perform a complete census of star formation (see also Bonato et al. 2021, who construct radio luminosity functions for the separate populations).Rest-frame 150 MHz luminosities are calculated using the following formula, where  is the observed-frame frequency, 150 MHz,  is the radio spectral index (we assume  = −0.7),and   is the flux density at the observed frequency. is the source redshift (we use the spectroscopic redshift, if it exists, and if not, the photometric redshift, as described in Section 2.4), and   is the corresponding luminosity distance: (1) We aim to calculate the space density of sources per luminosity bin, the 'luminosity function', as a function of redshift.First, we divide the sources into eleven wide redshift bins.We place the observed LOFAR sources into these redshift bins, and calculate a rest-frame 150 MHz luminosity for each source.We bin sources further, by restframe 150 MHz luminosity, and then for each luminosity bin at each redshift, we derive the 150 MHz luminosity function, Φ(, ), using the non-parametric 1/ max method (Schmidt 1968): where Δ log 10  is the width of the luminosity bin and  max is the volume within the redshift bin over which a source with a given restframe 150 MHz luminosity would be observable, given the sensitivity limits of the survey.This is particularly important for a radio survey such as ours, where the depth is not uniform across each individual field, nor between fields. max is calculated for each source using: where  ()  is the whole-sky co-moving volume in the redshift range [,  + ] and  (, ) is the fractional area over which a source of that flux density would have been detected with 5 signal-tonoise.We perform the integral numerically, by dividing the wide redshift bin (with edges  min and  max ) into narrow redshift slices of size Δ() = 0.0001.For each narrow redshift slice, we calculate the 150 MHz flux density that we would observe if the source were located at the centre of that redshift slice (using its known rest-frame luminosity).We then calculate  (, ) for that flux density using: where Ω[()] is the solid angle over which a source with flux density  can be detected at 5, and  radio [()] is the radio completeness as a function of flux density, as derived in Section 3.2. photometric () is an order-unity correction factor which accounts for inaccuracies (such as aliasing) in the photometric redshifts, as derived in Section 3.3.This enables us to fold in the spatially-varying radio depths over the fields, as well as uncertainties in photometric redshifts.
Repeating the process for each of the wide redshift bins enables us to construct luminosity functions from  ∼ 0 to  ∼ 5.

Radio completeness corrections
Completeness was calculated as a function of source flux density by simulating the source detection rate using the same techniques used to identify real sources.Approximately 100, 000 mock Gaussian sources of known angular extent and known 150 MHz total  As described in Section 3.2, the completeness curves were calculated using source extraction of mock input sources from the radio images, using the same PyBDSF parameters as were used for the extraction of real sources.A realistic source size distribution was applied, based on the observed size distribution of star-forming galaxies with flux densities in the range 1 − 5 mJy, where completeness is high.As expected based on total observing times, Elais-N1 has the deepest radio coverage, followed by the Lockman Hole, and then Boötes.Note that due to a very small number of injected mock sources overlapping with real sources and other injected sources, the completeness never quite reaches 100 per cent.
flux density were placed into regions of the LOFAR image that are covered by the optical imaging (accounting for masked regions).In practice, we insert 1, 000 mock sources at a time, a small number compared to the number of real sources in the LOFAR image, and repeat ∼ 1, 000 times; this is to avoid bias due to confusion.We used a continuous distribution of galaxy major axis size (from 6 ′′ to 30 ′′ ), and flux density (from 0.1 mJy to 40 mJy).For each source, minor axis size was drawn from the distribution of minor axes of observed sources with roughly the same major axis size.Each mock source was separated from neighbours by at least twice its major axis size, so that sources did not overlap.PyBDSF was run on the image with inserted sources, using the same parameter choices as were adopted for the real radio source extraction.A source was deemed to be 'recovered' if PyBDSF identified a source with flux density greater than 5 of the local rms, within 2 ′′ of the position of the inserted source.The fraction of recovered to input sources was then calculated to characterise the completeness as a function of flux density and source size.To derive a single completeness curve per field, we then folded in our best estimate of the 'true' source size distribution.This was derived from the size distribution of the star-forming galaxies in our sample where completeness reaches ∼ 100 per cent (source with flux densities in the range 1 − 5 mJy).We are careful not to use only the very brightest sources, which will be biased in favour of nearby, spatially extended sources.Instead, we use the size distribution of sources at the point where the completeness curve flattens; in this way, we obtain the most similar size distribution to our sample, which is dominated by compact sources (note that intrinsically compact sources can have sizes larger than the beam due to the calibration).The size distribution of star-forming sources at 0.5 1.0 1.5 2.0 2.5 3.0 Figure 3: Corrections applied to the number densities of radio sources derived in each redshift bin, for each field.To derive these corrections, a photometric redshift was drawn for each source from a flat redshift distribution between  1,MIN and  1,MAX .The number of sources within each redshift bin studied was then compared to the number of sources within the bin that is derived using  BEST .The shaded regions are derived using 10, 000 bootstrapped samples.We apply these corrections to the number densities derived for each field, up to and including the  = 2.5 − 3.3 bin.At higher redshifts, these corrections become less reliable due to long tails of photometric redshift probability distributions that peak at lower redshift.Hence, we do not apply corrections at  > 3.3.
1 − 5 mJy is dominated by compact sources, with 84 per cent having major axis sizes < 9 ′′ and a tail to larger values.We derive a single completeness curve, for each field (see Figure 2) and provide these in Table A1.The 50 per cent completeness values are 128 Jy for Elais-N1, 246 Jy for Boötes, and 180 Jy for the Lockman Hole.Equivalent curves were also derived assuming a source-size distribution appropriate for AGN, and are presented in (Kondapally et al. 2022).Radio completeness corrections were applied when constructing luminosity functions, as described in Section 3.1.We discard luminosity bins where applying the derived radio completeness correction results in a change of > 0.5 dex.

Incorporating uncertainties in photometric redshifts
When deriving luminosity functions, the redshift used for each source,  BEST , was the spectroscopic redshift where available, or otherwise the median of the primary photometric redshift solution.As discussed in Section 2.4, the derivation of these redshifts combines template fitting and machine learning techniques, which yield particularly reliable results for the star-forming galaxies studied in this paper.However, as seen in Figure 1, there is a dip in galaxy numbers at  ∼ 1.5, which may be driven by uncertainties in photometric redshifts.This effect is particularly marked in Elais-N1 and the Lockman Hole, which lack -band data.The -band is important in characterising the observed-frame wavelength of the 4000Å break; for galaxies in the wavelength range 1.3 <  < 2.0, the 4000Å break falls between the -band and the -band.For these sources, -band data constrains the flatter spectrum above the break, and hence enables more accurate redshifts to be derived.In the absence of the -band, some photometric redshift aliasing can occur, and uncertainties on the derived photometric redshifts are larger.Here, we characterise the impact of these uncertainties on derived number densities in each of the three fields.
Fully correcting for the uncertainties in photometric redshifts and their impact on the derived luminosity function would be complex.Assuming that the photometric redshift probability distribution of each source is robust, we could sample this finely and repeat the multi-code SED fitting, source classification and determination of a consensus SFR at each redshift for each source.When constructing the luminosity function, we could then draw bootstrapped samples, with the redshifts of each source drawn from its photometric redshift distribution.Given the substantial effort and computational expense required to derive source properties given just the best estimate redshift (see Best et al. 2023), this would be impractical.There is also little evidence that the photometric redshift uncertainties are dependent on the radio luminosity of the source (which could require us to apply a luminosity-dependent correction).Instead we aim to take a simpler approach that will provide an approximate correction factor to the entire LF at each redshift for each field,  photometric ().This method accounts for aliasing in the photometric redshifts, under the assumption that there are no radio luminosity-dependent effects, and that source classifications and physical properties are not systematically changed by redshift errors.
Our approach involves perturbing the redshift of each source according to its photometric redshift distribution, and calculating the change in numbers of galaxies that fall within each redshift bin studied.Instead of assuming  TRUE =  BEST , we draw a photometric redshift for each source from a flat redshift distribution between its  1,MIN and  1,MAX .We then compare the new number of sources within each redshift bin to the number of sources within the bin that is derived using  BEST .We repeat this process 10, 000 times, using different random values between  1,MIN and  1,MAX .This enables us to derive correction factors to the number densities of sources within each field and redshift range, as shown in Figure 3 and tabulated in Table B1.While most redshift bins have corrections ∼ 1 (i.e.no net gain or loss of sources, as the same number of sources are scattered into a given redshift bin as out of it), we see larger corrections required for Elais-N1 and the Lockman Hole (the fields that lack -band data) at  ∼ 1.5.More sources get corrected into the  ∼ 1.5 redshift bin than out of it; this is because sources have  BEST values that are pushed to higher/lower redshifts, but with large uncertainties that cover the  ∼ 1.5 bin.This leads to an upward correction to the luminosity function of ∼ 50 per cent for Elais-N1 and the Lockman Hole at  ∼ 1.5, with corresponding decreases for the neighbouring higher and lower redshift bins.We apply the derived corrections (as a multiplicative factor to ) up to and including the  = 2.5 − 3.3 bin.At higher redshifts, these corrections become less reliable due to long tails of photometric redshift probability distributions that peak at lower redshift and are therefore not applied.

Parametric fits to the radio luminosity function
To characterise the evolution of the 150 MHz luminosity function, we fit each of the derived LFs with a parameterised expression.Numerous studies of radio luminosity functions of star-forming galaxies (Saunders et al. 1990;Best et al. 2005;Novak et al. 2017) have adopted the following parametrisation, first used by Sandage et al. (1979): where  ★ provides the normalisation,  ★ is the luminosity at the turnover,  is the faint end slope, and  describes the steepness at the bright end.
We repeat this process to derive  and  self-consistently from our own LOFAR Deep Fields data.We fit the low redshift ( = 0.03−0.3)luminosity function, as shown in Figure 4. We exclude the first luminosity bin from the fit due to the larger (> 0.3 dex) completeness corrections for the faintest sources.We also exclude the final luminosity bin due to the increased potential importance of misclassification for the brightest sources; as shown in Figure 4, radio AGN dominate the whole radio sample here.We repeat this fitting process for 18 combinations of luminosity bin size and position, and average over the results to derive the following best fit parameters: log 10 ( ★ / Mpc −3 dex −1 ) = −2.46 ± 0.01,  = 0.49 ± 0.01,  = 1.12 ± 0.01 and log 10 ( ★ /W Hz −1 ) = 22.40 +0.02 −0.03 .We fix  = 0.49 and  = 1.12 for the remainder of this work.
We overplot several measurements from the literature in Figure 4. Kondapally et al. (2022) derived radio luminosity functions of radio-selected AGN in the three LoTSS Deep Fields as a function of redshift.Their  = 0.03 − 0.3 measurement is plotted here.We note the very different shapes of the luminosity functions for the two populations.As expected, AGN dominate the source counts at the highest luminosities.At  150 MHz ∼ 10 23.75 W Hz −1 , there are approximately equal contributions of AGN and SFGs to the source counts.We also show luminosity functions constructed for the starforming galaxy population by Best & Heckman (2012) and Sabater et al. (2019).While in good agreement with each other, these data lie below ours at all but the brightest luminosities.Partly, this is due to the lack of radio completeness corrections in the earlier work.Differing redshift distributions of the samples will also contribute: for our sample, the median redshift of sources included in the  = 0.03 − 0.3 subsample is  median = 0.20.The median redshift of star-forming sources using in both Best & Heckman (2012) and Sabater et al. (2019) will be lower, due to less sensitive radio imaging.Coupled with the strong redshift evolution seen for these samples (see Figure 6), this will naturally lead to a small offset in derived luminosity functions.The luminosity functions derived for each of the LOFAR Deep Fields, as well as for the combination of the three, are shown in Figure 5.These luminosity functions have been corrected for radio completeness and uncertainties in photometric redshifts.The fields display excellent agreement with each other, indicating that cosmic variance effects are minimal (as expected over such large areas).As a sanity check, we confirm that our estimates agree well with data from Bonato et al. ( 2021), who measured luminosity functions of star-forming galaxies and radio quiet AGN out to  = 2.8 using Deep Fields data in the Lockman Hole.Their measurements are consistent with modelled populations in the Tiered Radio Extragalactic Continuum Simulation (T-RECS; Bonaldi et al. 2019).As shown in Figure 5, we also find very good agreement with the luminosity functions derived for VLA-COSMOS by Novak et al. (2017), once luminosities are scaled from 3 GHz, under the assumption of a fixed radio spectral index (we use  = −0.7),apart from at the very highest redshifts studied, where number densities in COSMOS drop more significantly.In Figure 6, we show the derived luminosity functions at all redshifts on a single panel, for easier comparison.
At each redshift, we fit the radio luminosity function for each of the three fields, as well as for all fields combined, with the parametrisation presented in Equation 5, fixing  = 0.49 and  = 1.12.The inferred  ★ and  ★ values (displayed in Figure 7 and Table 2) are in good agreement between the three fields, except at the very highest redshifts, where  ★ falls well below the sensitivity limits of the radio data. ★ increases monotonically out to at least  ∼ 3, displaying an evolution of > 1.5 dex between  ∼ 0.25 and  ∼ 3.  ★ remains roughly constant back to  ∼ 0.8 but then falls steeply at higher redshifts, decreasing by > 1 dex between  ∼ 0.7 and  ∼ 3.

The star formation rate function constructed using different SFR estimates
In this section, we derive the star formation rate function (SFRF) using multiple SFR estimates.In principle, we can transform the radio LFs to SFRFs using a previously-calibrated relation between   We construct an SFR function using SFR estimates derived from both  150MHz and stellar mass as follows.For each galaxy, we input the 'consensus' stellar mass estimate provided by Best et al. (2023) and the radio luminosity into Equation 7 to obtain a new SFR estimate.We then construct the SFR function in a similar way to the luminosity function, applying radio luminosity and photometric uncertainty corrections on a source-by-source basis.We overlay the SFR function derived for each redshift bin in yellow in Figure 8.The SFR functions diverge from those derived from radio luminosity without stellar mass only at very high SFRs and at redshifts beyond the range that was used to derive the  150 MHz − SFR relation.
We also construct SFR functions using the consensus SFR estimates presented by (Best et al. 2023, blue lines); as described in Section 2.6, these were derived using SED fits to the multi-wavelength data, rather than a single wavelength indicator.We compare the SFRFs derived from the two radio luminosity calibrations (red and yellow lines) to those derived using consensus SFR estimates (blue lines).As seen most clearly at low redshifts, below the break of the function (SFR ★ ), the SFRF estimates are in reasonable agreement.However, they can differ substantially at the highest SFRs, in some cases by an order of magnitude.This is despite the  150 MHz − SFR conversion also being derived from the LOFAR Deep Fields data and SED fits.At first glance, it is worrying that we obtain such different SFRFs when using SFRs derived using different methods, given that the parent radio samples used are the same.The difference is largest at high SFRs, so is of particular concern for studies like ours, where most of the SFRF data points are above ∼ SFR ★ at all but the lowest redshifts studied.Since the star formation rate density is derived by integrating the SFR function, this could have implications for the normalisation and the shape of the inferred cosmic star formation rate density-redshift relation.In Section 4.2, we show that this effect arises due to scatter in the  150 MHz − SFR relation, which we have not accounted for thus far, and we develop a method to correct for this bias.

The impact of the calibration between radio luminosity and SFR
As discussed in Section 4.1, the derived star formation rate function depends strongly on the method used to estimate star formation rates.Above ∼ SFR ★ , the SFR function is highly dependent on the method used to infer SFRs.We expect that this is, at least in part, due to a combination of Eddington bias (Eddington 1913) and the scatter in the L 150 MHz − SFR relation.In this section, we explore the effects of different amounts of scatter on the derived star formation rate function.Using a simple simulation, we demonstrate the magnitude of the bias and derive correction factors.We begin by generating ∼ 300 million mock sources with star formation rates drawn from a Saunders et al. (1990) function with default values of SFR ★ = 0 and  ★ = 0. We set  = 0.38 (the average value of  from fits to the consensus estimate-derived SFR function at  ≲ 1).The modelled base SFR function is shown in black in the left-hand panel of Figure 9.We then simulate the radio luminosities of the mock sources using Equation 6, adding scatter drawn from a Gaussian distribution with Σ = [0.2,0.3, 0.4, 0.5] dex, but truncated at 0.7 dex.This enables us to mimic the star-forming galaxy selection applied in Best et al. (2023), where radio-excess sources with luminosities exceeding the 'ridgeline' value by > 0.7 dex are classified as radio-loud AGN and thus excluded from the SFG sample.Finally, we convert the modelled radio luminosities back to star formation rates, assuming Equation 6.This yields a sample of sources with estimates of star formation rates that are perturbed from their original values according to the modelled scatter on the L 150 MHz − SFR relation.
We construct SFR functions for each instance of modelled scatter.We plot these in colour in Figure 9, also showing the true input SFR function in black.The modelled SFR functions differ substantially above ∼ SFR ★ : larger scatter in the relation causes the observed number density of galaxies to increase at the bright end, causing a gentler fall-off of the exponential (e.g.see the exaggerated scatter values of 0.5 shown by the green lines).Because of the steepness of the original modelled SFR function, the number of sources 'scattered up' to those SFRs (coloured lines) can vastly exceed (by up to ∼ 1 dex at the highest SFRs) the genuine number of sources.This will have a significant impact on the SFR function derived from L 150 MHz measurements, which we need to correct for.The differences we see are qualitatively in line with the differences seen in our observational estimates of SFR functions in Figure 8: as described in Section 5, the SFR functions constructed using consensus estimates tend to lie below those inferred directly from  150 MHz SFR ≳ SFR ★ .
In the right-hand panel of Figure 9, we plot the offsets between the original modelled SFR function and those modelled using various values of scatter in the  150 MHz −SFR relation (i.e.deviation of each of the coloured lines from the black line).We use these modelled offsets, in combination with the differences between our multiple estimates of the SFR function, to constrain the true scatter on the  150 MHz − SFR relation.We use the SFR functions derived from the consensus estimates as the 'truth' (analogous to the black SFR function in the left-hand panel of Figure 9); of course, SFRs derived from SED fitting have their own uncertainties, but these are not dependent on the scatter in the  150 MHz − SFR relation.We compare them to the SFRs derived from the radio luminosity functions, using a single  150 MHz − SFR scaling with no scatter.Our data are shown in the right-hand panel of Figure 9 (grey points).They are broadly consistent with a scatter of ∼ 0.3 dex on the  150 MHz − SFR.This is in good agreement with the results of Smith et al. (2021), who found significant scatter only at SFR > 0. They estimated the scatter to be 0.31 ± 0.01 dex at 1 < log 10 (SFR/M ⊙ yr −1 ) < 2.
As shown, the measured radio luminosity function (and the SFR function derived by scaling this with a single  150 MHz − SFR calibration) will have an artificially shallower bright end slope due to the ∼ 0.3 dex scatter in the  150 MHz − SFR relation.This will lead to an overestimation of the cosmic star formation rate density.By integrating the modelled SFR function derived using different values of scatter, and comparing to the integrated base SFR function, we estimate the degree of boosting of the SFRD.We hence derive correction factors that can be applied to values of SFRD derived using a single  150 MHz − SFR calibration.We derive the multiplicative correction factor CORR SFRD = [1.0,0.96, 0.93] for scatter Σ = [0.2,0.3, 0.4] dex.Without the exclusion of sources classified as radio-loud AGN, this correction would need to be larger.In Section 5, we correct our estimates of cosmic star formation rate density using a scaling factor of 0.96 ± 0.04, assuming a ∼ 0.3 dex scatter in the  150 MHz − SFR relation.

THE COSMIC STAR FORMATION DENSITY HISTORY
The cosmic star formation rate density can be calculated at a given epoch, z, by integrating the radio luminosity function as follows: At each redshift, we derive estimates of the SFRD for each of the three fields and for all fields combined, by integrating the parameterised radio luminosity function presented in Equation 5with the appropriate best-fitting parameters.We adopt a lower luminosity limit  min = 0.03  ★ , and an upper luminosity limit of  max = 10 28 W Hz −1 (an order of magnitude brighter than our brightest luminosity bin).In Appendix E we present tests showing that the derived SFRD is robust to our choice of integration limits.CORR SFRD = 0.96 ± 0.04 is the correction factor derived in Section 4.2, which accounts for the potential bias in SFRD due to the scatter in the  150 MHz − SFR relation.We present our estimates for the star formation rate density at 0 ≲  ≲ 4 in Figure 10.Estimates using the three fields individually show good consistency, with larger differences between fields at  ∼ 1.5, where the lack of -band data in Elais-N1 and the Lockman Hole drives particularly large uncertainties in photometric redshifts (see Section 3.3 and Figure 3).Although we have attempted to resolve this by correcting for these uncertainties, the 'corrected' data   C1.points for the SFRD at  ∼ 1.5 for Elais-N1 and the Lockman Hole disagree with the estimate using Boötes and with the neighbouring redshift bins (see Figure 10 and Table C1).Because of this, we believe that these data are unreliable and we adopt Boötes data as our best estimate at  ∼ 1.5.There is also ∼ 1  discrepancy between the estimates derived from different fields at  ∼ 3 − 4.This arises due to the slightly different depths of the radio data.LOFAR coverage of Boötes is the shallowest, and fitting without the faintest luminosity data point leads to a slightly higher  ★ and lower  ★ being favoured.
In Figure 11 we plot our derived star formation rate density history from all three fields combined (see blue circles).We tabulate these estimates in Table 2.At most of the redshifts studied, our best estimate comes from integrating the luminosity function of the three fields combined: this gives the greatest numbers of sources, and enables us to average over any potential bias due to cosmic variance.At  ∼ 1.5, our best-estimate SFRD comes from measurements taken in the Boötes field, due to increased uncertainties on the photometric redshifts of sources in Elais-N1 and the Lockman Hole.

The functional form of the SFRD
Following Hopkins & Beacom (2006) and Behroozi et al. (2013), we use the emcee fitting code Foreman-Mackey et al. (2013) to fit the following functional form to our data, using the SFRD derived from all three LOFAR fields: We also include an additional radio-derived SFRD measurement at  = 0.043 (Mauch & Sadler 2007) to anchor the fit at  ∼ 0 (note that we apply a correction to account for differing assumed IMFs).We adopt this measurement because the LoTSS Deep Fields do not probe enough volume to provide secure cosmic SFRD measurements at  ∼ Mauch & Sadler 07 Figure 11: The cosmic formation rate density, derived using the full sample of star-forming galaxies from LOFAR, with, literature data overlaid.Our new results using the three LOFAR deep fields combined are shown in navy (solid circles).At  ∼ 1.5, our best estimate of the SFRD comes from measurements taken in the Boötes field, due to increased uncertainties on the photometric redshifts of sources in Elais-N1 and the Lockman Hole.The navy shaded region shows our best estimate for the evolution of the SFRD from  ∼ 0.2 to  ∼ 4, the 1  posterior of our fit to the SFRD derived from all three fields.We also include a  ∼ 0 measurement from Mauch & Sadler (2007), shown in dark grey, in the fit.Our data are bracketed by previously-derived fits to data by Hopkins & Beacom (2006) (above) and Madau & Dickinson (2014) (below); the Behroozi et al. (2013) fit shows best agreement.Coloured symbols represent a selection of estimates from the literature, derived using widely varying sample selections.These include UV-selected (Bouwens et al. 2015; grey squares), optical/NIR-selected (blue squares; Driver et al. 2018), FIR-selected (green hexagons and orange stars; Gruppioni et al. 2013 andDunlop et al. 2017, respectively), and radio continuum-selected (red, pink, purple and brown symbols); Karim et al. 2011;Novak et al. 2017;Leslie et al. 2020;Enia et al. 2022).Overall, there is considerable scatter in measurements across the literature, with disagreements of > 0.4 dex at any given redshift.Our data show good agreement with the FIR-based measurements from Dunlop et al. (2017) and radio continuum-based analyses of Leslie et al. (2020) and Enia et al. (2022).Measurements made by Driver et al. (2018) fall below our estimates at  ≲ 3.5.0; in comparison, the (Mauch & Sadler 2007) measurement is made over a ∼ 300 times larger sky area (7000 deg 2 ).Although too shallow for high-redshift studies like ours, their data are deep enough to constrain the SFRD at very low redshift.We derive  = −0.89+0.08 −0.07 ,  = 0.22 ± 0.04, log 10  = −0.76± 0.05 and  0 = 1.22 ± 0.15.
The best-fitting SFRD is overplotted and compared to previous fits in Figure 11.Our data and fit are broadly consistent with previous fits to older data, lying approximately at or below the estimate of with the fit presented by Behroozi et al. (2013); their fitted form lies within ∼ 1  of our fit at 0.7 ≲  ≲ 4.0.At lower redshifts, we measure a shallower evolution of the SFRD.

Comparison to literature data
In this section, we compare our SFRD measurements to estimates from individual studies of star-forming galaxies selected using a variety of methods.While there exist a vast number of low-z measurements in the literature, we focus on those that reach out to high redshift, including several that have been published since the compilation of Madau & Dickinson (2014), and those that use longer-wavelength SFR estimators.As shown in Figure 11, previous estimates show a consistent general form, with the SFRD increasing between  = 0 and  ∼ 2 and then declining towards higher redshift.However, exact measurements disagree by > 0.4 dex at any given redshift.Here, we compare our LOFAR results to several previous measurements in detail.Bouwens et al. (2015) identified galaxies at  ∼ 4 − 10 in the HST legacy fields using the Lyman break technique.At  = 3.8, their sample consists of −band dropouts.They estimated the SFRD from the raw UV luminosities and also using a correction for dust attenuation based on the IRX- relation.At  = 3.8, the mean dust extinction,  UV , is 2.4.We correct their estimates from a Salpeter (1955) to Chabrier (2003) IMF, and plot both corrected and uncorrected estimates in grey.Our dust-independent estimate lies in between the dust-uncorrected and dust-corrected values, which may indicate that the necessary correction for dust attenuation was overestimated.Driver et al. (2018) combined -band selected galaxies from GAMA (Driver et al. 2011;Liske et al. 2015), -band selected galaxies from G10-COSMOS (Davies et al. 2015;Andrews et al. 2017) and 1.6 m-selected galaxies from 3D-HST (Momcheva et al. 2016) in their analysis.They derived SFRs using MAGPHYS, with various combinations of multi-wavelength data.Their SFRD estimates (blue squares) fall below ours at all redshifts apart from  ∼ 3.5.Interestingly, there is a particular increase in the offset between  = 1.6 (offset 0.18 dex) and  = 1.975 (offset 0.31 dex).Between these two redshifts, the sample changes from including G10-COSMOS sources to a 3D-HST-only sample for which SFRs are derived without FIR data.The large offsets between our dust-independent estimate and theirs at  ∼ 2 and  ∼ 2.4 suggests that dust-obscured star formation is significant at these redshifts and that, in the absence of FIR data, they are under-predicting the SFRD.Our estimates are in nearperfect agreement at  ∼ 3.5 − 4.This implies that the contribution of dust-obscured star formation to the total SFRD is lower by then; this is broadly in agreement with Dunlop et al. (2017) and Zavala et al. (2021).
SFRD estimates using FIR-based studies are generally in better agreement with ours.We overplot estimates from Herschel-selected samples (Gruppioni et al. 2013;green hexagons).These agree with our best-fitting line to within ∼ 0.15 dex at  < 3.5.At  ∼ 3.6, their error bars are very large but remain consistent with our estimate.We also compare to measurements from Dunlop et al. (2017) (orange stars), who combined direct detections and stacking of ALMA imaging in the Hubble Ultra Deep Field with HST-derived SFR measurements from the rest-frame UV to estimate the total SFRD.Our measurements are consistent with theirs.
Finally, we compare to other work based on radio continuum emission.Karim et al. (2011) stacked 1.4 GHz data from VLA-COSMOS, at positions of a 3.6 m-selected sample of > 10 5 galaxies.Leslie et al. (2020) built on this work, stacking 3 GHz data within the same field and refining the source selection to a fully mass-selected sam-ple (using   -band data for galaxies at  < 2.5 and 3.6 m data at higher redshifts).Our measurements are fully consistent with those of Leslie et al. (2020), but are discrepant with those of Karim et al. (2011) at  ≲ 1.5.Leslie et al. (2020) note that Karim et al. (2011) use a different  radio −SFR calibration, which yields lower SFRs, but differences are also likely driven by the deeper parent catalogue used by Leslie et al. (2020).Enia et al. (2022) constructed a 1.4 GHz-selected sample using VLA observations in GOODS-N to measure the SFRD to  ∼ 3.5.These measurements are also fully consistent with ours and with those of Leslie et al. (2020).We note, though, that our considerably larger area (∼ 26 deg 2 compared to ∼ 2 deg 2 for VLA-COSMOS and 171 arcmin 2 for GOODS-N) enables much tighter constraints.The measurements of Novak et al. (2017) are below ours (by up to 0.4 dex) and also below those of Enia et al. (2022) and Leslie et al. (2020).This is perhaps surprising, given that like Enia et al. (2022) and this work, Novak et al. (2017) use a radio-selected sample and the luminosity functions they derive are in good agreement with ours (once scaled to the same rest-frame wavelength; see Figure 5).Enia et al. (2022) suggest that discrepancies between their measurements and those of Novak et al. (2017) might be due to the shallower faint end slope used by Novak et al. (2017), but our fitted faint end slope is actually slightly shallower than that derived by Novak et al. (2017) (note that Novak et al. 2017 derive their slope from other samples of radio data rather than their own VLA-COSMOS data; at the faint end, data in their fit is drawn from Condon et al. 2002).As noted by Leslie et al. (2020) and highlighted in Figure 8, the impact of different  radio −SFR calibrations is significant.This can be seen most clearly in the differences between our SFRD predictions and those derived by Novak et al. (2017).Given the consistency with our luminosity function measurements out to  ∼ 3, this discrepancy appears to stem from their different, redshift-dependent,  IR -based  radio −SFR conversion.
Our results highlight significant differences between SFRD measurements derived from UV/optical/IR data and those derived from FIR/radio data.As described above, these differences likely stem from a number of sources.Incomplete samples and uncertainties in dust corrections affect samples selected at shorter wavelengths, and may drive some of the differences between the SFRD estimated by Driver et al. (2018) and other estimates.Differences in the adopted SFR calibrations and values assumed for the faint-end slope of luminosity functions affect all estimates, and are most clearly seen from the different SFRD measurements derived using similar data (e.g.discrepancies between the 1.4 GHz-derived SFRD measurements constructed by Karim et al. 2011Novak et al. 2017, Leslie et al. 2020, and Enia et al. 2022).In this work, we fix the faintend slope of the radio luminosity function to the value derived at  = 0.03 − 0.30.Its true value is unconstrained by our data at higher redshifts, and evolution would lead to systematic errors in our SFRD estimates.For example, Yüksel et al. (2008) noted that a steeper faint-end slope at high redshift could help reconcile SFRD estimates made by integrating UV luminosity functions with those made using gamma ray bursts.

Comparison to models of galaxy formation
We show predictions for the SFRD from various models of galaxy formation in Figure 12.We note here that a proper comparison would involve making predictions for the multi-wavelength emission (including radio continuum) of the simulated sources, folding in source detection and classification based on the mock SEDs and repeating the analysis on analogously-selected galaxy samples.This is clearly beyond the scope of this paper; instead we present a brief comparison with some initial thoughts here.FIREbox (Feldmann et al. 2022) evolves a small cosmological volume (22.1 cMpc) 3 down to  = 0 using the models initially designed for zoom-in galaxies within the 'Feedback in Realistic Environments' (FIRE) project (Hopkins et al. 2014(Hopkins et al. , 2018(Hopkins et al. , 2023)).FIREbox represents a simulation with particularly high dynamic range, given its low baryonic particle mass ( baryon = 6 × 10 4 M ⊙ ) and medium box size.The total volume-averaged star formation rate density lies above our estimate at all redshifts, with particular deviations at  < 1 and  > 3.This could be due to a number of factors, including galaxy selection effects.When a stellar mass selection of  ★ > 10 9.3 M ⊙ (approximately corresponding to  > 0.03  ★ ) is applied to the simulated galaxies, the predicted SFRD changes significantly, showing better agreement with our data at  ∼ 1 − 2 and a more rapid decrease at  ≳ 2. Importantly, the difference between FIREbox predictions and our data at high- is largely associated with low mass, lower SFR objects that fall below the detection limit of current radio surveys.Deeper radio surveys, or potentially radio-stacking approaches (e.g.Leslie et al. 2020), may bring the data into better agreement with the simulation.At low redshifts, FIREbox underestimates the fraction of massive, quenched galaxies (Feldmann et al. 2022, see also Parsotan et al. 2021); this is likely responsible for the overestimation of the SFRD relative to our estimate at  ≲ 1.The inclusion of AGN feedback in the simulation may help suppress star formation and alleviate this (Cochrane et al. in preparation).
We also compare our estimates to data from the cosmological hydrodynamical simulation Simba (Davé et al. 2019).SFRD estimates from Simba include all the star formation in the box, with no cut on galaxy stellar mass.As noted by Davé et al. (2019), the Simba-predicted SFRD peaks slightly earlier than the compilation of Madau & Dickinson (2014); since our LOFAR estimate of the position of the peak is similar to that of Madau & Dickinson (2014), the Simba SFRD also peaks earlier than our LOFAR estimate.At all redshifts, Simba predicts a lower SFRD than measured from LOFAR, with the most substantial discrepancies (of more than a factor of 2) seen at  ∼ 2. This is consistent with the star-forming main sequence of Simba galaxies displaying a lower normalisation than is observed at this epoch.

Further work
In this work, we have provided important constraints on the cosmic history of star formation based on a statistical study of the deepest available low-frequency radio source counts.However, there are some areas in which our method might be improved in future work.We have adopted a simplified model of the separation of radio sources into those produced by AGN and those produced by stars; in reality, there will be sources where both mechanisms contribute to the radio emission.There will be a jet contribution to the radio emission in some of the radio-quiet AGN included in the sample in this paper, and a star formation contribution to the radio emission in some of the radio-loud AGN excluded from this work.It has long been known that synchrotron radio jets ejected by the AGN can induce star formation as they propagate outwards from their host galaxy nuclei into the galactic and intergalactic medium (e.g.Rees 1989;Gaibler et al. 2012).A prominent low-redshift example is 3C277.3/ComaA (Miley et al. 1981;Capetti et al. 2022).There are also indications that this mechanism could be important at high redshift; in both the spiderweb proto-cluster at  = 2.2 and 4C41.17 at  = 3.8, alignments seen between the radio, optical, CO and X-ray emission have been interpreted as star formation being induced by the radio jets (Bicknell et al. 2000;Miley & De Breuck 2008; see also Duncan et al. 2023).
LOFAR's unique combination of sensitivity and high resolution at low frequencies equips it well to detect and map radio-loud galaxies out to the highest redshifts.Imaging with the international baselines will help to distinguish radio jets from star formation morphologically (e.g.Morabito et al. 2022a).In addition, the new William Herschel Telescope Enhanced Area Velocity Explorer (WEAVE; Jin et al. 2022), a multi-object fiber-fed spectrograph that has just seen first light, will target all radio-detected sources within the LOFAR Deep Fields (WEAVE-LOFAR; Smith et al. 2016).This will provide a vastly larger number of spectroscopic redshifts for the radio sources, additionally enabling better source classifications.

CONCLUSIONS
In this paper, we have used data from the pioneering wide and deep LOFAR Deep Fields to study the cosmic star formation history in a dust-independent manner.The three fields studied, Elais-N1, Boötes and the Lockman Hole, all benefit from extensive UV-FIR coverage, enabling the reliable exclusion of AGN-dominated radio sources from our analysis.We derive 150 MHz luminosity functions for samples of galaxies with radio emission dominated by star formation, from  ∼ 0 to  ∼ 4. Our main conclusions are summarised here: • Out to  ∼ 3, our 150 MHz luminosity functions are in good agreement with the scaled 1.4 GHz luminosity functions derived by Novak et al. ( 2017) using VLA-COSMOS data (assuming a spectral index  = −0.7).Given the larger area spanned by the LOFAR Deep Fields (∼ 25 deg 2 , compared to the ∼ 2 deg 2 VLA-COSMOS survey), and the use of three fields to overcome cosmic variance, we can constrain radio luminosity functions to roughly an order of magnitude brighter luminosities, while reaching similar luminosities at the faint end.• Our derived 0 <  < 0.3 150 MHz luminosity function is wellfitted by a parametrisation of the form: with log 10 ( ★ / Mpc −3 dex −1 ) = −2.46 ± 0.01,  = 0.49 ± 0.01,  = 1.12 ± 0.01 and log 10 ( ★ /W Hz −1 ) = 22.40 +0.02 −0.03 .• Using the values of  and  derived to our low redshift data, we fit our higher redshift radio luminosity functions using the same parametrisation, to constrain the evolution of  ★ and  ★ .• We show that transforming a radio luminosity function to a star formation rate function is complicated by the scatter in the  radio − SFR relation.Star formation rate functions derived using this conversion tend to lie above those derived using SFRs obtained from SED fitting at SFR ≳ SFR ★ , with deviations of up to an order of magnitude at the highest SFRs.Using a simple model, we show that higher values of scatter in the  radio − SFR cause a gentler fall-off of the exponential of the measured radio luminosity function.This effect is most important where the luminosity (or star-formation rate) function is steepest, above its break.This can lead to an overestimation of the cosmic star formation rate density, which is generally derived by integrating the measured luminosity function.The magnitude of the correction depends on the the form of the luminosity function, the scatter in the  radio − SFR relation, and the details of the sample selection (i.e.whether sources that are particularly radio-loud for their SFR are excluded from the sample).By comparing the difference in the inferred SFRFs using the two methods to our model of the bias, we constrain the scatter in the  radio − SFR relation of star-forming galaxies to be ∼ 0.3 dex.Encouragingly, this value is in line with recent work that constrains the scatter using an independent method (Smith et al. 2021).We derive an appropriate correction factor to apply to the SFRD of ∼ 0.96 ± 0.04.• We constrain the cosmic star formation rate density from  ∼ 0 to  ∼ 4, by integrating our  150 MHz luminosity functions, in combination with a self-consistently derived  150 MHz − SFR relation, correcting for its scatter.Since the SFRD is constructed using radio-selected samples, our measurements are robust to the effects of dust.Our derived SFRD lies between previous compilations at all redshifts studied.Our measurements are in good agreement with those previously derived using smaller 1.4 GHzselected samples (e.g.Leslie et al. 2020;Enia et al. 2022) and from FIR-based studies (e.g.Gruppioni et al. 2013;Dunlop et al. 2017).Our derived SFRD is well-fitted by a model of the form SFRD() =  10 (− 0 ) +10 (− 0 ) , with  = −0.89+0.08  −0.07 ,  = 0.22 ± 0.04, log 10  = −0.76± 0.05 and  0 = 1.22 ± 0.15.
Prospects for future census studies of radio-selected star-forming galaxies are bright.The LOFAR Deep Fields survey continues to observe all three fields.EN1 has already been observed for 500 hr, with imaging reaching ∼ 12 Jy/beam.By mid-2023 (following LOFAR Cycles 18 and 19), we expect to reach 16 Jy/beam in Boötes with 312 hr of data, and 13 Jy/beam in the Lockman Hole with 352 hr of data.We are also observing in the NEP, where we expect to reach 13 Jy/beam in 400 hr.These deeper LoTSS radio data, alongside spectra from WEAVE-LOFAR, will build up large, dust-independent samples of star-forming galaxies for further study, at fainter star formation rates than previously possible.This will enable not only the characterisation of the global cosmic SFRD but also the investigation of the drivers of star formation and quenching in sub-populations over cosmic time.3.

APPENDIX A: RADIO COMPLETENESS CORRECTIONS
In Figure 2, we show radio completeness as a function of 150 MHz flux density for each LOFAR Deep Field.We tabulate these values in Table A1.

APPENDIX B: PHOTOMETRIC UNCERTAINTY CORRECTIONS
In Section 3.3, we describe corrections derived to account for uncertainties in photometric redshifts.These corrections are shown for each field in Figure 3 and tabulated in Table B1.

APPENDIX C: SFRD ESTIMATES FOR EACH FIELD
Estimates for the star formation rate density for each individual LO-FAR Deep Field are presented in Table C1.

APPENDIX E: THE IMPACT OF INTEGRATION LIMITS ON THE DERIVED SFRD
We have tested the impact of changing the range of radio luminosities over which we integrate to obtain the SFRD.In Figure E1, we plot the Boötes-derived SFRD for different choices of lower (left) and upper (right) luminosity limits.We find that our derived SFRD is robust to changes in the lower limit between 0.003  ★ and 0.1  ★ (we adopt 0.03  ★ in this work) and to changes in the upper limit between 10 27.5 W Hz −1 and 10 28.5 W Hz −1 (we adopt 10 28 W Hz −1 in this work).2021), is shown in blue, purple and red for Elais-N1, Boötes and the Lockman Hole, respectively.In orange, we show the SFRD that would be derived without correcting for uncertainties in the photometric redshifts (see Figure 3 and Table B1) or scatter in the  150 MHz − SFR relation (see Section 4.2).Corrections for photometric redshift uncertainties have minor effects on the majority of redshift bins, but are important at  ∼ 1 − 2 for Elais-N1 and the Lockman Hole, which lack −band data.Figure E1: The impact of the choice of lower (left) and upper (right) integration limit on the derived SFRD, shown here for Boötes.Small artificial x-axis offsets are added to display small differences most clearly.The impact of changing either the lower or the upper integration limit is small compared to the reported uncertainties.

Figure 1 :
Figure 1: The distribution of 150 MHz radio luminosity, star formation rate and stellar mass, for the main sample of galaxies studied in this paper (all galaxies within the three fields for which radio continuum emission at 150 MHz is dominated by star formation; see Section 2.5 for details of the SFG/AGN separation and sample selection).The redshifts plotted are  BEST -the spectroscopic redshift where this is available, and the median redshift of the preferred photometric redshift solution,  1,MEDIAN , where it is not.The stellar mass and SFR values plotted are the consensus estimates derived from the combination of four different SED fitting codes (see Section 2.6 and Best et al. 2023).

Figure 2 :
Figure2: Completeness as a function of 150 MHz flux density, for each of the three LoTSS Deep Fields.As described in Section 3.2, the completeness curves were calculated using source extraction of mock input sources from the radio images, using the same PyBDSF parameters as were used for the extraction of real sources.A realistic source size distribution was applied, based on the observed size distribution of star-forming galaxies with flux densities in the range 1 − 5 mJy, where completeness is high.As expected based on total observing times, Elais-N1 has the deepest radio coverage, followed by the Lockman Hole, and then Boötes.Note that due to a very small number of injected mock sources overlapping with real sources and other injected sources, the completeness never quite reaches 100 per cent.

Figure 4 :
Figure 4: Characterisation of the form of the low redshift 150 MHz luminosity function for star-forming galaxies.Left: the luminosity function of star-forming galaxies in the redshift range  = 0.03 − 0.30, constructed using sources in the three LOFAR Deep Fields that show no evidence of AGN-driven radio emission.The dark blue points show our estimate including the radio completeness corrections, with error bars calculated using 1000 bootstrapped samples.Data points without the completeness correction are shown by open squares (black).Our data are fitted with a Saunders et al. (1990) parametrisation; the pale blue shaded region shows the 16 th − 84 th percentile of the posterior distribution of one example of the binning scheme.For comparison, the local ( < 0.3) luminosity function of 150 MHz-selected star-forming galaxies from LoTSS DR1, derived by Sabater et al. (2019), is shown in green.The local ( < 0.3) luminosity function of 1.4 GHz-detected star-forming galaxies (scaled using a radio spectral index  = −0.7),derived by Best & Heckman (2012), is shown in orange.The  = 0.03 − 0.30 luminosity function for radio-selected AGN, derived by Kondapally et al. (2022), is shown in pink.Right: the posterior distributions of the four fitted parameters, which are well-constrained.The fitting was repeated for 18 combinations of luminosity bin size and position, and the results were then averaged.The fitted values depend little on the choice of luminosity binning.Shown in the right hand panel is just one example of the posteriors from one choice of binning; we use the average fitted values of  = 0.49 and  = 1.12 in our fits at other redshifts.

Figure 5 :
Figure5: The 150 MHz luminosity function for galaxies with radio emission dominated by star formation, for the three LOFAR Deep Fields: Elais-N1 (mid blue), Boötes (purple) and the Lockman Hole (red), as well as for the combined sample (navy).Note the change in x-axis scale for the different redshifts.Corrections for uncertainties in photometric redshifts have been applied as a fixed scaling for a given redshift for each field, as described in Section 3.3.Corrections for radio completeness have also been made; data are only plotted for bins where this correction is < 0.5 dex.The luminosity functions show excellent consistency between the three fields and also (except at the highest redshifts) good agreement with the estimates ofNovak et al. (2017), shown in green, which have been scaled from 3 GHz using a radio spectral index  = −0.7.

Figure 6 :
Figure 6: Radio luminosity functions derived using all three LOFAR Deep Fields and presented in Figure 5, plotted on a single figure to illustrate the strong redshift evolution.

Figure 7 :
Figure7: The evolution of best-fitting luminosity function parameters, for the three LoTSS Deep Fields individually and combined.Values of  ★ have been corrected for uncertainties in photometric redshifts, as described in Section 3.3.To guide more physical intuition, we show in the left-hand panel the characteristic star formation rate corresponding to the characteristic luminosity, derived using the  150 MHz − SFR relation given in Equation6.The derived parameter values are generally in good agreement across the fields.The three estimates of  ★ and  ★ diverge at the highest redshift ( ∼ 4), where  ★ falls well below the sensitivity limits of the radio data.
argued that the stellar mass dependence of the relation can introduce substantial systematic errors (of order 0.5 dex) on SFRs derived from  150MHz alone, particularly in cases where the sample for which SFRs are derived has a different stellar mass distribution to the sample from which the relation was derived.They noted that these offsets are potentially larger than the intrinsic scatter in the  150 MHz − SFR relation.

Figure 8 :
Figure8: The SFR function derived from the 150 MHz-selected sample in all three LOFAR Deep Fields, constructed using different SFR estimates.The SFR function derived using the consensus SFR estimates derived byBest et al. (2023) is shown in dashed blue.The solid red line shows that derived using the  150 MHz − SFR conversion derived bySmith et al. (2021).The dashed yellow line shows that derived using the mass-dependent  150 MHz − SFR conversion derived bySmith et al. (2021).Despite the  150 MHz − SFR conversion being derived from the same LOFAR Deep Fields data, the derived SFR functions differ above SFR ★ , in some cases by an order of magnitude.

Figure 9 :Figure 10 :
Figure 9: Constraining the scatter in the  150 MHz −SFR relation.Left: simulated deviation of the star formation rate function, as measured from  150 MHz , from the input star formation rate function (black), for different values of the scatter in the  150 MHz − SFR relation.Larger scatter in the relation causes the observed number density of galaxies to increase at the bright end, causing a gentler fall-off of the exponential.This leads to an overestimation of the cosmic star formation rate density.Right: the method used to constrain the scatter in the  150 MHz − SFR relation.Grey data points show offsets between the star formation rate function constructed using consensus SFR estimates, and that constructed by applying a single  150 MHz − SFR relation to the 150 MHz luminosity function.Coloured lines show the modelled offsets for different values of scatter in the relation.From this, we calibrate the scatter to be ∼ 0.2 − 0.4 dex.Integrating these coloured lines, and comparing to the integrated black line (no scatter), we obtain scaling factors to correct values of the cosmic star formation rate density derived using a single  150 MHz − SFR conversion.These are: CORR SFRD = [1.0,0.96, 0.93], for scatter Σ = [0.2,0.3, 0.4] dex.
, we show the impact of applying corrections for uncertainties in photometric redshifts and the scatter in the  150 MHz −SFR relation on the derived SFRD for each of the three fields.Corrections for photometric redshift uncertainties lead to small changes in the estimated SFRD for the majority of redshift bins, but are important at  ∼ 1 − 2 for Elais-N1 and the Lockman Hole, the two fields lacking -band data.The correction for the scatter in the  150 MHz − SFR relation serves to move the whole derived relation to slightly lower SFRD values.
FigureD1: The cosmic star formation rate density, derived by integrating the radio luminosity function using the  150 MHz − SFR relation ofSmith et al. (2021), is shown in blue, purple and red for Elais-N1, Boötes and the Lockman Hole, respectively.In orange, we show the SFRD that would be derived without correcting for uncertainties in the photometric redshifts (see Figure3and TableB1) or scatter in the  150 MHz − SFR relation (see Section 4.2).Corrections for photometric redshift uncertainties have minor effects on the majority of redshift bins, but are important at  ∼ 1 − 2 for Elais-N1 and the Lockman Hole, which lack −band data.

Table 2 :
Saunders et al. (1990)cosmic star formation history 11Redshift  sources log 10 ( ★,150MHz /W Hz −1 ) log 10 ( ★ /Mpc −3 log L −1 ) SFRD/M ⊙ yr −1 Mpc −3 Parameters derived from fits of theSaunders et al. (1990)luminosity function to the data shown in Figure5.For all redshift bins apart from  = 1.3 − 1.6, the luminosity function is derived using data from all three LOFAR Deep Fields.At  = 1.3 − 2.0, our best estimates of the luminosity function are from the Boötes field alone, due to the larger photometric uncertainties in Elais-N1 and the Lockman Hole (see Section 3.3).★represents our best estimate of the characteristic number density, after correcting for photometric redshift uncertainties, as described in Section 3.3.We also provide our estimate of the cosmic star formation rate density integrated down to 0.03  ★ at each redshift, calculated using a  150 MHz − SFR calibration derived from the same dataset, in combination with a correction for the scatter in this relation (see Section 4.2).SRD estimates for each individual field are provided in TableC1.
Gürkan et al. (2018)l intrinsic scatter on the relation (∼ 0.3 dex) at SFR = 10 − 100 M ⊙ yr −1 .Adding an additional parameter such as stellar mass can decrease the scatter on this relation.Using shallower LOFAR data in the Herschel-ATLAS NGP field, and focusing on star-forming galaxies at  < 0.4,Gürkan et al. (2018)found evidence of a dependence of the radio luminosity on both SFR (the primary driver) and stellar mass (Davé et al. 2019)022erived cosmic star formation rate density presented in Figure11(navy stars, with fitted functional form in shaded navy), with predictions from hydrodynamic simulations overlaid.The total SFRD predicted by FIREBox (dark red;Feldmann et al. 2022) exceeds our estimate, with the greatest deviations at  < 1 and  > 3. When the FIREbox SFRD is calculated using only galaxies with  ★ > 10 9.3 M ⊙ (pale red; this limit corresponds approximately to 0.03  ★ , which we integrate down to in this work), agreement is better at  ∼ 1 − 2. At lower redshifts, the FIREbox SFRD estimate is up to a factor of a few higher than our data, likely due to the lack of AGN feedback in the simulations.Simba(Davé et al. 2019)under-predicts the SFRD at all redshifts, with the most substantial discrepancies at  ∼ 2.

Table A1 :
Kondapally et al. (2022)ux densities in the range 0.11 − 40 mJy, for each field studied.The source size distribution used was defined using the size distribution of star-forming galaxies with flux densities in the range 1 − 5 mJy.A similar table for a source size distribution appropriate for AGN is presented byKondapally et al. (2022).

Table B1 :
Photometric correction factors, for each field, as described in Section 3.3 and shown in Figure