Exploring the faintest end of mid-infrared luminosity functions up to $z\simeq 5$ with the JWST CEERS survey

Mid-infrared (MIR) light from galaxies is sensitive to dust-obscured star-formation activities because it traces the characteristic emission of dust heated by young, massive stars. By constructing the MIR luminosity functions (LFs), we are able to quantify the overall dusty star formation history and the evolution of galaxies over cosmic time. In this work, we report the first rest-frame MIR LFs at 7.7, 10, 12.8, 15, 18, and 21 $\mu$m as well as the total IR LF from the James Webb Space Telescope (JWST) Cosmic Evolution Early Release Science (CEERS) survey. We identify 506 galaxies at $z=0-5.1$ in the CEERS survey that also have optical photometry from the Hubble Space Telescope. With the unprecedented sensitivity of the JWST, we probe the faintest end of the LFs at $z=0-1$ down to $L^* \sim 10^7 L_\odot$, $\sim 2$ orders of magnitude fainter than those from the previous generation of IR space telescopes. Our findings connect well with and continue the faint end of the MIR LFs from the deepest observations in past works. As a proxy of star formation history, we present the MIR-based luminosity density up to $z\simeq4.0$, marking the first probe of the early Universe by JWST MIRI.


INTRODUCTION
The galaxy luminosity function (LF) is a statistic that profiles the number density of luminosities of galaxies in a specific volume, redshift range, and population.Serving as one of the fundamental quantities in extra-galactic astronomy, LF provides crucial and direct information for studying the evolution and properties of galaxies.The star-forming (SF) activities of galaxies play a significant role in shaping and evolving LF, as the star formation rate (SFR) and the ultra-violet (UV) / infrared (IR) luminosity are linked by well-known empirical relations (Kennicutt 1998;Madau et al. 1998).
While UV light provides direct evidence of SF, IR-based LF is vital for a comprehensive understanding of galaxy evolution because it can trace the dust-obscured star formation that is opaque to UV.This makes the IR-based LF particularly crucial for studying the high- galaxy population, where the dust-obscured star formation is expected to be more significant (e.g., Hopkins et al. 2001).Specifically, we focus on the mid-infrared (MIR) wavelength due to its unique ability to detect polycyclic aromatic hydrocarbon (PAH) emission from 3 to 20 m.The emission, along with its characteristic features at approximately 6.2, 7.7, and 8.6 m, is mainly due to the heated dust by young stars.The MIR/PAH emissions dominate the spectral energy distribution (SED) of SF galaxies and have already been proven as a good indicator of total IR (TIR) luminosity  TIR (which is estimated by integrating the galaxy SED from 8 to 1000 m), e.g., Caputi et al. 2007;Goto et al. 2011 for 8 m, Pérez-González et al. 2005 for 12 m, and SFR (Wu et al. 2005;Calzetti et al. 2005Calzetti et al. , 2007)).
Studies on the IR LF can be dated back to the first-generation IR space telescopes IRAS (Neugebauer et al. 1984) and ISO (Kessler et al. 1996), and continue with AKARI (Murakami et al. 2007), Spitzer (Werner et al. 2004), and Herschel (Pilbratt et al. 2010).These works (e.g., Saunders et al. 1990;Rowan-Robinson et al. 1997;Babbedge et al. 2006;Caputi et al. 2007;Goto et al. 2010Goto et al. , 2011Goto et al. , 2019) ) have proven that IR LF is a powerful tool to unravel the evolution in dusty galaxies and have greatly improved our understanding of the cosmic star formation history (CSFH).The evolution of IR LFs has been extended to the whole  = 0 − 4 range since Herschel era (Gruppioni et al. 2013;Magnelli et al. 2013), and is recently pushed to  ≃ 6 (Gruppioni et al. 2020) with ALMA.With the launch of the James Webb Space Telescope (JWST, Gardner et al. 2006;Kalirai 2018), we are now able to expand our understanding of IR sources to an unprecedented degree.For instance, Donnan et al. (2023) and Harikane et al. (2023) suggest UV LFs at  = 8 − 15 and  = 8 − 17, respectively, with JWST Near-Infrared Camera (NIRCam).On the other hand, source counts studied at 7 − 21 m MIR wavelengths covered by the Mid-Infrared Instrument (MIRI, Rieke et al. 2015) have demonstrated exceptional sensitivity of MIRI to below Jy level (e.g., Ling et al. 2022;Wu et al. 2023).Kim et al. (2024) interpret the cosmic star-formation history (CSFH) and black hole accretion history (BHAH) by applying the backward evolution of local LFs (Gruppioni et al. 2011) to their source counts (Ling et al. 2022;Wu et al. 2023).
This study seeks to further the efforts from the preceding works by constructing the first IR LF from JWST MIRI observations.The precise MIR photometry, combined with NIR and optical data from Hubble Space Telescope (HST) and large ground-based telescopes, can provide a complete picture of SED for faint and distant galaxies.We target to illustrate the MIR LF as well as its evolution to the faintest and oldest end that has not been entirely explored previously.This paper is organised as follows: We introduce the JWST MIR data, completeness, and the multi-wavelength merged catalogue we produced in §2.The procedure to derive the LF is explained in §3.In §4, we present our LFs and discuss their evolution.The conclusion is given in §5.We adopt the Planck18 cosmology (Planck Collaboration et al. 2020), i.e., Λ cold dark matter cosmology with (Ω  , Ω Λ , Ω  , ℎ) = (0.310, 0.689, 0.0490, 0.677) throughout the paper.

