Cosmic star-formation history and black hole accretion history inferred from the JWST mid-infrared source counts

With the advent of the James Webb Space Telescope (JWST), extra-galactic source count studies were conducted down to sub-microJy in the mid-infrared (MIR), which is several tens of times fainter than what the previous-generation infrared (IR) telescopes achieved in the MIR. In this work, we aim to interpret the JWST source counts and constrain cosmic star-formation history (CSFH) and black hole accretion history (BHAH). We employ the backward evolution of local luminosity functions (LLFs) of galaxies to reproduce the observed source counts from sub-microJy to a few tens of mJy in the MIR bands of the JWST. The shapes of the LLFs at the MIR bands are determined using the model templates of the spectral energy distributions (SEDs) for five representative galaxy types (star-forming galaxies, starbursts, composite, AGN type 2 and 1). By simultaneously fitting our model to all the source counts in the six MIR bands, along with the previous results, we determine the best-fit evolutions of MIR LFs for each of the five galaxy types, and subsequently estimate the CSFH and BHAH. Thanks to the JWST, our estimates are based on several tens of times fainter MIR sources, the existence of which was merely an extrapolation in previous studies.


INTRODUCTION
Extragalactic source counts (defined as the number count of sources as a function of flux density) have been used as one of the basic cosmological tests because it is simple and straightforward.It is a practical tool to check the variation of the source distribution with increasing luminosity distance (or redshift).In the early studies in the infrared (IR), such as Elbaz et al. (1999), Serjeant et al. (2000), and Pearson (2005), the detection limit at 15m was at the sub-milliJansky (sub mJy) level based on the Infrared Space Observatory (ISO) observations in the 1990s.Thanks to the AKARI (Murakami et al. 2007) and Spitzer (Werner et al. 2004) space telescopes, the detection limit became deeper, reaching a few tens of Jy in subsequent studies by Pearson et al. (2010), Takagi et al. (2012), and Pearson et al. (2014).
Much deeper IR science data from the James Webb Space Telescope (JWST, Gardner et al. 2006;Kalirai 2018, which was launched in late December 2021), has been publicly available, e.g., Stephan's ★ E-mail: seongini@gmail.comQuintet 1 and the Cosmic Evolution Early Release Science Survey 2 (CEERS; Finkelstein et al. 2017) since July 2022.This availability allowed for new source count (SC) tests using very faint (  < Jy) mid-infrared (MIR) sources.Ling et al. (2022Ling et al. ( , 2023) ) presented the source counts at the F770W, F1000W and F1500W bands of the Mid-IR Instruments (MIRI, Rieke et al. 2015) using early data from the Stephan's Quintet field.Additionally, Wu et al. (2023) utilized data from the CEERS field to provide additional source counts at the F1280W, F1800W, and F2100W bands.In the previous SC studies with the former IR telescopes (e.g., Elbaz et al. 1999;Serjeant et al. 2000;Malkan & Stecker 2001;Lagache et al. 2005;Pearson 2005;Pearson et al. 2010), one of the key issues was understanding the properties of the sources contributing to the cosmic infrared (IR) background (CIB, Puget et al. 1996;Lagache et al. 2005).The recent JWST source counts have reached the flux limits below the Jy level (e.g., at 7.7, 10, and 12.8 𝜇m;Ling et al. 2022Ling et al. , 2023;;Wu et al. 2023), revealing galaxies several tens of times fainter.The JWST achieves 0.3 arcsec resolution at 7.7m (F770W), while the Spitzer has a1 https://webbtelescope.org/contents/media/images/2022/034/01G7DA5ADA2WDSK1JJPQ0PTG4A 2 https://ceers.github.io/overview.htmlpoint spread function (PSF) of 2 arcsecs at 8 m (IRAC 4) (Gardner et al. 2006;Werner et al. 2004).Therefore, the JWST has seven times better resolution compared to the Spitzer space telescope.The sources contributing to the IR background in previous works are now resolved into discrete sources.
Half, at least, of the energy in the Universe, is emitted in the IR wavelengths, with a peak in the far-IR, providing insight into the overall bulk of the IR energy (i.e., total IR luminosity, Lutz 2014; Casey et al. 2014).In the mid-IR (MIR) regime, where emission is dominated by rather warmer/hotter dust heated by strong radiation, various emission features exist, such as polycyclic aromatic hydrocarbons (PAHs) or silicate features, which offer clues to the detailed dust properties associated with star-formation (SF) activity (Tielens 2008;Desai et al. 2007;Takagi et al. 2010;Kim et al. 2019;Sirocky et al. 2008, Ohyama et al. 2018).Additionally, MIR observations probe the hot dust heated by active galactic nucleus (AGN) activity (Gruppioni et al. 2008;Wang et al. 2020;Lacy & Sajina 2020).Importantly, the advantage of IR observation is that the data is less affected by dust obscuration.
Luminosity function (LF) is a statistical description of the distribution of galaxy luminosities and is defined as the number of galaxies in a specific volume as a function of luminosity (Johnston 2011).It is a measure of the number density of galaxies at different luminosities, and provides important information about the populations of galaxies and their evolution (e.g., Elbaz et al. 1999;Rowan-Robinson 2009;Goto et al. 2010;Gruppioni et al. 2011).One of the key features of the luminosity function is its shape, which is affected by a variety of evolutionary factors such as star formation activity, galaxy evolution and growth, the role of supermassive black holes (SMBHs), and so on (e.g, Lagache et al. 2003, Gruppioni et al. 2013).More accurate measurement of LFs has been one of the main concerns in observational cosmology, to obtain more accurate cosmic star-formation history and galaxy evolution (e.g., Gruppioni et al. 2010;Goto et al. 2015, Kim et al. 2015;Kilerci Eser & Goto 2018;Goto et al. 2019).The source count (SC) is closely related to the luminosity function of galaxies, which describes the distribution of galaxies in terms of observed brightness (e.g., Gruppioni et al. 2011;Pearson et al. 2014).If redshift information is obtained (measured or assumed), astronomers can construct multi-wavelength luminosity functions by combining information from source counts at different wavelengths.This can provide a more complete picture of the galaxy populations.Therefore, the source count is an important tool in studying galaxy evolution and luminosity functions.
Previous source count studies have shown that IR luminosity functions have to change dramatically with redshift and the "no-evolution model" can not describe the bump feature around 0.1-1mJy in the source counts (Elbaz et al. 1999;Serjeant et al. 2000;Lagache et al. 2003;Pearson et al. 2010).Given the crucial role of evolution, modelling IR source counts becomes essential to understand the impact of galaxy evolution on counts and gain insights into the dusty star formation history (Rowan-Robinson 2009;Gruppioni et al. 2011).In this work, we take much fainter galaxies into account to understand the connection between galaxy source counts and cosmic evolution history.Our primary goal is to reproduce the observed source counts across the contiguous MIRI bands of the JWST, and derive constraints on the cosmic star formation history (CSFH) and black hole accretion history (BHAH).We aim to develop a model that effectively describes the faint end (∼ sub-Jy) of the MIR source counts obtained by the JWST, while also considering previous results from the former IR space telescopes that cover up to a few tens of mJy.Utilizing the model we obtain, we interpret source count results in terms of the evolution of different galaxy populations contributing to the source counts across a wide range of flux densities.Additionally, the CSFH and BHAH will be updated using new parameters obtained including very faint IR galaxies.
This paper is structured as follows.In Section 2, we summarise the source counts from the JWST as well as from previous works.Section 3 shows how to derive our source count model in the MIR and compares them with observed source counts.In section 4, we discuss the model we obtain and infer CSFH and BHAH.We summarise our work in section 5, In the model computation, we assume the Salpeter initial mass function (IMF, Salpeter 1955) and Λ cold dark matter (ΛCDM) cosmology with H 0 = 70km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7.
Their source count results showed roughly reasonable agreement with the extension/prediction of the previous models, which tells us that our understanding of the IR Universe has been consistent with previous expectations.However, there are slight deviations in their comparisons, especially in the fainter range, which a simple extension of the old models can not give sufficient explanations for.A large error bar of source count in a certain flux bin comes from a small number of sources contributing to the flux bin, which is determined based on the Poisson error for small-number statistics (Gehrels 1986).Their observed number counts show occasional scatter and fluctuation of the data points in a certain range of flux (Wu et al. 2023), however, it is not certain if it is because of statistical variance/fluctuation between the small patches of the different sky fields (cosmic variance).We attempt to see how the evolutionary parameters have to be changed/adjusted to describe the new source count results obtained from the JWST observations.

Source counts from the previous IR telescopes
Since Oliver et al. (1997) provided the 15m ISOCAM source counts from the HDF field, many MIR source count studies have been carried out based on the ISOCAM (Elbaz et al. 1999;Gruppioni et al. 2002;Pearson 2005, etc.).AKARI's large area survey (e.g., North Ecliptic Pole, Wada et al. 2008;Kim et al. 2012) and Spitzer (e.g., EGS, SWIRE, and CDF-S field, Fazio et al. 2004;Shupe et al. 2008;Papovich et al. 2004) data were used to derive the MIR source counts (Takagi et al. 2010;Pearson et al. 2010Pearson et al. , 2014;;Murata et al. 2014;Davidge et al. 2017).In this work, we also use the MIR source counts from previous works when we derive new evolution parameters.The previous source counts from the AKARI bands provide good comparison samples because wavelength coverage is similar to that of MIRI, and some of them have almost the same effective wavelengths (e.g., 15 and 18 m).In this work, we accept the previous source counts data only where the completeness for the source detection is higher than 80%.

Local luminosity functions and evolution models
To obtain physical motivations on the parameters of LFs, we attempt to model/reproduce the mid-IR source counts by combining the JWST/MIRI source counts (Ling et al. 2022(Ling et al. , 2023;;Wu et al. 2023) with those obtained by previous works (Oliver et al. 1997;Serjeant et al. 2000;Pearson et al. 2010;Takagi et al. 2012;Pearson et al. 2014;Davidge et al. 2017).
We use one of the simple/popular ways, taking the backward evo-lution method by Gruppioni et al. (2011), where evolution models are described by evolution parameters given for different galaxy types.The evolution of the local luminosity functions (LLFs) towards high redshifts has been one of the frequently used methods (e.g., Rowan-Robinson 2009;Béthermin et al. 2011;Gruppioni et al. 2011, etc.) because this empirical approach reasonably fits the observed galaxy number counts and is consistent with theories.Gruppioni et al. (2011) present the model LLFs for mid-IR bands (e.g., 12 and 15 m): Φ(), described by the following function (Saunders et al. 1990, modified Schechter function): which behaves as a power law for  <<  ★ and as a Gaussian in log for  >>  ★ .One of the main reasons we use the methods from Gruppioni et al. ( 2011) is that we can take advantage of testing different types of galaxies (see Fig. 1) because Gruppioni et al. (Fig. 1 and Table 1 of 2011) provide the LLFs for five representative populations: normal star-forming galaxies (SFGs), starbursts (SB), composites -a mixture of SFG and low-luminosity AGN (LLAGN), AGN type-2 (obscured AGN; AGN2), and AGN type-1 (unobscured AGN; AGN1).First, we obtain LLFs for six MIR bands of the JWST/MIRI (F770W, F1000W, F1280W, F1500W, F1800W, and F2100W)based on the given information concerning 15m LLFs for each galaxy population (i.e., the LF parameters such as  ★ , Φ ★ , , and  in Table 1 of Gruppioni et al. 2011), the LLFs for the other wavelength bands are derived using the filter bands of the MIRI and model spectral energy distributions (SEDs).We can reproduce the model 15m LLFs presented in Gruppioni et al. (2011) using the filter response of the F1500W and galaxy model SEDs (e.g., Polletta et al. 2007).While doing this, we can also use the other MIRI filters (F770W, F1000W, F1280W, F1800W, and F2100W) to obtain LLFs at these bands at the same time.
For this procedure, we take the SED model templates from Polletta et al. (2007), as shown in Fig. 1.As in Gruppioni et al. (2011), we applied the template variations for SFG and SB according to the luminosity range, where the slope of the luminosity function changes rapidly.For the SFG population, we take the assumption that an SED evolves with luminosity from Sa ( 15m < 10 9 L ⊙ ) to Sdm ( 15m ≥ 10 10 L ⊙ ).For the SB population, the SED varies from a moderate NGC6090 ( 15m ≤ 10 10 L ⊙ ) up to extreme Arp 220 ( 15m > 10 11.9 L ⊙ ).The template variation is discrete as shown in Fig. 1.For the composites, AGN1, and AGN2 populations, the templates of Seyfert-2, TQSO1, and Mrk231 are taken, respectively.Using these templates, we obtained local luminosity functions (LLFs).In Fig. 2, we show the model LLFs at the six JWST/MIRI bands (7.7, 10, 12.8, 15, 18, and 21𝜇m) for five different galaxy types.

Evolution of LLFs and number counts
For the evolution of the LLFs for galaxy populations, we follow the evolution model with the parameterization that regulates the shapes of the evolution functions as shown in Gruppioni et al. (2011): where erf d is the "error function".  and  Φ are the normalisations for luminosity evolution  ★ () and density evolution Φ ★ (), respectively. regulates the shape of the function, whose combinations with  (scale factor) define the evolution peak and skewness (asymmetric shape) of the function towards high redshift.
In this work, the model LLFs derived in the six bands of the JWST evolve in terms of luminosity and density up to  ∼8 according to equations ( 2) and (3).Gruppioni et al. (2011, in their Table 2) present these 4 parameters (  ,  Φ , , and ) for 5 different types of galaxies, i.e., 20 parameters in total.Here, we aim to update/renew these evolution parameters by fitting the faint end of the MIR source counts obtained from the JWST.We attempt to encompass a broader redshift range, reaching up to =8 (which corresponds to the reionisation era) in our calculation, while Gruppioni et al. (2011) covered up to =5.
Therefore, in this work, the mixture of the five types of galaxies is given in terms of their LFs (Gruppioni et al. 2011).They evolve differently according to equations ( 2) and ( 3) with 20 parameters.For example, we show the evolution of the 15 m LLFs for five galaxy types (SFG, SB, composite, and AGN2 and AGN1) in Fig. 3.The redshift interval is 0.02, and the flux density interval is given in 0.1 dex in the calculation.
If we convert the luminosity to the observed flux density, then each panel will give the number of galaxies as a function of flux density (S  ) in the observed frame at 15 m as shown in Fig. 4 (Δ log   = 0.1).The sum of all the individual SC is given in a broken line on the top of each panel: green dashed (SFG), cyan dashed (starburst), red dot-dashed (composite), magenta dot-dot-dashed (AGN2) and blue dashed (AGN1).These plots clearly show that the brightest galaxies tend to be local ones (given in black/blue) whereas the fainter galaxies tend to be at higher redshift (in red/orange).Especially for SB and AGN2, high redshift galaxies are dominant in number at the faint end.On the bottom right panel, these five summed LFs for each galaxy type are finally integrated into a thick grey line.This grey line shows the distribution of all galaxies which are accumulated and projected to the observed frame.Therefore, counting the number of galaxies along this grey curve is the same as counting the number of galaxies projected on the night sky.These number counts of the model galaxies are shown in Euclidean normalised values (dN/dS S 2.5 ).We will compare these SC models with the observed source counts in each band (in Sec 4.1).

Model fit to the observed source counts
We perform a comparison between our model and the observed source counts by utilizing a curve-fitting procedure.Initially, we define functions that generate source count models for each JWST/MIRI band, incorporating galaxy spectral energy distributions (SEDs) and initial LLFs that follow the backward evolution of luminosity and density (as described in equations 2 and 3, in Sec 3.2).In total, we prepare six user-defined functions for the six MIRI bands.These source count models, along with the evolution parameters, are then fitted simultaneously to the observed source counts in each band to obtain the best-fit results.To facilitate this process, we employ the 'MPFIT' package (Markwardt 2009(Markwardt , 2012)), which utilizes the Levenberg-Marquardt technique for solving the least-squares problem.
As an initial guess/input in our model, we adopt the evolution parameters provided by Gruppioni et al. (2011) -their best estimates for the parameters as listed in Table 2. To prevent potential instabilities and catastrophic runaways of parameters during the fitting, we .impose limits on the parameter ranges based on preliminary trials (e.g., testing a single parameter variation with the others fixed at the initial values) using the following -larger values for the parameters   ,  Φ , and  results in the SC model shifting upward, while increasing  causes the overall height of the SC model to decrease.Furthermore,  should not be zero, as it appears as the denominator in the error functions in equations ( 2) and ( 3).Additionally, it should be also noted that a single population alone, without considering other population components, cannot fit the observed data points, highlighting the importance of multiple components in the model, although these are somewhat arbitrary.
For the data points to fit, we take the observed source counts from the literature mentioned in Sec.2.The errors in the observed source counts are also taken as input, which are used as weight in  3) and new parameters obtained from the fitting in this work.We show 15m LFs at redshifts  = 1.0, 2.0, 3.0, ..., 7.0, and 8.0 only, for better visibility.The dotted curves indicate the LFs with initial parameters, which are from Figure 3.The solid curves are the LFs with newly derived parameters.The inset (a small box) in each panel shows luminosity evolution (solid curve) and density evolution (dot-dashed curve) as a function of redshift.The y-axis 'log 10 (Evol)' represents the evolution terms as a function of , except for the first terms (log 10 L ★ (0) and log 10 Φ ★ (0)) in equations 2 and 3. Therefore, luminosity evolution is log 10 L ★ (z)-log 10 L ★ (0).In this small box, black curves denote evolution with initial parameters, and magenta curves show evolution with new parameters.Therefore, LFs in dotted curves are based on the evolution in black, while LFs in solid curves are based on the magenta curves in the insets.
the curve fitting procedure.To construct a comprehensive dataset, we combine the JWST source counts with those from previous works mentioned in Sec.2.Specifically, the AKARI 11m (S11) source counts are combined with the source count at 10m (F1000W of the JWST).The ISO 12m source counts (Rocca-Volmerange et al. 2007) are combined with the JWST 12.8m (F1280W) source count.The previous 15 and 18m source counts are merged with the F1500W and F1800W source counts of the JWST, respectively.Fazio et al. (2004, in Fig. 1 and Fig. 2) demonstrated that most of the detected sources in the bright ranges (8m) are stars.Pearson et al. (2010, in Table 2) showed that even at the MIR bands (e.g., 15 m), the stellar sources occupy more than 50 % (at 18.6 mJy) of the source detection.Davidge et al. (2017) pointed out that Fazio et al. (2004) estimated stellar contribution statistically by modelling because the observed area was not fully covered by optical bands at the time.It should be noted that the reliability of the stellar subtraction method can impact the SCs at the brightest flux ranges.Consequently, we did not include these data points in our fitting process, as indicated by grey data points (in Fig. 9).After the fitting, we see that the data points presented in grey exceed the best-fit model in the brightest flux range.