MIRI observations
We utilise images from JWST Cosmic Evolution Early Release Science (CEERS; Finkelstein et al. 2017) survey to construct a MIR source catalogue.The CEERS survey is one of the Early Release Science programs of JWST, providing the first observations from NIRCam, Near-Infrared Spectrograph (NIRSpec), and MIRI in the Extended Groth Strip (EGS) legacy field.This work focuses on the MIRI pointings in the observation that were observed using 6 continuous broad-band filters (F770W, F1000W, F1280W, F1500W, F1800W and F2100W), covering a wavelength range from 7.7 m to 21.0 m.Unlike the proposed observation (Finkelstein et al. 2017), only two pointings (observation ID: o001_t021 and o002_t022, hereafter "the first two") have been observed using all 6 filters.F770W and F2100W observations are missing in the other two (o012_t026 and o015_t028, hereafter "the last two") pointings, leaving 4 available filters for the pointings.To obtain a larger sample, we use all the four pointings in this work, regardless of the waveband coverage in 7.7 and 21 m.The total sky coverage for the 4 pointings is 31380.56arcsec 2 (∼ 8.7 arcmin 2 ), and 15694.27arcsec 2 for F770W and F2100W which have only 2 pointings.
We obtain level-3 image product of the pointings from the Mikulski Archive for Space Telescopes (MAST) and conduct source extraction and photometry following Wu et al. (2023).Two photometry software, Photutils (Bradley et al. 2022) and Source-Extractor (Bertin & Arnouts 1996) are used in the procedure.Wu et al. (2023) found that applying two software for background estimation and photometry separately is more effective.Please refer to Wu et al. (2023) for the technical details of the photometry, including the effectiveness and the result compared to the public source catalogue from MAST.In order to keep a sufficient sample size, we do not select sources by a specific band.Instead, any sources with at least one detection in MIRI bands will be included in the catalogue.A total of 1210 sources are identified through this selection.

Completeness
Completeness correction is essential to source statistics.The completeness assesses the reliability of the source within an image, which is affected by the sensitivity and exposure time of the image.Wu et al. (2023) has presented an effective approach to probe the completeness of JWST images.Here, we measure the completeness function for all 4 CEERS pointings by performing Monte Carlo (MC) simulations with the same method as in Wu et al. (2023) and Takagi et al. (2012).The method is briefly described as follows.First, we randomly implement 20 artificial MIRI PSF sources with the same flux into the target image for each flux bin.The width of the flux bin is Δ log (   /Jy) 0.1 dex.The same photometric analysis is subsequently applied to these mock sources to determine their recovery rate as the final completeness.The simulation was repeated 500 times to reduce the uncertainty.We build the completeness as a function of flux for each pointing and filter.Figure 1 shows completeness functions of the 4 CEERS pointings.The 80% completeness limits for each pointing and MIRI band are concluded in Table 1.These limits are taken by a linear interpolation between our data points in Figure 1.We find that the completeness and the 80% limit of the last two pointings are generally higher than the first two.This is due to the shorter exposure time in the last two pointings.The average exposure time of the last two pointings is 1243 seconds (F1000W and F1800W) and 932 seconds (F1280W and F1500W), while for the first two is 1673 seconds (except for F2100W, which has a much longer exposure time of 4811 seconds).In Figure 2, we show the cumulative areal coverage of images of the sum of the 4 pointings as a function of 1 sigma flux error of each pixel, obtained from the fits images.The vertically rising line in Figure 2 indicates an identical pixel error for most individual pointings, thus ensuring the uniform exposure time.Meanwhile, the stepwise increase illustrates the difference in exposure time between the first two and last two pointings.For F1280W and F1500W bands, we notice a more gradual rise because of variations in exposures within the last two pointings.We should note that no additional correction in completeness is applied to account for these variations because the random distribution of artificial sources in the MC simulation already includes and averages for such effects.
We inspect false/spurious detections in the images by evaluating their reliability, following Wu et al. (2023).The reliability is assessed as (N pos -N neg ) / N pos , where N pos is the number of sources detected in the original images and N neg is the number in the negative images.The same photometry procedures are applied to the negative images (i.e., the original MIRI pointing multiplied by -1) to extract sources.We show the reliability for each pointing and band beyond their 80% completeness levels (refer to Table 1) in Table 2.All the pointing have a high level of reliability > 90%, except for F1280W and F1500W band in pointing o015_t028.This is expected because of the nonuniform and shorter exposures in these pointings, as suggested above.The completeness functions and reliability are accounted for in the derivation of LF to correct possible incompleteness and spurious detection of our sources.We elaborate on the correction procedure in §3.3.

Multi-wavelength merged catalogue
A comprehensive photometric catalogue covering multiple wavelengths is necessary to accurately model the SEDs of our MIR sources.This is achieved by matching our source catalogue with the CANDELS-EGS Multi-wavelength catalogue (Stefanon et al. 2017) EGS field.The CANDELS-EGS sources are selected by detections in HST WFC3 F160W with a depth of 26.62 AB (90% completeness).An accompanying redshift catalogue (Kodra et al. 2023) of CANDELS-EGS sources is also combined for analysis in §3.1.We apply a matching radius of 0.5 arcsec, considering the average size of our sources is ∼ 0.25 arcsec.94% of identified pairs have a separation of less than 0.25 arcsec.
The compiled catalogue contains 573 sources with 15 + 6 band photometry from the CANDELS-EGS catalogue and our JWST MIR source catalogue.Table 3 summarises the depths in the catalogue.We should emphasise that the first two MIRI pointings do not include NIRCam photometry due to the lack of NIRCam observations in the pointings.Instead, we take the data from the CANDELS-EGS catalogue for NIR wavelengths, which are obtained by CFHT WIRCam, HST WFC3, and Spitzer IRAC (refer to Table 3).
637 (52% of the 1210 sources) MIR sources have not been identified in the CANDELS-EGS catalogue.We exclude these sources from the compiled catalogue because their absence in optical/NIR bands could result in poor SED fit.As shown in Figure 3, the majority of these sources (which come from F1280W, F1500W, and F2100W filters) mostly lie outside the 80% completeness limit of the filters.For further check, we have eye-balled those EGS-undetected sources and found that most of them in F2100W are spurious among the noisy background, justifying the exclusion.This is also noted by the CEERS team (Yang et al. 2023b), which shows the measured depth for F2100W is shallower than the Exposure Time Calculator prediction by ∼ 0.9 magnitude.In F1500W, about half of the EGS-undetected sources are close to the edge of the field, where the dithering of CANDELS-EGS and other MIRI filters just missed.Here we forewarn readers for the treatment of these spurious sources, especially in F2100W filter.