Best-fit model and new evolution parameters
We perform a simultaneous fitting of the source counts across all six MIRI bands, enabling us to derive new evolution parameters (see   ,  Φ , , and  in Table 1) and obtain the best-fit models for the source count in each band.In Fig. 5, we show new evolution models for 15m LLFs and compare them with the initial LLFs to illustrate how the new parameters have changed the initial evolution model.The solid curves represent new LFs, and dotted curves represent the initial LFs, at redshifts  = 1, 2, 3, ..., 7, and 8 (see colour bar on the bottom right).Additionally, to aid understanding, we add insets to show how luminosity (i.e., equation 2, solid line) and density (i.e., equation 3, dot-dashed) evolutions are changed explicitly -the evolution with initial parameters is given in black, and new evolution models are shown in magenta.In particular, in Fig. 6, we showcase a comparison of our model LFs (at different redshifts), against the observed LFs, which are measurements from the AKARI's NEP-Wide field (Kim Similarly to Fig. 5 for LLFs, we present Fig. 7 to illustrate how the number counts have changed with the new parameters.The dotted curves represent the initial number counts while the solid curves represent the new ones based on the new parameters (presented in Table 1).The new parameters exhibit slight deviations from the original values.These minor changes result in slightly different models for each galaxy population.Although these differences may not lead to a significant alteration in the source count model, they enable us to have a slightly different description of the JWST source counts.
In Fig. 8, we present a detailed comparison of our source counts models (red curve) against the model from Gruppioni et al. ( 2011).The SC model from Gruppioni et al. (2011, Fig. 5) is presented in Fig. 8, using a thick light brown curve.Overall, the shapes of both models fit well to the characteristic bump near a few 0.1 mJy and show a good agreement, although they show a difference at   > 2mJy.For comparison, the slope of the light brown curve (Gruppioni et al. 2011) is extended towards the faint end.This extension (indicated by a dashed brown line) is formed by summing up each extrapolation for five populations.It is slightly higher than our model.Between 2-5Jy, both models are consistent, but below Jy, extension of Gruppioni et al. (2011) goes higher again.The updated parameters do not cause a dramatic variation in the final SC model, but they introduce slight Table 1.Parameter changes after the fitting procedures -the evolutionary parameters as the initial guess in this work (the best estimates for the parameters in Gruppioni et al. (2011)) have changed as shown in the table.Fitting results/errors are from a software package, 'MPFIT'.Most of them are fitted within the parameter ranges we set, except for   for SB which reached the imposed limit (8.50) to prevent catastrophic behaviour.Gr+11 denotes Gruppioni et al. (2011)  variations in each component/population.The combination of these slight differences eventually enables a slightly better description of the JWST source counts.The slight change means that the previous parameters have also been remarkably good.Hence, the parameter fitting with the JWST data serves as a fine-tuning procedure for the evolution parameters.In Fig. 9, we show the results of the simultaneous fits to all observed source counts in six MIRI bands.The resultant best-fit models are given by thick red curves, which are the integrated results of five galaxy populations (based on the best-fit parameters).The distribution of the observed source counts from previous works mainly describes the shape of the bump (see especially 15 and 18 m), while the JWST SCs show how the behaviour of the faint-end slope has to be fitted.Overall, our source counts models peak around 0.1-0.5 mJy.Towards the faint end, they monotonically decrease describing the source counts in the JWST bands.From the bump (∼ 0.3 mJy) towards the bright end (∼ 20 mJy), our models slowly decrease or become flat (e.g., 7.7 and 12.8 m).We, however, noticed some of the observed source counts (for example, the square and plus symbols in 7.7, 10 and 15m panels, presented in grey) rather increase.We regard this related to stellar subtraction.This is the reason we ruled out the SC data from Fazio et al. (2004, Spitzer 8𝜇m SC) for 7.7m, ADF-S field data (Davidge et al. 2017) for the 7.7 and 10 m fitting, and ones from NEPW field (Pearson 2005) for 15 m.In panel (d), the triangles are excluded because their completeness is lower than 80%.The previous 24m source counts (e.g., Papovich et al. 2004) might be available, but we do not use them in the 21m fitting procedure because the effective wavelength is quite different.
Our model shows that the bump shapes are strongly dependent on the composite population, shown by the brown dot-dot-dashed curve in Fig. 9. From the bump towards the faint end where the JWST data points are distributed, the shape of the best-fit curve is shaped mostly by the composite and SFG populations.Near the faint end (80% completeness limit is indicated by vertical dot-dashed grey line), SFGs dominate by a few factors on average (the SFG has a dominant influence on the faint and bright end).The other populations (i.e., SB, AGN2, and AGN1) do not significantly influence the simultaneous variation of the bump shape and faint end slope (compared to the composite).This is because the relative ratios originated from their different LFs as well as the different evolution styles.A larger number of data points around the bump feature in 15m (or 18m) source counts may exert stronger constraints in the fitting procedure.A possible factor that slightly restricts the flexibility of the fitting process is the imposed constraints on the parameter ranges to prevent catastrophic parameter runaway during the fitting.There is scatter in the JWST SCs.We regard that one of the reasons for this is the small number of sources in the bins as well as the faint sources near the detection limit (80% completeness limit).This may be attributed, in part, to cosmic variance.Drawing a conclusive remark on the 21m source count is not straightforward owing to the lack of data points in the brighter flux range, compared to 15 and 18 m counts.An important point of this work is not the difference from previous work, but we used detected Jy sources -just an extrapolation in previous work -to estimate the LFs and CSFH.

CSFH and BHAH
In order to estimate the cosmic evolution of the SF density, we derive the total IR luminosity and IR luminosity functions (IRLFs).Total IR luminosity ( tot 8−1000 ) is derived by integrating galaxy SED from 8m to 1000m, based on Polletta et al. (2007).This contains  (2000).Yellow asterisks (10, 15, and 18 m SCs) and red crosses (18 m SCs) indicate the results on the AKARI's NEP-Deep and NEP-Wide field (Takagi et al. 2012;Pearson et al. 2010Pearson et al. , 2014)).The grey colour indicates that they are excluded in the curve-fitting procedure (see Sec. both SF and AGN contributions.To estimate the SF contribution ( SF 8−1000 = frac × L tot 8−1000 ), we used the fraction (percentage) provided by Gruppioni et al. (2011, Table 3).The SED of SFG and SB comes from purely SF activities.For the composite, AGN2 and AGN1, SF contributions are 96%, 64%, and 46%, respectively.The other parts originate from the AGN activity.Total IR LFs are ob-tained with the updated evolution parameters, following equations (2) and (3).In Fig. 10, we show the total IR LFs of the five populations of galaxies.Based on these new IR LFs, we obtain comoving IR luminosity density,  SF IR (), obtained by multiplying luminosity excluding the AGN contribution to the total IR energy.This can be directly converted to the star formation rate (SFR) density as a function of , using the conversion of Kennicutt (1998),  SFR = 1.7 × 10 −10  SF IR (M ⊙ yr −1 Mpc −3 ) based on the Salpeter initial mass function (Salpeter 1955).This is a crucial tool to understand galaxy evolution in terms of cosmic star formation history (CSFH) (Goto et al. 2010(Goto et al. , 2019)).
We derive the IR luminosity density as a function of redshift, as estimated from our model and compare it with those from other works.In Fig. 11, We show a comparison with a model from Gruppioni et al. (2011, yellow meshed area).The grey dashed curve indicates the cosmic star formation history from Madau & Dickinson (2014).The Blue dotted line indicates the model curve from Harikane et al. (2022), which compares with their high- ( ∼ 13) galaxy candidates.Our model (red curve) shows a rapid increase up to  ∼1, followed by a peak at 1 <  < 2. It decreases sluggishly from  ∼ 2 toward higher-, becoming flat at  > 6.From the current epoch to around the  ∼1, all works show a good agreement.While the pink cross (Barrufet et al. 2023) at  = 5 and the blue dotted line (Harikane et al. 2022) around  = 4.5 show the largest differences (about factor three), our model is consistent with the other recent observational constraints in error ranges.
The black hole accretion rate (BHAR) is described as where BC is the bolometric correction to the 1-1000 m IR luminosity depending on the SED type, and  AGN 1−1000 indicates the 1-1000 m IR luminosity due to the AGN.For the bolometric correction (BC), we take the same values as Gruppioni et al. (2011), and use BC ∼ 1.5 for the AGN1 and AGN2, and BC ∼ 2 for the composite (LLAGN).These bolometric corrections are first-order empirical estimates derived in Pozzi et al. (2007), where the broad-band SEDs of a sample of X-ray-selected AGNs have been studied.Fig. 12 shows the black hole accretion rate derived from our model and comparison with other results.In Merloni & Heinz (2008), the hard X-ray luminosity function is taken as a tracer of the AGN growth rate distribution.They estimated black hole accretion density and presented a synthesis model for AGN evolution based on the hard X-ray LF.Hopkins et al. (2007) derived the quasar bolometric LF from X-ray data.From  = 0 to  = 1, all works are consistent and roughly agree.However, at  > 1.5 differences appear.Compared to other works, our model and Gruppioni et al. (2011) are higher by a few factors.This might be because of the different methods using different data sets.
Our estimates incorporate MIR sources that are several tens of times fainter, thanks to the JWST, which were only predicted or extrapolated in previous studies.As the slight parameter adjustments show, this work demonstrates the remarkable accuracy of predictions and estimations made in the previous studies.