ANALYSIS
There are several procedures involved in obtaining LF.We describe them step by step in the following subsections.In §3.1, we introduce the SED fitting results of our galaxies and evaluate their performance.In §3.2, we present the conversion from the observed flux in the SED to the rest-frame luminosity with -correction.After the conversion, we show how we construct the rest-frame LFs in §3.3, where the correction to the completeness of our galaxy sample is also addressed.
Figure 3.The magnitude histogram of all MIR source detection in each filter.Sources that are also identified in the CANDELS-EGS catalogue are marked as "match" (blue), otherwise, they are marked as "unmatched" (orange).The black dashed (dash-dotted) line indicates the 80% completeness of the filter for pointing o002_t022 (o015_t028), converted to AB magnitude.Note that sources may be counted more than once here, as they are detected in multiple filters.model SEDs and physical properties of galaxies with observations across from far-UV to radio spectrum.cigale offers users various modules and parameters to optimise the fit.For our study, we mainly follow the fitting configuration described in Yang et al. (2023a), which uses the same sample as us and analyses SEDs to investigate the Active Galactic Nuclei (AGN) population at high-.The modules used for the fitting are presented in Table 4. Unlisted parameters remain the default value from cigale.

Module
The main difference between our parameters and Yang et al. ( 2023a) is that we adopt dustatt_modified_CF00 module for dust attenuation, which reduces the systematically lower photo- from cigale (Figure 2 in Yang et al. 2023a).We also utilise a denser redshift grid that ranges from 0.01 to 6.0 with 600 linearly spaced steps, yielding a step of 0.01 in redshift.This is preferred as the luminosity scales with the square of the distance, and a denser grid can reduce potential biases.
The fitted galaxy SEDs are separated into SF and AGN types based on cigale parameter frac AGN .frac AGN is the AGN contribution to IR luminosity ratio in Table 4, defined by  AGN  AGN + galaxy within 3 − 30 m.Galaxies with frac AGN exceeding 20% are categorised as AGN hosts (following Wang et al. 2020, in which they used catalogues from Kim et al. 2021 andHo et al. 2021), otherwise as an SF galaxy.133 (26%) of our sample are identified as AGN host galaxies.Examples of best-fit SED results from cigale for each redshift bin and SED type are presented in Figures A1-A4.Additionally, the fitted far-IR SED is included in these figures for reference.
It is crucial to note that our far-IR SED fitting is dependent on observations at shorter MIR wavelengths due to the absence of far-IR detection; most galaxies are too faint to be observed with Herschel.Despite the potential impact this may have in deriving  TIR , we would like to point out that  TIR is constrained to a certain degree by  MIR1 or  PAH , especially for the SF galaxy, which is the majority of our sample.The empirical relation between  MIR and  TIR has been reported in the literature (e.g., Caputi et al. 2007;Houck et al. 2007;Goto et al. 2011;Lin et al. 2024), where  MIR has been demonstrated as a reliable proxy for  TIR given that both luminosities are primarily attributed to SF activity.On this consideration, the derivation of  TIR should remain valid and in line with the current understanding.
We remind readers, however, that the assumption can only hold if both  MIR and  TIR are attributed to SF activity.If the AGN dominates the mid-IR emission, the contamination from AGN in the mid-IR can be as high as 80 − 100% (e.g., the right panel of Figure A4).Therefore, in the presence of an AGN, the AGN component must be subtracted from the total SED fit, as we demonstrate in §4.2.
To examine the overall quality of the fit, we plot the distribution of reduced  2 of each source's fit in Figure 4.The median of the distribution is 1.39.A criterion of reduced  2 < 5 is then set to exclude 16 poor fits.We also present an analysis similar to Yang et al. (2023a) for redshift.In Figure 5, we compare the photo- estimated by cigale with the redshift catalogue (Kodra et al. 2023) which provides spectroscopic redshift (spec-) measurements for 107 sources in our sample.The assessment of photo-, i.e., the normalised median absolute deviation  NMAD (defined as 1.48 × median{|Δ − median(Δ)|/(1 +  spec )}) and spec- outlier fraction  (defined as Δ/(1 +  spec ) > 15%), can be derived by comparing these 107 galaxies with spec- (blue open circles in Figure 5) to their photo-, Δ = | photo −  spec |.In result, we obtain  NMAD = 0.029 and  = 0.00%.These values are comparable to those reported in Kodra et al. (2023) ( NMAD = 0.0227;  = 6.7%) and Yang et al. (2023a) ( NMAD = 0.031;  = 0.00%).We note, the 0% outlier fraction (as in Yang et al. 2023a) simply reflects the fact that all the spec- samples have a photo- (from cigale) with < 15% deviation.For consistency, we excluded 51 (9.16%) galaxies that are outside the 15% redshift uncertainty range from our sample, i.e., the photo- outliers shown in the upper panel of Figure 5.In summary, we apply below criteria for the reduction: • Reduced  2 of SED fitting > 5 • Photo- outliers with > 15% deviation from Kodra et al. (2023) This leaves 506 available photo- and SEDs for the following procedures.To compensate for the removal by the criteria, We multiply the luminosity functions in all the -bins by a factor of 573/506 ≃ 1.13 in the following derivation (see §3.3), with the underlying assumption that their redshift distribution is the same as the rest.

𝐾-correction
Due to the expansion of the Universe, the observed fluxes are redshifted from the rest/emitted-frame.The difference in transmission efficiency between the two filters at different redshifts has to be considered to obtain accurate measurements for luminosity, especially since we are interested in the MIR spectra related to star formation.
We apply the -correction (Oke & Sandage 1968) to correct the effect, with the correction factor  () in terms of luminosity for the AB magnitude system. o and  e are the filter wavelength in the observed and rest/emitted-frame respectively, and  () is the transmission curve for the specific filter.  is the luminosity density per unit wavelength (e.g.,  ⊙ m −1 ).The observed flux   and the rest-frame luminosity   are related by ) where  L is the luminosity distance at redshift .For each filter, we convolve our best-fit galaxy SED in the rest-frame to the transmission curve of the filter to obtain the corrected luminosity.