SUMMARY
To fully understand galaxy evolution, it is ideal to have redshift information and investigate galaxy luminosity functions.However, at the time of writing, we do not have photo-z measurements for all the JWST data we used in this work.Some fields lack enough optical data needed for accurate photo-z.Therefore, we took a model-based approach to interpret galaxy number counts and extract information on galaxy evolution.This approach allows us to obtain a physical interpretation of the MIR source counts (Ling et al. 2022, Wu et al. 2023) in a wide flux range using extremely faint galaxies, and constrain cosmic star-formation history and black hole accretion history.Following Gruppioni et al. (2011), we used the backward evolution of the parameterized LLFs for five representative galaxy populations at the JWST/MIRI bands.We fit our model to the source counts from the JWST as well as previous works from ISO, AKARI, and Spitzer (Oliver et al. 1997;Serjeant et al. 2000;Rocca-Volmerange et al. 2007;Pearson et al. 2010Pearson et al. , 2014;;Davidge et al. 2017).By simultaneously fitting the source counts at six mid-infrared bands, we obtained the best-fit evolutions of MIR LFs for each of the five types of galaxies.These parameters gave us inferred CSFH and BHAH.From the current epoch to 1 <  < 2, the obtained CSFH shows good agreement with previous works, at higher-, however, uncertainties become large and the overall decreasing trends are not the same, yet consistent within error ranges.BHAH estimates are consistent with previous estimates in terms of peak density; however, our estimate is slightly higher compared to others at  > 1.These differences seem to be due to the different data and methods used.Our results take advantage of IR data that is less affected by dust extinction.Thanks to the JWST, this work is based on infrared galaxies that are several tens of times fainter than those detectable by previous IR telescopes, shedding light on the accuracy of previous studies.