Luminosity function
The LF describes the number density of galaxies () as a function of their intrinsic brightness  in a specific volume .This volume depends on the redshift distribution of the galaxies, i.e., We should be careful that the bin may include a wide range of variations in evolution.
For a volume-limited sample,  is just the survey volume.However, the survey is often limited by the sensitivity of the telescope.In practice, this results in a flux-limited sample, which can be incomplete in terms of volume.An example is that the survey volume of redshift bin  = [1, 2] cannot be applied to faint galaxies with luminosities below a certain level, because they can only be detectable within  < 1.5 for the filter.
The 1/ max method (Schmidt 1968) is thus introduced to address the volume incompleteness from a flux-limited sample.This method utilised the  max , the maximum redshift where a given object would still be detectable by the telescope, by redshifting the SED of the object toward the flux limit of a specific filter.The final  max is taken as the redshift that reaches the flux limit.If the  max goes out of the redshift bin, we use the upper limit of the redshift bin as  max .Then, the effective survey volume  max for the object can be calculated using equation 4. The exact formula for deriving the LF based on waveband  is given by where Δ log  = 0.5 dex is the width of our luminosity bin, and  , is the correction factor for -th galaxy at a waveband .Specifically, area  is the sum of the sky coverage in every pointing of a waveband  (see §2.1).To obtain the completeness for a given galaxy , we first select the completeness function corresponding to the specific waveband  and pointing position of the galaxy, as shown in Figure 1.Subsequently, the completeness can be computed from the 's flux at waveband  ( , ).A similar approach is applied for the reliability correction, where the correction factor is also set according to the galaxy's pointing and waveband (Table 2).Notably, we derive TIR LF using corrections for F1000W band, because we find the largest sample of galaxies at 10 m (see Figure 6).In addition to the correction for individual galaxies, we add a global compensation factor of ∼ 1.13 (13%) to the LF in each redshift bin for removed galaxies, as indicated in §3.1.
Identifying limiting luminosities to filter out galaxies with insufficient completeness is crucial when interpreting LFs.To achieve this, we have converted the 80% flux completeness limits from Table 1 into limiting luminosities.Two types of SED, SF and AGN galaxies, are assumed in the calculation based on the classification stated in §3.1.The SED of NGC6090 and Sey2 (representative of average Seyfert 2 galaxies) from SWIRE Template Library (Polletta et al. 2007) are utilised as template SEDs for SF and AGN galaxies, respectively.These templates are selected according to Yang et al. (2023a) which suggests that the median SEDs of CEERS MIRI SF and AGN galaxies resemble them respectively.The limiting luminosity is tailored to the redshift, completeness limit, and SED type for each galaxy by the assumed template SED.In practice, we first obtain the 80% completeness limit specific to the LF of the given waveband and the pointing of the galaxy, similar to the correction  procedure above.Next, an appropriate template SED based on the galaxy's type (SF or AGN) is selected and redshifted to the desired redshift.We then integrate for the luminosity (either monochromatic or TIR, depending on the LF) at which the flux of SED at the specific band would be equal to the completeness limit as the limiting luminosity.Galaxies with luminosities that fall below these limiting thresholds are regarded to be unreliable and are thus excluded.
Figure 7 shows the criteria for the TIR LFs (Figure 11) which uses the F1000W band limit.We note the small zigzags in these completeness lines simply reflect the slight difference in flux completeness limits for the galaxies at each CEERS pointing.
The final size of our galaxy samples for monochromatic and TIR LFs after the 80% completeness selection is listed in Table 5.We derived luminosity limits for all redshift bins  =

Monochromatic LFs
We present the monochromatic rest-frame LFs at 7.7, 10, 12.8, 15, 18 and 21 m, based on 6 available MIRI filters from JWST.In §4.1.1 we show the LFs of all galaxy populations, and in §4.1.2we discuss the LFs for SF and AGN host galaxies.To estimate the error, we resample the LF distribution 100 times.This is done by perturbing the redshift (i.e.photo-) of each galaxy with associated probability distributions from cigale and then recalculating the luminosities (monochromatic and TIR), as well as corresponding LFs.Subsequently, we add a Poisson error to each luminosity bin.Bins with less than 3 galaxies are removed to avoid fluctuations from small statistics.
For comparison, we overplot the LFs from previous observations with IR space telescopes.We provide a brief summary of them below.Goto et al. (2019) and Kim et al. (2015) are based on AKARI observations in the 5.4 deg 2 -wide NEP field.Goto et al. (2019) measures the rest-frame 8 and 12 m LFs at 0.35 <  < 2.2 using 18 bands mid-IR photometry (Kim et al. 2012) and 5 optical bands from the Subaru Hyper Suprime-Cam (Goto et al. 2017) 2010) presents rest-frame 8, 12, 15, and 24 m LFs up to  ≃ 2.5.We only plot LFs from the literature whose redshift ranges are similar to the median of our redshift bins.
We first remark on the superior sensitivity of the JWST.In Figure 8, we find the luminosities of galaxies can be probed down to  * ∼ 10 7 − 10 8  ⊙ for the lowest redshift bin  = 0 − 1, while their luminosity limits are at  * ∼ 10 8 − 10 9  ⊙ at higher redshifts.Still, This is about 1 to 2 orders of magnitude fainter than the limiting luminosities found by AKARI and Spitzer at similar redshifts (∼ 10 9 − 10 10  ⊙ ).It is no surprise as such improvement has been shown in early studies on MIR data from JWST (e.g., Ling et al. 2022;Wu et al. 2023).