Figure 1 .
Figure 1.Model spectral energy distribution (SED) templates at  = 0 for five representative galaxy populations from Polletta et al. (2007).Green/orange: normal star-forming galaxy (SFG) population represented by Sa and Sdm with intermediate SEDs (shown in grey) between them.Intermediate SEDs are generated by the combination of the two templates.Cyan/dark red: starburst (SB) population represented by NGC 6090 and Arp 220 with intermediate SEDs (shown in pink).Red: composite galaxy population represented by Seyfert 2 -a mixture of SFG and low luminosity AGN (LLAGN).Magenta: AGN type-2 (obscured), and blue: AGN type-1 (unobscured).In a small box, we display the MIR range (3-40m) along with the JWST/MIRI filter bands.

Figure 2 .
Figure 2. Model local luminosity functions (LLFs) in the six mid-IR (MIR) bands of the JWST/MIRI, for five representative galaxy types, which are derived based on the parameters given for 15m LLF (Table1ofGruppioni et al. 2011).Green: normal star-forming galaxy (SFG) given by spiral, cyan: starburst (SB), red: composite (a mixture of SFG and low luminosity AGN), magenta: type-2 AGN, and blue: type-1 AGN.The thick salmon line indicates the integrated LF of these five types.For 12m and 15 m bands, we compare the model LFs with observed LFs (grey regions) based on AKARI's NEP survey data(Kim et al.  2015).