LFs of all galaxy populations
In the four panels (7.7, 12.8, 15 and 21 m) in Figure 8 that have previous results overplotted, our LFs agree with and well extend the faint end from previous works to more than one order of magnitude.The luminosity bins in our study show little overlap with those reported in the literature, which implies that most parts of our LFs are observed for the first time, especially for  > 2. Among all the monochromatic LFs, a gradual evolution in luminosity can be seen clearly among all redshift bins.At higher redshift ( > 3), the change in LF curves is dominated by the density evolution.The overall evolution between  = 1 − 2 and  = 2 − 3 appears more moderate.At high- (> 1) our data are incomplete at low luminosities, as indicated by coloured arrows.Therefore, the faint end of the LF is not sampled and cannot be estimated.(Goto et al. 2019;Le Floc'h et al. 2005;Pérez-González et al. 2005;Babbedge et al. 2006;Rodighiero et al. 2010)

LFs of SF and AGN galaxies
While MIR emission serves as a good proxy for star-formation activity, it may be contaminated by the radiation from heated dust due to AGN activity.To study the impact of AGN on LF from different galaxy populations, we obtain LF for typical SF and AGN host galaxies separately by frac AGN (refer to §3.1).The resulting LFs are presented in Figure 9 (SF) and 10 (AGN).LF for all galaxies (from Figure 8) are also overplotted here for reference.
We find that LFs from SF galaxies generally follow those from all galaxies, exhibiting only slight differences (≲ 0.2 dex), as the majority of galaxies in our sample are SF galaxies.The main features of LFs in Figure 8 (all galaxies) can be also found in Figure 9 (LFs of SF galaxies).The disappearance of the high- bright end in Figure 9 (due to the removal of luminous AGN) implies a potentially stronger evolution in SF galaxies compared to what is suggested in Figure 8, because fewer SF galaxies are expected to be found at  > 3.
Figure 10 suggests that AGN host galaxies are less significant at lower redshift within our main luminosity coverage  * ∼ 10 8 − 10 11  ⊙ .AGNs are more dominant at the highest redshift bin at  = 3 − 5.1, as presented by previous analyses on IR populations (e.g., Gruppioni et al. 2013).Nevertheless, it is difficult to draw a solid conclusion on the AGN evolution within the error bar, because of the insufficient sample size.

TIR LF
We now show the total (bolometric) IR LF in Figure 11.The procedure to construct the TIR LF is the same as what has been described in §4.1, and we use the F1000W limit-selected sample (refer to Figure 7 and Table 5).409 galaxies in total are used to obtain the TIR LF regardless of galaxy type, i.e., all galaxies are used including AGN components, indicated by circles.For comparison, we also show the TIR LFs which only use star-forming components of the SEDs.AGN components of the SEDs (torus, polar dust, and disk emissions) are excluded in computing those LFs shown with diamond symbols, in order to examine the possible bias imposed by AGNs in interpreting star-forming activities.No significant difference is found in comparing TIR LFs that include all the SED components (circles) with those excluding AGN components (diamonds).Their medians are only < 6% different on average.This suggests that AGN emissions do not dominate the SED in the far-IR wavelengths.Note that it differs from our frac AGN definition, which only depends on emission within 3 − 30 m.As a result, the impact from the AGN components is expected to be marginal.
Comparing Figure 11 to Figure 9, we find the evolution of TIR LF traces the trend of 7.7 and 12.8 m SF galaxy LFs.The correlation between 7.7 and 12.8 m luminosity and TIR luminosity results from the strong PAH emission at these wavelengths, which plays a significant role in shaping the IR spectrum of SF galaxies.We note that this has been investigated in previous studies (e.g., Caputi et al. 2007;Bavouzet et al. 2008;Goto et al. 2011;Lin et al. 2024).
While we have pushed the TIR LF to a notably high- ∼ 5 with the brand new JWST data, we should caution that the TIR luminosity of galaxies at higher redshift ( > 3) could be undetermined.The dust emission in MIR which we trace for deriving TIR luminosity will start to be replaced by stellar emission at high-, thus the FIR (and partially MIR) part of SED will mostly rely on extrapolation of the model.However, with the reddest MIRI bands in the CEERS field (F1800W and F2100W), 3.3 m PAH features from star-forming activity can still be detected at  = 4 − 5, as discussed in Yang et al. (2023a).Furthermore, the majority of such high- galaxies are AGN (refer to §4.1.2).When considering MIRI coverage, the SEDs of these galaxies (Figure A3 and A4) are typically dominated by AGN emission rather than stellar emission, which indicates a clue to the FIR SED.Previous studies (e.g., Elbaz et al. 2011;Dai et al. 2018) have shown a redshift-independent relation between AGN emission and FIR / TIR luminosity.Besides these points, we caution readers about possible uncertainty in estimating TIR luminosity based on the MIRI flux.2013) concentrated on the deep pencil beam GOODS-S field.Both of them have similar luminosity limits with respect to their redshift range, as the markers suggest.Note that we derive the limits for Gruppioni et al. (2013) with the nominal 100 m limiting flux of SF-AGN galaxy in the GOODS-S field (1.2 mJy), assuming the SED of NGC6090.Gruppioni et al. (2020) (open diamonds) utilise the sub-mm observations from the ALMA ALPINE survey.By tracing the rest-frame FIR continuum emission, they first extend the IR LF to a wide range  = 0.5 − 6.

Comparison with the literature
Similar to the monochromatic LFs, our results are consistent with previous studies and further advance them more than one order of magnitude fainter.The depth of JWST enables us to explore the faintest MIR objects ever seen, with a luminosity limit roughly corresponding to  * ∼ 10 9.7  ⊙ ( = 0 − 1), 10 10.6  ⊙ ( = 1 − 2), 10 10.9  ⊙ ( = 2 − 3), and 10 11  ⊙ ( = 3 − 5.1).To put this into context, these limits are approximately 1.5 dex fainter compared to the previous MIR space telescope AKARI (Goto et al. 2019).In all the redshift bins, we extend the limits of Herschel works by an order of magnitude.Additionally, for higher redshifts, we push the luminosity bins to ∼ 0.5 dex fainter than those derived from Gruppioni et al. (2020).
The shape of JWST TIR LF is similar to the Herschel or AKARI works.Still, we would like to point out a deviation compared to Gruppioni et al. (2020) at  = 1 − 2 and  = 2 − 3, where our LFs are higher and steeper in the faint end.The deviation has been previously suggested in the comparison with Gruppioni et al. (2013) and Magnelli et al. (2013) in Gruppioni et al. (2020), though it was not significant according to their luminosity limits.Given that we now constrain the LF to a much fainter luminosity (for their  = 1.5 − 2.5 and  = 2.5 − 3.5 bins), this deviation is not negligible, although we must note that the redshift bins are not exactly the same (i.e., those of Gruppioni et al. 2020 are shifted up by about  = 0.5).On the other hand, we observe a good agreement between our data points and Gruppioni et al. (2020) for the highest redshift.

MCMC analysis
To quantify the evolution of the TIR LFs, we adopt a modified-Schechter function (Saunders et al. 1990) to fit the LFs: We also provide luminosity limits for SF galaxies (longer solid arrows) and all galaxies (shorter dashed arrows, obtained from Figure 8).
Figure 10.Same as Figure 9 but for AGN host galaxies (shown in triangle / solid line) only.We also provide luminosity limits for AGNs (longer solid arrows) and all galaxies (shorter dashed arrows, obtained from Figure 8).which is a power law (determined by ) for  <  * and a Gaussian in log  (determined by ) for  >  * . * is the normalised factor for density.We apply the Markov chain Monte Carlo (MCMC) method to fit our LFs.This is conducted by the Python package emcee (Foreman-Mackey et al. 2013), with 100 walkers initialised and 1000 iteration steps.Data points are allowed to move within their error bars in each iteration, assuming a Gaussian distribution.Given that our data points mostly probe the faint end of the whole LF, it is not sufficient to fit the four free parameters  * ,  * ,  and  simultaneously.To overcome the issue, we use a technique similar to Babbedge et al. (2006) and Gruppioni et al. (2020) for fitting.First, we fix the slope of the bright end  to 0.5.The value is obtained from the Herschel LF (Gruppioni et al. 2013) that has been well-constrained on bright galaxies  * = 10 11 − 10 13  ⊙ .Then, we fit the lowest redshift bin  = 0 − 1 to the remaining three parameters.A flat prior range of log( * / ⊙ ) = [8, 13], log( * /Mpc The median and 1-sigma uncertainty for fitted parameters are summarised in Table 6, where we use the 16th-and 84th percentile from the MCMC results for the uncertainty range.We overplot the fit curves within 1-sigma uncertainty (grey lines) and the median fit curve (black lines) in Figure 12-15 for each redshift bin.The median fit curve for pure SF TIR LF (dashed line) is also provided.Faint luminosity bins that are further away from the luminosity completeness limit are excluded to avoid bad fits.These bins are shown in grey.In the right panel of Figure 12-15, we show the probability distributions of the fit parameters.
The MCMC analysis suggests a strong degeneracy in  * and  * .While the slope of the faint end  is relatively constrained (right panel of Figure 12), parameters for the knee of the LF seem to be highly undetermined, especially for  * , due to the limited samples with high luminosity.We further address this issue in §4.3.From the weakly constrained  * , we also notice a clear scaling relation between  * and  * .The density ( * ) and luminosity ( * ) evolution from our medianfit TIR LF are shown in Figure 16.In addition, we overplot the best-fit curves, which parameterised the evolutions as  * ∝ (1 + ) 0.90±0.98 and  * ∝ (1 + ) −1.73±0.42 .These curves are obtained by fitting the data points with non-linear least squares, assuming a function form of (1 + )  , where  is the normalisation factor and  is the slope.The associated errors are obtained from the derived covariance matrix.
Our findings indicate that our fit broadly follows the fitted evolution presented in Gruppioni et al. (2013) but with a noticeably large error bar.While  * indeed exhibits a decreasing trend with redshift, the overall evolution of  * appears to be minor, as it falls within the uncertainty range caused by significant degeneracy.Despite this, it is important to note that the global shape of fitted LF does not show obvious changes within the 1-sigma range (grey lines) due to the strong association between  * and  * .Hence, we stress that such degeneracy is unlikely to impact the subsequent integration of luminosity density.

Luminosity density evolution
The IR luminosity density, denoted by Ω, depicts the average IR energy emitted by all galaxies within a specific age of the Universe.This is obtained by integrating the LF over the luminosity range, i.e., based on the modified-Schechter function fit () for the LF.We integrate the LF down to  * = 10 8  ⊙ .This is the same as Gruppioni et al. (2013), which has similar redshift coverage for comparison.
To derive the 1-sigma uncertainty of Ω, we compute Ω for every individual model fit from the MCMC results (Figure 12-15).Since a range of parameter sets are considered based on their probability distribution, we can ensure that the calculated distribution of Ω agrees with our LF fit.From this Ω distribution, the median, the 16th-and 84th percentile are obtained as the final 1-sigma uncertainty of Ω.
The horizontal error bar shows the redshift range of that bin, same as Figure 16.The derived values of Ω for each redshift bin are given in Table 7.
Figure 17 presents the IR luminosity density (Ω) evolution.We obtained the probability distribution of Ω by integrating MCMC results.For better visual interpretations, we overplot an associated "violin" (purple shades, centred at the bin) for each redshift bin, which displays the kernel density of Ω in a symmetric shape.The width of the violin at a specific Ω value indicates the corresponding probability density.We observe that the evolution of Ω well agrees with previous observational works (Magnelli et al. 2013; Gruppioni     al. 2013;Goto et al. 2019;Gruppioni et al. 2020) within their 1-sigma uncertainties, except for the  = 1 − 2 bin.Our data points are generally higher than the model fit from Kim et al. (2024), but the two error bars just touch each other.The long tail distribution of Ω suggested by the violin plot apparently affects the estimation of its median and 1-sigma, because of the degeneracy in parameter fitting.If their mode, instead of median, is considered, the results would be more aligned with the literature.We attributed the observed excess of Ω at  = 1 − 2 to issues in the model fitting process that caused deviations in the fitted LF shape and consequently estimate of Ω.This is an example of fitting's  inability to determine the break luminosity  * accurately.As Figure 16 shows, an unusual bump in  * evolution can be identified at  = 1 − 2. Compared to Magnelli et al. (2013) and Gruppioni et al. (2013) which have been well constrained at this range, the bright end of our LF appears to be overestimated (refer to Figure 13).Given our limited information on the bright end of the LF,  * is highly sensitive to the last luminosity bin that typically lies near the knee of our LF.In the case of  = 1 − 2, the last luminosity bin is relatively high with respect to the literature, potentially contributing to the observed phenomenon.We also note this effect is not significant for LF at other redshifts because their fitted bright end is generally consistent with the literature.
Our work pushes the MIR-based IR SF history to  = 4.05 (∼ 1.52 Gyr after the Big Bang).While IR SF history has already probed to  ∼ 6 with sub-mm data (Gruppioni et al. 2020), it is still a sig-nificant improvement compared to the previous-generation IR space telescopes, where the highest redshifts of TIR LF they reached are  ∼ 1.4 (Goto et al. 2019, AKARI) and  ∼ 3.6 (Gruppioni et al. 2013, Herschel).Our results highlight the SFRD peak and the turnover at  ∼ 1.5, and further reveal the dust-obscured SF history before the cosmic noon.
We next investigate the luminosity density evolution in 7.7 m due to its characteristic PAH emission.Although  7.7 is expected to trace the star-formation rate and also  TIR , testing how far this correlation holds is of importance in understanding the relation between the far-IR dust emission and those from PAH molecules.Also, quantifying the density evolution of 7.7m emission is crucial to plan future larger mid-IR surveys with JWST and future telescopes such as the Far-IR Spectroscopy Space Telescope (FIRSST) and PRobe far-infrared mission for astrophysics (PRIMA, Moullet et al. 2023).The density evolution of this specific luminosity has been previously studied by Goto et al. (2010).By fitting 7.7 m LFs in Figure 8 and 9 with the same method described in §4.2.2 and above, we compute the Ω 7.7  for all galaxy and SF-only galaxy samples.The result is shown in Figure 18.First, we note our Ω 7.7  evolution agrees with the trend of Goto et al. (2010) in their redshift range  = 0 − 2. A peak, similar to Ω, has been found at  ∼ 1.5.Nevertheless, our Ω 7.7  is severely affected by the insufficient data points in 7.7 m LF that give huge degeneracy and uncertainty in the fitting.While the alignment of the Ω evolution of the two samples is expected as §4.1.2discussed, no quantitative conclusions can be drawn on the differences or ratios between the two considering the overlapping error bars.Since PAH emission has served as an important indicator of SF activity, studying PAH evolution and its contribution from SF / AGN galaxies is crucial for a comprehensive understanding of its behaviour and robustness across cosmic time.
Here, we present a first step towards higher redshifts and encourage future works to investigate the PAH evolution with better statistics.
It certainly requires a more extensive mid-IR survey to follow up on and confirm our results.We suggest utilising future mid-IR observations by JWST with a larger survey area such as Rieke et al. (2017).From a sample of thousands of galaxies, we expect to obtain the LF with higher resolution in redshift and reduced error.

CONCLUSIONS
In this work, we show the first JWST rest-frame total IR LFs and MIR LFs at 7.7, 10, 12.8, 15, 18 and 21 m at  = 0 − 5.1 with MIR images from the CEERS survey.While our LFs are consistent with previous studies on MIR LFs, we also present the LFs up to the highest redshift  > 4 and to 2 orders of magnitude fainter in terms of IR luminosity.For TIR LF, we found an overall density evolution of ∝ (1 + ) −1.73±0.42 .Due to the degeneracy in parameter fitting, the luminosity evolution shows a flatter shape of ∝ (1 + ) 0.90±0.98 .The turnover of SFRD occurred at  ∼ 1.5.With these results, we demonstrate the potential that JWST can bring to future research on IR galaxy evolution.Still, we should exercise caution, as our LF may be undetermined at high- because MIRI only sampled the near-IR continuum, mainly contributed by stellar emission.While scenarios that could reduce this effect have been discussed, supporting observations are surely needed for a more precise determination of JWST galaxy IR LFs at high-.Given that our LFs are based solely on photo- measurements, spec- from the JWST NIRSpec and MIRI MRS will be crucial to minimise the redshift uncertainty in the calculation of IR luminosities and LFs for those faint JWST galaxies.
Since there is currently no far-IR space telescope with similar sensitivity to JWST, counterpart observations in the sub-mm with ALMA bands will be essential to constrain the rest-frame far-IR SEDs.This is particularly effective for high- galaxies because the lowest luminosity bin of our LF at  = 3−5.1 is already the same order of magnitude (by only ∼0.5 dex) as the LF reported by Gruppioni et al. (2020) from ALMA data.Furthermore, utilising the reddest MIRI band, F2550W (25.5 m) in future surveys will also help to extend the redshift range for probing rest-frame MIR emission.These findings motivate further MIRI and corresponding sub-mm observations in order to deepen our view of the obscured star-forming history of the Universe.(Magnelli et al. 2013;Gruppioni et al. 2013;Goto et al. 2019;Gruppioni et al. 2020;Kim et al. 2024) are also provided.SFR density ( SFR ) converted from Ω is also shown for reference, assuming the relation from Kennicutt (1998).MNRAS 000, 1-18 (2024)

Figure 2 .
Figure 2. The cumulative areal coverage of the sum of all 4 pointings as a function of 1-sigma pixel error for each band.

Figure 4 .
Figure 4.The reduced  2 histogram of the fitting result.The black line is the median (1.39), and the red dashed line shows the criterion of  2 = 5.

Figure 5 .
Figure 5.Estimated photo- from cigale ( CIGALE ) versus redshift from the CANDELS-EGS catalogue ( CANDELS−EGS ), where  CANDELS−EGS can be either spec- or photo-.The photo- value from the catalogue is used if spec- is unavailable.The two redshift types are presented separately in the upper panel, and outliers for photo- are shown in grey crosses.1-to-1 line (black line) and 15% uncertainty (grey region) are also plotted.The lower panel shows only sources with spec- available in the catalogue and is magnified to the proper range (indicated by the red square in the upper panel) for clarity.

Figure 6 .
Figure 6.The redshift ( CIGALE ) histogram of the sources in each band, divided into our redshift bins.The colour regions indicate the span of each bin, and the dot-dashed lines are the median.
Figure 8 (monochromatic LFs for all galaxies) b Figure 9 (monochromatic LFs for SF galaxies) c Figure 10 (monochromatic LFs for AGN) d Figure 11 (TIR LFs for all galaxies)Table5.The final sample size for monochromatic and TIR LFs based on 80% completeness limit selection in each filter and galaxy population.
[0, 1], [1, 2], [2, 3], [3, 5.1] in each LF.These limits are determined by the brightest among all limiting luminosities of different SED types (i.e., SF and AGN) assumed in the given redshift bins, as illustrated in Figure 7.The luminosity limits are shown along with the corresponding LFs in §4.1 and §4.2 to help readers confirm at what luminosity the LFs are confident.

Figure 7 .
Figure7.TIR luminosity ( TIR ) as a function of redshift ( CIGALE ).The solid line is the completeness luminosity limit converted from 80% completeness flux limits for the F1000W band, assuming SF (blue) and AGN (orange) SED template.Galaxies below the line (depending on their SED type) are excluded from the analysis.The coloured arrows on the left indicate the TIR luminosity limits for each redshift bin used in this work, which is the brightest among all limits for different SED types in that redshift bin.

Figure 12 -
Figure12-15 compare our LFs to the literature, separated in each redshift bin.The luminosity completeness limits for our work and literature are denoted by different markers.Fitted curves of our LF are also presented, refer to Section 4.2.2.BothGruppioni et al. (2013) (open/closed coloured circles) andMagnelli et al. (2013) (coloured stars) use the FIR observations from Herschel to produce TIR LF.Specifically,Gruppioni et al. (2013) investigated an extensive field containing the GOODS, ECDFS and COSMOS area, andMagnelli et al. (2013) concentrated on the deep pencil beam GOODS-S field.Both of them have similar luminosity limits with respect to their redshift range, as the markers suggest.Note that we derive the limits forGruppioni et al. (2013) with the nominal 100 m limiting flux of SF-AGN galaxy in the GOODS-S field (1.2 mJy), assuming the SED of NGC6090.Gruppioni et al. (2020) (open diamonds) utilise the sub-mm observations from the ALMA ALPINE survey.By tracing the rest-frame FIR continuum emission, they first extend the IR LF to a wide range  = 0.5 − 6.Similar to the monochromatic LFs, our results are consistent with previous studies and further advance them more than one order of magnitude fainter.The depth of JWST enables us to explore the faintest MIR objects ever seen, with a luminosity limit roughly corresponding to  * ∼ 10 9.7  ⊙ ( = 0 − 1), 10 10.6  ⊙ ( = 1 − 2), 10 10.9  ⊙ ( = 2 − 3), and 10 11  ⊙ ( = 3 − 5.1).To put this into context, these limits are approximately 1.5 dex fainter compared to the previous MIR space telescope AKARI(Goto et al. 2019).In all the redshift bins, we extend the limits of Herschel works by an order of magnitude.Additionally, for higher redshifts, we push the luminosity bins to ∼ 0.5 dex fainter than those derived fromGruppioni et al. (2020).The shape of JWST TIR LF is similar to the Herschel or AKARI works.Still, we would like to point out a deviation compared toGruppioni et al. (2020) at  = 1 − 2 and  = 2 − 3, where our LFs are higher and steeper in the faint end.The deviation has been previously suggested in the comparison withGruppioni et al. (2013) andMagnelli et al. (2013) inGruppioni et al. (2020), though it was not significant according to their luminosity limits.Given that we now constrain the LF to a much fainter luminosity (for their  = 1.5 − 2.5 and  = 2.5 − 3.5 bins), this deviation is not negligible, although we must note that the redshift bins are not exactly the same (i.e., those ofGruppioni et al. 2020 are shifted up by about  = 0.5).On the other hand, we observe a good agreement between our data points andGruppioni et al. (2020) for the highest redshift.

Figure 9 .
Figure 9. Same as Figure 8 but for star-forming (SF) galaxies (shown in star / solid line) only.Galaxy LFs for all populations from Figure8are also plotted (shown in circle / dashed line) for reference.We also provide luminosity limits for SF galaxies (longer solid arrows) and all galaxies (shorter dashed arrows, obtained from Figure8).

Figure 11 .
Figure 11.The rest-frame total infrared (TIR) luminosity functions.The four redshift bins  = 0 − 1 (blue),  = 1 − 2 (green),  = 2 − 3 (orange), and  = 3−5.1 (red) are shown.TIR LFs derived from all galaxies (including AGN components in the SEDs) are marked by circles.For comparison, TIR LFs derived from using only star-forming components of the SED (i.e., excluding the contribution from AGN components of the SEDs such as torus, polar dust, and disk emissions) are marked by diamonds, with an offset of 0.125 dex in  TIR for clarity.The coloured arrows at the top show the luminosity completeness limits, adopted from Figure 7.The data points below the limits are shown as open markers and are not connected.
−1], and  = [−1, 3] are set.We determined and fixed  from the fitting result, which means only  * and  * are fit for the rest of the redshift bins.