Figure 3 .
Figure 3.The evolution of the local luminosity functions (LLFs) at 15m band (F1500W), according to equations (2) and (3), with redshift represented by a colour gradient from dark blue through red to yellow (refer to the colour bar on the bottom right).Luminosity bin size is 0.1 dex, and the redshift interval is 0.02, in our computation.Normal star-forming galaxies (SFGs, top-left panel), starbursts (SB, top-right panel), composite (a mixture of SFG and LLAGN, middle-left panel), AGN type-2 (middle-right panel) and AGN type-1 (bottom-left panel) are shown.

Figure 4 .
Figure 4.The evolution of number counts.Here, the x-axis is given in unit of flux density (Jy) at 15 m, which is converted from luminosity (the x-axis in Fig. 3) based on the K-correction using the 15 m filter (F1500W).Each panel shows the evolution of each galaxy type, with redshift represented by a colour bar on the bottom.The broken line in each panel shows integrated values over all redshifts -SFG (top-left panel), starbursts (SB, top-right panel), composite (a mixture of SFG and LLAGN, middle-left panel), AGN type-2 (middle-right panel), and AGN type-1 (bottom left).The bottom right panel shows the sum of all the types.

Figure 5 .
Figure 5.Comparison of LLF evolutions between the initial parameters (shown in Fig.3) and new parameters obtained from the fitting in this work.We show 15m LFs at redshifts  = 1.0, 2.0, 3.0, ..., 7.0, and 8.0 only, for better visibility.The dotted curves indicate the LFs with initial parameters, which are from Figure3.The solid curves are the LFs with newly derived parameters.The inset (a small box) in each panel shows luminosity evolution (solid curve) and density evolution (dot-dashed curve) as a function of redshift.The y-axis 'log 10 (Evol)' represents the evolution terms as a function of , except for the first terms (log 10 L ★ (0) and log 10 Φ ★ (0)) in equations 2 and 3. Therefore, luminosity evolution is log 10 L ★ (z)-log 10 L ★ (0).In this small box, black curves denote evolution with initial parameters, and magenta curves show evolution with new parameters.Therefore, LFs in dotted curves are based on the evolution in black, while LFs in solid curves are based on the magenta curves in the insets.

Figure 6 .
Figure 6.Panel (a): Comparison of 15m LFs -our model LFs at different redshifts (curves), against the observed 15m LFs.Our model LFs are based on the new evolution parameters obtained through the fitting to the SCs.A grey shaded area indicates observed LLF from Kim et al. (2015, AKARI/NEP-Wide field), and the meshed area is from Le Floc'h et al. (2005, Spitzer/CDFS field), and open squares from Ling et al. 2023 (submitted, JWST/CEERS field).Blue, green, orange and red lines indicate the LF at redshift  = 0.5, 1.5, 2.5, and 4.0, respectively.They are the sum of the LFs for each population at each redshift -see panels (b), (c), (d), and (e), which show how the LFs on the top panel are generated.Single LFs for each population in lower panels are a subset of LFs explained in Fig.5.

Figure 7 .
Figure 7.Comparison of the 15 m galaxy number count models at different redshifts between the initial parameters (as given in Fig.4) and the new parameters obtained from the fitting in this work.The dotted curves indicate the number counts models with initial parameters (as given in Fig.4).Solid curves are number counts based on the new parameters.The solid curves represent the source count based on the newly obtained parameters.These curves show slight differences due to the changes in parameters.

Figure 8 .
Figure 8.Comparison of 15 m SC models between the initial and resultant parameters.Individual populations based on the initial parameters (i.e., Gruppioni et al. 2011) are presented by thin lines and new models with updated parameters are presented by thick lines, as specified in the legend.The open red circles represent the recent MIR SC (Ling et al. 2022; Wu et al. 2023) obtained from the early release science (ERS) data of the JWST.Green boxes and blue diamonds indicate the source counts on the AKARI Deep Field South (ADF-S), and Spitzer/IRAC dark field (IRAC-DF) observed by the AKARI, respectively (Davidge et al. 2017).Crosses (+) and yellow asterisks indicate the results from the AKARI's NEP-Deep and NEP-Wide fields(Takagi et al. 2012;Pearson et al. 2010).The extension is indicated by the dashed line from the light brown curve (model fromGruppioni et al. 2011).This extension is created by summing up each extrapolation for five populations.