Figure 12 .
Figure 12.Left panel: the rest-frame total infrared (TIR) luminosity functions, same as Figure 11 but for  = 0 − 1 bin only.The TIR LF are marked in black open circles.Grey open circles show luminosity bins not used in fitting.The median model fits are plotted in black solid lines, and the grey lines are the fits within the 1-sigma uncertainty range of the parameters.TIR LFs from previous works (Magnelli et al. 2013; Gruppioni et al. 2013; Goto et al. 2019; Gruppioni et al. 2020) with similar redshift range are overplotted.Markers on the top indicate the luminosity completeness limits of literature with the same marker, where the black arrow is for our work.Right panel: the corner plot showing the probability distribution of fit parameters from MCMC analysis.The median (blue cross), as well as the 16th-and 84th percentile of fit parameters, are provided.

-
Table 7. IR and 7.7 m luminosity density evolution, taken from Figure 17 and 18, respectively. et

Figure 16 .
Figure16.The redshift evolution of  * and  * from the fitted TIR LF, taken from Table6.The horizontal error bar is the redshift range of the bin.The model fit fromGruppioni et al. (2013) is provided (grey dashed line), and our best fits for the data points are shown in orange dashed lines.

Figure 17 .
Figure 17.The redshift evolution of Ω (open circle).The purple shades are the violin plots that illustrate the probability distribution of Ω. Ω evolution from the literature(Magnelli et al. 2013;Gruppioni et al. 2013;Goto et al. 2019;Gruppioni et al. 2020;Kim et al. 2024) are also provided.SFR density ( SFR ) converted from Ω is also shown for reference, assuming the relation from Kennicutt (1998).

Figure 18 .
Figure 18.The redshift evolution of Ω 7.7  for all galaxy (blue open circle) and only SF galaxy (orange open star).Comparison to Goto et al. (2010) is provided.

Figure A1 .
Figure A1.Example fitted galaxy SEDs (in observed frame) from cigale for  = 0 − 1 bin, with their photo- and reduced  2 .The left panel shows a star-forming galaxy and the right panel shows an AGN host galaxy.The IAU designation from the CANDELS-EGS catalogue is provided as id.

Table 1 .
The 80% completeness limit for each pointing and band.

in the First two pointings Last two pointings Band o001_t021 o002_t022 o012_t026 o015_t028
. Built from the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011), the CANDELS-EGS catalogue provides broad coverage of wavelengths from near-UV to MIR (0.4 m − 8 m)

Table 2 .
The reliability beyond the 80% completeness limits for each pointing and band.

Table 3 .
to better estimate the photometric redshifts (photo-) and SEDs for sources in the compiled catalogue by incorporating JWST mid-IR photometry.cigale is a Python code that can Depths in each filter in the compiled catalogue.The 5 depths of the CANDELS-EGS catalogue band are taken fromStefanon et al. (

Table 4 .
The modules and parameters used in cigale.