Figure 9 .
Figure 9. Observed source counts (SCs) from various works and best-fit models at the six bands (F0700W, F1000W, F1280W, F1500W, F1800W, and F2100W) of the JWST/MIRI.Open red circles show the recent MIR SC (Ling et al. 2022, Wu et al. 2023) using early release science (ERS) data of the JWST.The vertical dot-dashed grey line indicates 80% completeness limit.Magenta triangles (in the panels for 10 and 18 m SCs), green boxes (in the panels for 15 m SCs), and blue diamonds (in the panels for 10, 15, and 18 m SCs) show the source counts on the ELAIS field, the AKARI Deep Field South (ADF-S), and Spitzer/IRAC dark field (IRAC-DF) observed by the AKARI, respectively (Davidge et al. 2017).ISO data in 7 and 15 m are from Oliver et al. (1997) and Serjeant et al.(2000).Yellow asterisks (10, 15, and 18 m SCs) and red crosses (18 m SCs) indicate the results on the AKARI's NEP-Deep and NEP-Wide field(Takagi et al. 2012;Pearson et al. 2010Pearson et al. , 2014)).The grey colour indicates that they are excluded in the curve-fitting procedure (see Sec. 3.3).The thick red curve in each panel shows the model fitted to all source count results, which is the sum of all the contributions from five different types of galaxies.The green dashed line indicates the SFG, the cyan dot-dashed line indicates the starburst (SB), the dark-red dot-dashed line indicates the composite (LLAGN), the magenta line indicates the AGN type-2, and the blue dashed line indicates the AGN type-1.
Figure 9. Observed source counts (SCs) from various works and best-fit models at the six bands (F0700W, F1000W, F1280W, F1500W, F1800W, and F2100W) of the JWST/MIRI.Open red circles show the recent MIR SC (Ling et al. 2022, Wu et al. 2023) using early release science (ERS) data of the JWST.The vertical dot-dashed grey line indicates 80% completeness limit.Magenta triangles (in the panels for 10 and 18 m SCs), green boxes (in the panels for 15 m SCs), and blue diamonds (in the panels for 10, 15, and 18 m SCs) show the source counts on the ELAIS field, the AKARI Deep Field South (ADF-S), and Spitzer/IRAC dark field (IRAC-DF) observed by the AKARI, respectively (Davidge et al. 2017).ISO data in 7 and 15 m are from Oliver et al. (1997) and Serjeant et al.(2000).Yellow asterisks (10, 15, and 18 m SCs) and red crosses (18 m SCs) indicate the results on the AKARI's NEP-Deep and NEP-Wide field(Takagi et al. 2012;Pearson et al. 2010Pearson et al. , 2014)).The grey colour indicates that they are excluded in the curve-fitting procedure (see Sec. 3.3).The thick red curve in each panel shows the model fitted to all source count results, which is the sum of all the contributions from five different types of galaxies.The green dashed line indicates the SFG, the cyan dot-dashed line indicates the starburst (SB), the dark-red dot-dashed line indicates the composite (LLAGN), the magenta line indicates the AGN type-2, and the blue dashed line indicates the AGN type-1.

Figure 10 .
Figure 10.Total IR luminosity functions for five galaxy populations at the redshift =1,2,3, ..., 7, which are obtained based on the best-fit parameters to the observed number counts at six MIRI bands.See sec 4.2.
Purple and green shaded areas show the empirical reconstruction of extragalactic background light (EBL) and physical EBL model, respectively (Fermi-LAT Collaboration et al. 2018).We also compare with observed cosmic SFR densities from the literature.The red open circles represent the observed SFR density from Ling et al. (2023, based on the CEERS field IRLFs).Open diamonds indicate the observational constraints from the ALMA (Algera et al. 2023, REBELS).The pink cross indicates the measurement from the HST-dark, JWST/NIRcam galaxies (Barrufet et al. 2023).Open triangles are from Gruppioni et al. (2020, ALPINE-ALMA survey).Navy dashed region is from Magnelli et al. (2013, Herschel-PACS survey on GOODS field).Upward arrows are from Gruppioni et al. (2010), which used dusty galaxies detected by Herschel ( ∼3).Filled blue boxes indicate Rodighiero et al. (2010), based on a combination of Spitzer surveys of VVDS-SWIRE (de la Torre et al. 2007) and GOODS (Giavalisco et al. 2004) fields.Orange-filled boxes are taken from Madau & Dickinson (2014), Finkelstein et al. (2015), McLeod et al. (2016), Bouwens et al. (2020), which indicate the constraints

Figure 12 .
Figure 12.(a) Redshift evolution of the black hole accretion rate (BHAR).Composite, AGN-2, and AGN-1 populations are used whereas SFG and SB are not used to exclude pure SF activity.The red curve indicates our estimation.The red meshed area represents the fitting error reflected on this plane.A small panel on the top-right shows the contribution from three populations.Gruppioni et al. (2011, yellow region) is compared.Open squares indicate the measurements from Yang et al. (2023).The green solid line indicates Ananna et al. (2019, 2020).The Cyan dashed line indicates Ueda et al. (2014) The filled red area with black data point, green hatched area, and blue solid line indicate BHAHs presented by Delvecchio et al. (2014), Merloni & Heinz (2008), and Hopkins et al. (2007), respectively.(b) The overall fraction of the cosmic AGN activity, the ratio of  AGN 1−1000 with respect to the total energy density ( SF 1−1000 +  AGN 1−1000 ). .