Modeling the Extragalactic Background Light and the Cosmic Star Formation History

We use a flat ΛCDM cosmology, where in most of our models we fix the parameters H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m =0.30 (where Ω_Λ = 1 − Ω_m = 0.70 due to the flatness assumption). This cosmology is the most common one used for luminosity density measurements found in the literature, with which we compare our model results, and is thus convenient for our purposes. These cosmological values are also close to those independently measured with a variety of methods. There is a small but statistically significant tension between the value of H₀ found by using measures in the "late" universe (such as Type Ia supernovae and lensed quasars) and in using anisotropy of the cosmic microwave background (CMB) measured by Planck and other experiments (e.g., Riess et al. 2019, 2021, 2022; Wong et al. 2020). We make no attempt to resolve this tension here; however, we note that absorption of γ-rays by EBL photons can be used to constrain the expansion rate of the universe (e.g., Salamon et al. 1994; Mannheim et al. 1996; Blanch & Martinez 2005; Barrau et al. 2008; Domínguez & Prada 2013; Fairbairn 2013; Biteau & Williams 2015; Domínguez et al. 2019; Zeng & Yan 2019). In several of our model fits, we allow H₀ and Ω_m to be free parameters.

We will make frequent use of the cosmological function

$$\begin{eqnarray}&&\displaystyle \frac{{dt}}{{dz}}=\displaystyle \frac{-1}{{H}_{0}(1+z)\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}},\end{eqnarray} \tag{ 2 }$$

which relates a cosmological time interval to a redshift interval.

The "classic" initial mass function (IMF) is that of Salpeter (1955), given by

$$\begin{eqnarray}&&\xi (m)=\displaystyle \frac{{dN}}{{dm}}={N}_{0}{m}^{-2.35}\end{eqnarray} \tag{ 3 }$$

where m is the stellar mass in M_⊙ units. Although still widely used for convenience, this IMF is now disfavored by observations, especially at the low-mass end (e.g., Kroupa 2001; Chabrier 2003). We use it in one of our calculations, primarily to explore the effect of different IMFs on our results.

Baldry & Glazebrook (2003, hereafter BG03) fit a stellar population synthesis model with an IMF that was allowed to vary during their fit to a collection of luminosity density data. We primarily use their best-fit IMF, given by

$$\begin{eqnarray}\xi (m)=\displaystyle \frac{{dN}}{{dm}}={N}_{0}\times \left\{\begin{array}{ll}{\left(m/{m}_{c}\right)}^{-1.5} & m\leqslant {m}_{c}\\ {\left(m/{m}_{c}\right)}^{-2.2} & m\gt {m}_{c},\end{array}\right.\end{eqnarray} \tag{ 4 }$$

where m_c = 0.5.

In both the Salpeter and BG03 cases, the constant N₀ is determined by normalizing the IMF to a total mass of M_⊙, i.e.,

$$\begin{eqnarray}&&1={\int }_{{m}_{\min }}^{{m}_{\max }}{dm}\ m\ \xi (m).\end{eqnarray} \tag{ 5 }$$

In a future publication we will explore models with IMF parameters free to vary in the fit. We use ${m}_{\min }=0.1$ and ${m}_{\max }=120$.

To compute the evolution of the mean metallicity in the interstellar medium (ISM) and the stellar mass density one must know a number of quantities. One is the mass of the stellar remnant (m_r (m, Z)) at the end of a star's life, whether a white dwarf, neutron star, or black hole. This quantity depends on both the progenitor star's mass (m) and metallicity of the gas from which it was born (Z). One must also know the yield of new metals returned to the ISM for a star of a given mass and birth metallicity, p(m, Z). From these we can compute the quantities below.

The fraction of stellar mass returned to the ISM is calculated from (e.g., Madau & Dickinson 2014; Vincenzo et al. 2016)

$$\begin{eqnarray}&&R(Z)={\int }_{{m}_{\mathrm{ret}}}^{{m}_{\max }}{dm}\ [m-{m}_{r}(m,Z)]\ \xi (m),\end{eqnarray} \tag{ 6 }$$

where ξ(m) is the IMF. The yield of heavy elements (the mass of new heavy elements produced) returned to the ISM as a fraction of total mass that is not returned to the ISM is

$$\begin{eqnarray}&&{y}_{\mathrm{tot}}(Z)=\displaystyle \frac{1}{1-R}{\int }_{{m}_{\mathrm{ret}}}^{{m}_{\min }}{dm}\ m\ p(m,Z)\ \xi (m).\end{eqnarray} \tag{ 7 }$$

Here m_ret is the stellar mass above which metals are returned to the ISM. Both R(Z) and y_tot(Z) are fractions returned to the ISM in each stellar generation.

The quantities R(Z) and y_tot(Z) were previously calculated by Vincenzo et al. (2016) using a Salpeter IMF using m_ret = 1.0 and m_r (m, Z) and p(m, Z) as calculated and tabulated by Nomoto et al. (2013). We use these same method to compute R(Z) and y_tot(Z) for the BG03 IMF used here. Our results can be found in Table 1. For values of Z not in Table 1, we did a linear interpolation between these values to determine R(Z) and y_tot(Z). These quantities are only weakly dependent on Z so there should be little uncertainty from this interpolation.

To determine the ISM mean metallicity $\bar{Z}(z)$ and comoving stellar mass density ρ(z) of the universe, we use a one-zone, closed box, instantaneous recycling model (Tinsley 1980). This model assumes that stars with masses m < m_ret will live forever and never return matter to the ISM, and stars with masses m > m_ret will instantly return mass and newly formed metals to the ISM (e.g., Madau & Dickinson 2014; Vincenzo et al. 2016). This model is a good approximation for metals produced primarily by massive, short-lived stars (such as oxygen) but less good for metals produced primarily by low-mass, long-lived stars. Since oxygen is the most abundant heavy element by mass, this model should be a reasonable approximation for mean metallicity (Vincenzo et al. 2016).

The quantities $\bar{Z}(z)$ and ρ(z) are found by solving the equations (e.g., Madau & Dickinson 2014)

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\bar{Z}}{{dz}} & = & \displaystyle \frac{{y}_{\mathrm{tot}}(\bar{Z})[1-R(\bar{Z})]}{{\rho }_{b}-\rho (z)}\psi (z)\left|\displaystyle \frac{{dt}}{{dz}}\right|,\\ \displaystyle \frac{d\rho }{{dz}} & = & (1-R(\bar{Z}))\psi (z)\left|\displaystyle \frac{{dt}}{{dz}}\right|,\end{array}\end{eqnarray} \tag{ 8 }$$

respectively, where $R(\bar{Z})$ and y_tot are described in Section 2.3, ψ(z) is the comoving SFRD and ρ_b is the total comoving mass density of baryons (in stars and gas) in the universe. It can be found from ρ_b = ρ_c Ω_b where

$$\begin{eqnarray}&&{\rho }_{c}=\displaystyle \frac{3{H}_{0}}{8\pi G},\end{eqnarray} \tag{ 9 }$$

and G = 6.673 × 10⁻⁸ dyn cm² g⁻² = 6.673 × 10⁻¹¹ N m² kg⁻² is the gravitational constant. The parameter Ω_b can be determined from anisotropy in the CMB. We use Ω_b = 0.045 consistent with results from the Wilkinson Microwave Anisotropy Probe (Hinshaw et al. 2013).

The coupled ordinary differential Equations (8) were solved with a fourth-order Runge–Kutta numerical scheme (e.g., Press et al. 1992) with the initial conditions $\bar{Z}(z={z}_{\max })=0$ and $\rho (z={z}_{\max })=0$.

We use the PEGASE.2 code ⁸ (Fioc & Rocca-Volmerange1997) to generate simple stellar population spectra (SSPS) of a population of stars with masses between m_min and m_max as a function of metallicity (Z) and age (t). This assumes all stars are born instantly at the same time at t = 0, and evolve passively with no further star formation. The SSPS L_λ (t, Z) are computed in units ergs⁻¹ Å⁻¹ and are normalized to 1 M_⊙. The PEGASE.2 code allows the user to generate SSPS for user-defined IMFs. As discussed in Section 2.2, we use the BG03 IMF, Equation (4), and the Salpeter IMF, Equation (3). The primary effect of metallicity on the SSPS is that at higher metallicity, there is more IR emission and less UV emission (Figure 1).

Zoom In Zoom Out Reset image size

The PEGASE.2 code does not produce accurate SSPS for high-mass Population III (low-Z) stars. Therefore, for SSPS at Z = 0, we follow Gilmore (2012) and use the results from Schaerer (2002) for high-mass stars. For Z = 0 stars with m < 5, we use the PEGASE.2 results; for m > 5, we assume the stars are blackbodies with temperatures and luminosities given by Table 3 of Schaerer (2002) for lifetimes given by Table 4 of Schaerer (2002).

Once SSPS are produced for a grid of stellar population ages and metallicities, we interpolate to find the SSPS L_λ (t, Z) for any t and Z.

A fraction f_esc,d (λ, z) of photons will escape absorption by dust. Driver et al. (2008) calculated the average of this fraction at z ≈ 0 from fitting a dust model to 10,000 nearby galaxies. This fraction at z ≈ 0 was fit by Razzaque et al. (2009), resulting in

$$\begin{eqnarray}&&\begin{array}{l}{f}_{\mathrm{esc},d}(\lambda ,z=0)\\ =\,\left\{\begin{array}{ll}0.688+0.556{\mathrm{log}}_{10}\lambda ; & \lambda \leqslant 0.167\\ 0.151-0.136{\mathrm{log}}_{10}\lambda ; & 0.167\lt \lambda \leqslant 0.218\\ 1.000+1.148{\mathrm{log}}_{10}\lambda ; & 0.22\lt \lambda \leqslant 0.422\\ 0.728+0.422{\mathrm{log}}_{10}\lambda ; & 0.422\lt \lambda \end{array}\right.,\end{array}\end{eqnarray} \tag{ 10 }$$

where λ is wavelength in μm. The dust extinction evolves with redshift; Abdollahi et al. (2018) fit a collection of measurements with a curve as a function of z. We assume the normalization of f_esc,d (λ, z) evolves with z following this curve, so that

$$\begin{eqnarray}&&{f}_{\mathrm{esc},d}(\lambda ,z)=\displaystyle \frac{{f}_{\mathrm{esc},d}(\lambda ,z=0)}{{f}_{\mathrm{esc},d}(\lambda =0.15\ \mu {\rm{m}},z=0)}\times {10}^{-0.4A(z)},\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&A(z)={m}_{d}\displaystyle \frac{{\left(1+z\right)}^{{n}_{d}}}{1+{\left[(1+z)/{p}_{d}\right]}^{{q}_{d}}},\end{eqnarray} \tag{ 12 }$$

with the Abdollahi et al. (2018) fit results m_d = 1.49, n_d = 0.64, p_d = 3.4, q_d = 3.54. A similar fit was done by Puchwein et al. (2019) using a slightly different functional form.

We also use a dust model where

$$\begin{eqnarray}\begin{array}{ll}{f}_{\mathrm{esc},d}(\lambda ,z=0) & =\,\left\{\begin{array}{ll}{f}_{\mathrm{esc},2}+({f}_{\mathrm{esc},2}-{f}_{\mathrm{esc},1})/({\mathrm{log}}_{10}{\lambda }_{2}-{\mathrm{log}}_{10}{\lambda }_{1}){\mathrm{log}}_{10}\lambda ; & \lambda \leqslant {\lambda }_{2}\\ {f}_{\mathrm{esc},3}+({f}_{\mathrm{esc},3}-{f}_{\mathrm{esc},2})/({\mathrm{log}}_{10}{\lambda }_{3}-{\mathrm{log}}_{10}{\lambda }_{2}){\mathrm{log}}_{10}\lambda ; & {\lambda }_{2}\lt \lambda \leqslant {\lambda }_{3}\\ {f}_{\mathrm{esc},4}+({f}_{\mathrm{esc},4}-{f}_{\mathrm{esc},3})/({\mathrm{log}}_{10}{\lambda }_{4}-{\mathrm{log}}_{10}{\lambda }_{3}){\mathrm{log}}_{10}\lambda ; & {\lambda }_{3}\lt \lambda \leqslant {\lambda }_{4}\\ {f}_{\mathrm{esc},5}+({f}_{\mathrm{esc},5}-{f}_{\mathrm{esc},4})/({\mathrm{log}}_{10}{\lambda }_{5}-{\mathrm{log}}_{10}{\lambda }_{4}){\mathrm{log}}_{10}\lambda ; & {\lambda }_{4}\lt \lambda \end{array}\right.\end{array}\end{eqnarray} \tag{ 13 }$$

where λ₁ = 0.15 μm, λ₂ = 0.167 μm, λ₃ = 0.218 μm, λ₄ = 0.422 μm, λ₅ = 2.0 μm. Here {f_esc,k , k = 1...5} are free parameters.

We use two parameterizations for the comoving SFRD as a function of redshift. We primarily use the formulation from Madau & Dickinson (2014):

$$\begin{eqnarray}&&{\psi }_{\mathrm{MD}}(z)={a}_{s}\displaystyle \frac{{\left(1+z\right)}^{{b}_{s}}}{1+{\left[(1+z)/{c}_{s}\right]}^{{d}_{s}}},\end{eqnarray} \tag{ 14 }$$

where a_s , b_s , c_s , and d_s are free parameters. We also use a version similar to the piecewise function from Hopkins & Beacom (2006),

$$\begin{eqnarray}{\psi }_{\mathrm{piece}}(z)={10}^{{a}_{s}}\left\{\begin{array}{ll}{\zeta }^{{b}_{s}} & z\lt {z}_{1}\\ {\zeta }_{1}^{{b}_{s}-{c}_{s}}{\zeta }^{{c}_{s}} & {z}_{1}\leqslant z\lt {z}_{2}\\ {\zeta }_{1}^{{b}_{s}-{c}_{s}}{\zeta }_{2}^{{c}_{s}-{d}_{s}}{\zeta }^{{d}_{s}} & {z}_{2}\leqslant z\lt {z}_{3}\\ {\zeta }_{1}^{{b}_{s}-{c}_{s}}{\zeta }_{2}^{{c}_{s}-{d}_{s}}{\zeta }_{3}^{{d}_{s}-{e}_{s}}{\zeta }^{{e}_{s}} & {z}_{3}\leqslant z\lt {z}_{4}\\ {\zeta }_{1}^{{b}_{s}-{c}_{s}}{\zeta }_{2}^{{c}_{s}-{d}_{s}}{\zeta }_{3}^{{d}_{s}-{e}_{s}}{\zeta }_{4}^{{e}_{s}-{f}_{s}}{\zeta }^{{f}_{s}} & {z}_{3}\leqslant z,\end{array}\right.\end{eqnarray} \tag{ 15 }$$

where a_s , b_s , c_s , d_s , e_s , and f_s are free parameters, and we use the notation ζ = (1 + z) and {ζ_k = (1 + z_k ), k = 1...4}. We fix z₁ = 1.0, z₂ = 2.0, z₃ = 3.0, z₄ = 4.0.

The comoving stellar luminosity density (i.e., luminosity per unit comoving volume) is given as a function of comoving dimensionless energy = hc/(λ m_e c²) by

$$\begin{eqnarray}\begin{array}{rcl}\epsilon {j}^{\mathrm{stars}}(\epsilon ;z) & = & {m}_{e}{c}^{2}{\epsilon }^{2}\displaystyle \frac{{dN}}{{dtd}\epsilon {dV}}\\ & = & {f}_{\mathrm{esc},{\rm{H}}}(\epsilon ){f}_{\mathrm{esc},d}(\lambda ,z)\\ & & \times {\displaystyle \int }_{z}^{{z}_{\max }}{{dz}}_{1}\ {L}_{\star }({t}_{\star }(z,{z}_{1}),Z({z}_{1}))\ \psi ({z}_{1})\left|\displaystyle \frac{{dt}}{{{dz}}_{1}}\right|,\end{array}\end{eqnarray} \tag{ 16 }$$

where L_⋆(t, Z) = λ × L_λ (t, Z). We assume all photons at energies greater than 13.6 eV are absorbed by H i gas, so that

$$\begin{eqnarray}{f}_{\mathrm{esc},{\rm{H}}}(\epsilon )=\left\{\begin{array}{ll}1 & {m}_{e}{c}^{2}\epsilon \lt 13.6\ \mathrm{eV}\\ 0 & {m}_{e}{c}^{2}\epsilon \geqslant 13.6\ \mathrm{eV}\end{array}\right..\end{eqnarray} \tag{ 17 }$$

The age of the stellar population t_⋆(z, z₁) can be found by computing the integral

$$\begin{eqnarray}&&{t}_{\star }(z,{z}_{1})={\int }_{z}^{{z}_{1}}{{dz}}^{{\prime} }\left|\displaystyle \frac{{dt}}{{{dz}}^{{\prime} }}\right|,\end{eqnarray} \tag{ 18 }$$

which has an analytic solution (Razzaque et al. 2009).

We assume that the total energy absorbed by dust is re-emitted in the IR in three blackbody dust components. These three components are a large-grain component with temperature T₁ (left as a free parameter), a small-grain component with temperature T₂ = 70 K, and a component representing polycyclic aromatic hydrocarbons as a blackbody with temperature T₃ = 450 K. The comoving dust luminosity density is then

$$\begin{eqnarray}\displaystyle \begin{array}{rcl}\epsilon {j}^{\mathrm{dust}}(\epsilon ;z) & = & \displaystyle \frac{15}{{\pi }^{4}}\displaystyle \int d\epsilon \left[\displaystyle \frac{1}{{f}_{\mathrm{esc},d}(\lambda ,z)}-1\right]{j}^{\mathrm{stars}}(\epsilon ;z)\\ & & \times \sum _{n=1}^{3}\displaystyle \frac{{f}_{n}}{{{\rm{\Theta }}}_{n}^{4}}\displaystyle \frac{{\epsilon }^{4}}{\exp [\epsilon /{{\rm{\Theta }}}_{n}]-1},\end{array}\end{eqnarray} \tag{ 19 }$$

where Θ_n = k_B T_n /(m_e c²) is the dimensionless temperature of a particular component and f_n < 1 is the fraction of the emission absorbed by dust that is reradiated in a particular dust component. We use f₁ and f₂ as free parameters, with f₃ constrained by f₁ + f₂ + f₃ = 1. The temperatures of each component are fixed to the values given above.

Once the stellar and dust luminosity densities are calculated, the total luminosity density is simply

$$\begin{eqnarray}&&\epsilon {j}^{\mathrm{tot}}(\epsilon ;z)=\epsilon {j}^{\mathrm{stars}}(\epsilon ;z)+\epsilon {j}^{\mathrm{dust}}(\epsilon ;z).\end{eqnarray} \tag{ 20 }$$

The proper frame energy density as a function of the proper frame dimensionless energy _p is

$$\begin{eqnarray}&&{\epsilon }_{p}{u}_{\mathrm{EBL},{\rm{p}}}({\epsilon }_{p};z)={\left(1+z\right)}^{4}{\int }_{z}^{{z}_{\max }}{{dz}}_{1}\displaystyle \frac{{\epsilon }^{{\prime} }{j}^{\mathrm{tot}}({\epsilon }^{{\prime} };{z}_{1})}{(1+{z}_{1})}\left|\displaystyle \frac{{dt}}{{dz}}\right|,\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&{\epsilon }^{{\prime} }=\displaystyle \frac{1+{z}_{1}}{1+z}{\epsilon }_{p}.\end{eqnarray} \tag{ 22 }$$

EBL intensity is usually measured in units of nW m⁻² sr⁻¹; the energy density can be converted to intensity via

$$\begin{eqnarray}&&{\epsilon }_{p}{I}_{\mathrm{EBL}}({\epsilon }_{p};z)=\displaystyle \frac{c}{4\pi }{\epsilon }_{p}{u}_{\mathrm{EBL},{\rm{p}}}({\epsilon }_{p};z).\end{eqnarray} \tag{ 23 }$$

The absorption optical depth is

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\gamma \gamma }({\epsilon }_{1},z) & = & \displaystyle \frac{c\pi {r}_{e}^{2}}{{\epsilon }_{1}^{2}{m}_{e}{c}^{2}}\ {\displaystyle \int }_{0}^{z}\displaystyle \frac{{{dz}}^{{\prime} }}{{\left(1+{z}^{{\prime} }\right)}^{2}}\ \left|\displaystyle \frac{{{dt}}_{* }}{{{dz}}^{{\prime} }}\right|\ \\ & & \times {\displaystyle \int }_{\tfrac{1}{{\epsilon }_{1}(1+{z}^{{\prime} })}}^{\infty }d{\epsilon }_{p}\displaystyle \frac{{\epsilon }_{p}{u}_{\mathrm{EBL},p}({\epsilon }_{p};{z}^{{\prime} })}{{\epsilon }_{p}^{4}}\ \bar{\phi }({\epsilon }_{p}{\epsilon }_{1}(1+{z}^{{\prime} })),\end{array}\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&\begin{array}{l}\bar{\phi }({s}_{0})=\displaystyle \frac{1+{\beta }_{0}^{2}}{1-{\beta }_{0}^{2}}\mathrm{ln}{w}_{0}-{\beta }_{0}^{2}\mathrm{ln}{w}_{0}-\displaystyle \frac{4{\beta }_{0}}{1-{\beta }_{0}^{2}}\\ +2{\beta }_{0}+4\mathrm{ln}{w}_{0}\mathrm{ln}(1+{w}_{0})-4L({w}_{0}),\end{array}\end{eqnarray} \tag{ 25 }$$

${\beta }_{0}^{2}=1-1/{s}_{0}$, w₀ = (1 + β₀)/(1 − β₀), and

$$\begin{eqnarray}&&L({w}_{0})={\int }_{1}^{{w}_{0}}dw\,{w}^{-1}{\rm{l}}{\rm{n}}(1+w)\end{eqnarray} \tag{ 26 }$$

(Gould & Schréder 1967; Brown et al. 1973).

Ackermann et al. (2012) developed a method to determine the EBL absorption optical depth from γ-ray observations of a large sample of blazars. This involved a joint likelihood fit to a large statistical sample of blazars, using the LAT spectrum in the 1.0–10 GeV (where the γ-ray absorption should be negligible) extrapolated to >10 GeV as the intrinsic spectrum of each source, and allowing the normalization of the absorption optical depth to be an additional free parameter for all sources within a redshift and energy bin. They applied their method to 150 blazars in 3 redshift bins, using ≈4 yr of data. Abdollahi et al. (2018) applied this technique to a much larger sample and much more data: 9 yr of LAT data from 739 blazars, now divided into 12 redshift bins. They also included a constraint from the LAT detection of GRB 080916C (see also Desai et al. 2017). Desai et al. (2019) applied this method to measure the absorption optical depth using 106 IACT VHE γ-ray spectra from 38 blazars from the sample compiled by Biteau & Williams (2015). They divided their sample into two redshift bins. We include in our fit here the absorption optical depth measured with the 739 blazars and one GRB measured by Abdollahi et al. (2018), and measured with the 106 VHE spectra by Desai et al. (2019). The reported errors for all of the γ-ray absorption data include both statistical and systematic uncertainties.

Generally, γ-ray telescopes are only able to measure values of the EBL absorption optical depth in the range 10⁻² ≲ τ_{γ γ} ≲ 5. This is because, for low values of τ_{γ γ} , e.g., exp(−0.01) ≈ 0.99, absorption will be negligible, and there will be nothing to measure. At high values of τ_{γ γ} , e.g., exp(−5) ≈ 0.007, almost all of the γ-rays will be absorbed, and there will be no γ-ray signal to measure. Outside of this range, however, it is possible to put upper limits on τ_{γ γ} at lower optical depths, and lower limits on τ_{γ γ} at higher optical depths.

One should also note that Lorentz invariance violation could modify the γ γ pair-production threshold, reducing the absorption optical depth to γ-rays at ≳10 TeV (e.g., Kifune 1999). In this case, γ-rays at these energies could be potentially detectable with the upcoming Cherenkov Telescope Array (CTA; Abdalla & Böttcher 2018; Abdalla et al. 2021). If axion-like particles exist, γ-ray photons could convert to them and in principle avoid γ γ pair production (de Angelis et al. 2007; Sánchez-Conde et al. 2009). However, there is no evidence that this is occurring (Buehler et al. 2020).

We have searched through the literature to compile a large amount of luminosity density measurements from galaxy surveys to fit with our EBL model. We particularly found useful the compilations by Helgason & Kashlinsky (2012), Stecker et al. (2012, 2016), Scully et al. (2014), and Madau & Dickinson (2014). Our luminosity density data compiled from the literature can be found in Table 2. The redshift and wavelength coverage can be seen in Figure 2. The majority of the data come from three sources: Tresse et al. (2007), Andrews et al. (2017), and Saldana-Lopez et al. (2021). Coverage is fairly complete at z ≲ 1, mainly due to the work of Andrews et al. (2017) using data from the Galaxy and Mass Assembly and Cosmic Origins Survey projects. They included data from a wide variety of instruments across the electromagnetic spectrum. There is less complete coverage at longer wavelengths moving to higher z. Saldana-Lopez et al. (2021) made use of the CANDELS programs to get fairly complete wavelength coverage of the luminosity density up to z = 6. At z ≳ 6, only UV luminosity density data are available. These high-z UV data are due to deep exposures with the WFC3 instrument on the Hubble Space Telescope (HST; Finkelstein et al. 2015; Bouwens et al. 2016; McLeod et al. 2016). However, see Oesch et al. (2018). Much of the luminosity density data we use here were published after our last model (Finke et al. 2010), giving the work here a different view of the luminosity density, particularly at high redshift. Occasionally we needed to convert luminosity density from L_⊙ Mpc⁻³ to W Mpc⁻³, which we did using L_⊙ = 3.826 × 10²⁶ W.

Zoom In Zoom Out Reset image size

Table 2. Luminosity Density Data Used in Fits

z	λ (μm)	j(; z) (W Mpc−3)	Reference
0.0	0.1535	${0.84}_{-0.07}^{+0.07}\times {10}^{34}$	a
0.0	0.2301	${0.84}_{-0.07}^{+0.07}\times {10}^{34}$	a
0.0	0.3557	${1.47}_{-0.07}^{+0.07}\times {10}^{34}$	a
0.0	0.4702	${3.29}_{-0.07}^{+0.07}\times {10}^{34}$	a
0.0	0.6175	${4.41}_{-0.14}^{+0.14}\times {10}^{34}$	a
0.0	0.7491	${4.62}_{-0.14}^{+0.14}\times {10}^{34}$	a
0.0	0.8946	${4.69}_{-0.14}^{+0.14}\times {10}^{34}$	a
0.0	1.0305	${4.20}_{-0.21}^{+0.21}\times {10}^{34}$	a
0.0	1.2354	${3.78}_{-0.21}^{+0.21}\times {10}^{34}$	a
0.0	1.6458	${3.78}_{-0.21}^{+0.21}\times {10}^{34}$	a

Notes.^a Driver et al. (2012). ^b Andrews et al. (2017). ^c Tresse et al. (2007). ^d Arnouts et al. (2007). ^e Beare et al. (2019). ^f Budavári et al. (2005). ^g Cucciati et al. (2012). ^h Dahlen et al. (2007). ⁱ Marchesini et al. (2012). ^j Schiminovich et al. (2005). ^k Stefanon & Marchesini (2013). ^l Bouwens et al. (2016). ^m Finkelstein et al. (2015). ⁿ McLeod et al. (2016). ^o Dai et al. (2009). ^p Babbedge et al. (2006). ^q Takeuchi et al. (2006). ^r Saldana-Lopez et al. (2021). ^s Yoshida et al. (2006).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

When necessary we converted the luminosity density measurements to a cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 0.30 assumed in most of our models (Section 2.1). In several of our model fits (Models A.c and D.c), we allow these cosmological parameters to be free parameters. In these cases, the model luminosity density is converted from the model cosmology to the "standard" (H₀, Ω_m ) = (70 km s⁻¹ Mpc⁻¹, 0.30) cosmology used in the measurements with

$$\begin{eqnarray}&&j{\left(\epsilon ,z\right)}_{\mathrm{std}}=j{\left(\epsilon ,z\right)}_{\mathrm{model}}\times \displaystyle \frac{{\left|{dt}/{dz}\right|}_{\mathrm{model}}}{{\left|{dt}/{dz}\right|}_{\mathrm{std}}}.\end{eqnarray} \tag{ 27 }$$

Madau & Dickinson (2014) have compiled a large number of cosmological stellar mass density measurements. Stellar mass density measurements are model-dependent, particularly on the IMF. Madau & Dickinson have scaled their mass density calculation to a Salpeter IMF.

We use stellar mass density data from the compilation by Madau & Dickinson (2014). For models that use the BG03 IMF, we adjusted the stellar mass density as follows:

$$\begin{eqnarray}&&{\rho }_{\mathrm{BG}03}={\rho }_{\mathrm{Sal}}\times \displaystyle \frac{1-{R}_{\mathrm{BG}03}}{1-{R}_{\mathrm{Sal}}},\end{eqnarray} \tag{ 28 }$$

where ρ_BG03 and ρ_Sal are the mass density data with the BG03 and Salpeter IMFs, respectively, and R_BG03 and R_Sal are the return fractions (Equation (6)) for the BG03 and Salpeter IMFs, respectively. We use R_BG03 = 0.49 and R_Sal = 0.3 for all z and $\bar{Z}$; although the actual values are metallicity-dependent, these are roughly the midpoint of the different values (Table 1). To take into account the errors from assuming all stars have the same R, we increased the error on all ρ_BG03 by 10%, added in quadrature; that is, the resulting error in ρ_BG03 is

$$\begin{eqnarray}&&{\sigma }_{\mathrm{BG}03}=\sqrt{{\sigma }_{\mathrm{Sal}}^{2}+{\left(0.1{\rho }_{\mathrm{BG}03}\right)}^{2}},\end{eqnarray} \tag{ 29 }$$

where σ_Sal is the error on ρ_Sal reported by Madau & Dickinson (2014).

Similar to the luminosity density, when the cosmological parameters are allowed to vary in the model, we convert the model stellar mass density to the standard cosmology used in the stellar mass density measurements with

$$\begin{eqnarray}&&\rho {\left(z\right)}_{\mathrm{std}}=\rho {\left(z\right)}_{\mathrm{model}}\times \displaystyle \frac{{\left|{dt}/{dz}\right|}_{\mathrm{model}}}{{\left|{dt}/{dz}\right|}_{\mathrm{std}}}.\end{eqnarray} \tag{ 30 }$$

We include in our model fit the z = 0 EBL intensity-integrated galaxy light lower limits from Driver et al. (2016). This integrated galaxy light represents the resolved component of the EBL. The lower limits span a wavelength range 0.15 μm < λ < 500 μm and were derived from galaxy number count data from a variety of ground- and space-based sources.

We use the far-UV (FUV; 0.15 μm) dust extinction as measured by Burgarella et al. (2013) and Andrews et al. (2017). Burgarella et al. derived the extinction from the ratio of the FUV to integrated IR luminosity functions. Andrews et al. determined these FUV extinction values by fitting a dust model to a variety of multiwavelength data.

We use the open-source python-implemented MCMC code emcee (Foreman-Mackey et al. 2013) to determine our best-fit model parameters and their posterior probabilities. This code implements the affine-invariant ensemble sampler of Goodman & Weare (2010). It generates several MCMC chains ("walkers") that evolve through the model space simultaneously, with different model parameters (represented by $\overrightarrow{\theta }$ below) each time step. We typically use 86 walkers here. We do not use the first ≈2000 states for each walker (slightly different for each model) to avoid the "burn-in" period.

The MCMC algorithm uses a likelihood function, ${ \mathcal L }\,\propto \exp (-{\chi }^{2})$. For most of our models, we include fits to the γ-ray opacity (τ), luminosity density (j), stellar mass density (ρ), EBL intensity (I_EBL), and FUV dust extinction (A_FUV) data described above. Thus,

$$\begin{eqnarray}\begin{array}{rcl}{\chi }^{2} & = & \displaystyle \sum _{i}\displaystyle \frac{{\left[{\tau }_{\mathrm{obs}}({E}_{i},{z}_{i})-{\tau }_{\mathrm{model}}(\overrightarrow{\theta }| {E}_{i},{z}_{i})\right]}^{2}}{{\sigma }_{\tau ,i}^{2}}\\ & & +\displaystyle \sum _{i}\displaystyle \frac{{\left[{j}_{\mathrm{obs}}({\lambda }_{i},{z}_{i})-{j}_{\mathrm{model}}(\overrightarrow{\theta }| {\lambda }_{i},{z}_{i})\right]}^{2}}{{\sigma }_{j,i}^{2}}\\ & & +\displaystyle \sum _{i}\displaystyle \frac{{\left[{\rho }_{\mathrm{obs}}({z}_{i})-{\rho }_{\mathrm{model}}(\overrightarrow{\theta }| {z}_{i})\right]}^{2}}{{\sigma }_{\rho ,i}^{2}}\\ & & +\displaystyle \sum _{i}\displaystyle \frac{{\left[{I}_{\mathrm{EBL},\mathrm{obs}}({\lambda }_{i},{z}_{i})-{I}_{\mathrm{EBL},\mathrm{model}}(\overrightarrow{\theta }| {\lambda }_{i},{z}_{i})\right]}^{2}}{{\sigma }_{{I}_{\mathrm{EBL}},i}^{2}}\\ & & +\displaystyle \sum _{i}\displaystyle \frac{{\left[{A}_{\mathrm{FUV},\mathrm{obs}}({z}_{i})-{A}_{\mathrm{FUV},\mathrm{model}}(\overrightarrow{\theta }| {z}_{i})\right]}^{2}}{{\sigma }_{{A}_{\mathrm{FUV}},i}^{2}}.\end{array}\end{eqnarray} \tag{ 31 }$$

In some cases the upper and lower errors differ. In these cases, if ${\tau }_{\mathrm{obs}}({E}_{i},{z}_{i})\gt {\tau }_{\mathrm{model}}(\overrightarrow{\theta }| {E}_{i},{z}_{i})$ we use the lower error for σ_τ,i ; otherwise we use the upper error. We similarly discriminate between upper and lower errors for the other data points.

We use flat priors for all our parameters with appropriately large limits. Our results are generally independent of the limits on the priors except for the fit to the τ_{γ γ} only (Model B) as described by Abdollahi et al. (2018).

We consider Model A to be our fiducial model. This model uses the Madau & Dickinson (2014) SFRD parameterization, Equation (14), and the Baldry & Glazebrook (2003) IMF, Equation (4), and has all of the SFRD and dust extinction and emission parameters free in the fit, but keeps the cosmological parameters fixed. It fits all of the data described in Section 3: the γ-ray opacity, luminosity density, stellar mass density, and FUV dust extinction data, as well as the EBL intensity lower limits.

The luminosity densities as a function of z from our model are plotted in Figure 3, and compared with data. The model is clearly a good fit to the most recent data at these wavelengths. It is interesting to note the divergence at high z with the older data and model. In the FUV band, the 0.17 μm older data from Sawicki & Thompson (2006) used in the previous model by Finke et al. (2010) is plotted. Sawicki & Thompson (2006) used deep imaging with the Keck telescope. The updated deep HST observations (Bouwens et al. 2016) have found much more light at these wavelengths. Since the model from Finke et al. (2010) was designed to reproduce these observations (among others), it is not surprising that it is lower than Model A of this work. One should note that the models from Finke et al. (2010) did not extend beyond z = 6. At 0.28 μm and 0.44 μm, the older model, current model, and all data are in reasonably good agreement.

Zoom In Zoom Out Reset image size

As one moves to longer wavelengths, our Model A becomes lower than the older model from Finke et al. (2010), particularly at z > 2. The model reproduces most of the data at most wavelengths shown here. For instance, the current Model A follows the data from Stefanon & Marchesini (2013) at 1.22 μm, which was obviously not available to Finke et al. (2010). However, our Model A is considerably below the luminosity density measurements from Saldana-Lopez et al. (2021) at 2.2 μm, particularly at z ≳ 2, despite the fact that these data were included in the fit. We think this is because of our poor modeling of the polycyclic aromatic hydrocarbon (PAH) component (see below). The model is generally consistent with the Saldana-Lopez et al. (2021) data at z ≲ 2 where the emission is dominated by stellar emission, rather than PAH emission. The model also closely follows the lower luminosity densities of Stefanon & Marchesini (2013) at 1.2 μm and 1.6 μm at z ≲ 2. The stellar models seem to be highly constraining in this wavelength range. However, this is the wavelength range identified by Conroy et al. (2009) where assumptions about all stars having the same metallicity, rather than a distribution of metallicities, can create uncertainty in emission. This is illustrated in Figure 1 where different metallicities can create significant differences at the redder wavelengths. At the longer IR wavelengths (3.6–8.0 μm), most of the data are from the work of Babbedge et al. (2006) and Dai et al. (2009), and our model fits them well. However, at these wavelengths, the data do not extend to very high redshift, and the data that do exist have large uncertainties; thus this wavelength range is not strongly constrained. This paucity of data should be rectified in the near future with the James Webb Space Telescope Mid-Infrared Instrument. At 3.6–8.0 μm and z ≳ 2 our model is considerably lower than the older model from Finke et al. (2010).

The luminosity density as a function of wavelength at four redshifts is shown in Figure 4. At low redshifts (z ≤ 0.5) the model follows the data from Andrews et al. (2017) closely. At z = 0.10, older measurements at 3.6–8 μm by Babbedge et al. (2006) are generally higher but also have greater uncertainty compared with more recent measurements by Andrews et al. (2017; see also large uncertainties in these data in Figure 3). At this redshift, our updated Model A here and the model from Finke et al. (2010) are in agreement at all wavelengths except the λ ≈ 3.0–20 μm range, mainly due to the influence of the Babbedge et al. (2006) data on the older model. For the ≥12 μm measurements by Takeuchi et al. (2006), the measurements are somewhat lower than the model and other measurements, but still mostly consistent with them, considering the larger uncertainties. At z = 0.5, essentially all the data and models are in agreement.

Zoom In Zoom Out Reset image size

At z = 0.9, some separation between models and data can be seen. At this redshift, at λ < 1.2 μm, the measurements from Tresse et al. (2007) are below those from Andrews et al. (2017). Our Model A follows the Tresse et al. (2007) measurements here more closely, which is quite a bit lower than Model C from Finke et al. (2010) in the λ = 0.8–3.0 μm range. At λ > 30 μm, the updated model is higher than the older model from Finke et al. (2010). This is due to the influence of the Andrews et al. (2017) data at these wavelengths. This redshift (z = 0.9) is the highest redshift for which Andrews et al. have measurements.

The luminosity density at z = 3 is a particularly interesting case. The older model from Finke et al. (2010) is quite divergent from our updated model here. In the model itself, this seems to be related to the updated stellar spectra. The data we fit are also influential, and the updated model reproduces the data at λ < 3 μm well. Most of the optical-near-IR measurements used at this redshift (Marchesini et al. 2012; Stefanon & Marchesini 2013; Saldana-Lopez et al. 2021) did not exist when Finke et al. (2010) created their model. The only data point at z = 3 that is not reproduced by the model is the λ = 2.2 μm measurement by Saldana-Lopez et al. (2021). This could be due to the approximation of the PAH component by a blackbody. A more realistic treatment of this component might reproduce the high-redshift 2.2 μm data better, where there is considerable disagreement between data and model (see Figure 3). Another possibility is that our model is underestimating the metallicity (see Figure 8 below), and that a higher metallicity would lead to more intense luminosity density at this wavelength (Figure 1). Finally, it could be that AGNs, neglected here, could contribute to the 2.2 μm luminosity density observed by Saldana-Lopez et al. (2021). The AGNs would be required to make up 70% of the 2.2 μm luminosity density at z = 3.75, and 90% at z = 5.75. Estimates of the AGN contribution to the EBL indicate that it likely makes up no more than 10% at any relevant wavelength (Domínguez et al. 2011; Abdollahi et al. 2018; Andrews et al. 2018; Khaire & Srianand 2019). We therefore believe it unlikely that AGNs could be the reason for the discrepancy here.

The 68% confidence interval for the EBL intensity at z = 0 for our Model A is shown in Figure 5, along with a variety of lower limits and other models from the literature. Our model result has little uncertainty at λ ≲ 10 μm, while the uncertainty grows at longer wavelengths. Our updated model's EBL intensity is lower than most other models at λ ≲ 10 μm, and quite close to being in conflict with some of the lower limits from Driver et al. (2016), and below those from Biteau & Williams (2015). The only lower EBL model in this wavelength range is the one from Kneiske & Dole (2010), which was designed to be the lowest possible model consistent with lower limits. At longer wavelengths, uncertainty is greater. The far-IR peak at ≈200 μm is greater than the older model from Finke et al. (2010), considerably lower than the one from Domínguez et al. (2011), and mostly consistent with the others.

Zoom In Zoom Out Reset image size

The EBL absorption optical depth from our Model A for four redshift ranges is shown in Figure 6. These are compared with the EBL absorption measurements with the Fermi-LAT (Abdollahi et al. 2018) and IACTs (Desai et al. 2019), and the previous best fit model from Finke et al. (2010). In all cases, the 68% model errors on the absorption optical depth are quite small, typically ≲10%. The updated model is in quite good agreement with the data and with the previous 2010 model in almost all cases. In the lowest LAT redshift bin, 0.03 < z < 0.23, the new model is ≈20% lower than the old one at most energies. At the highest LAT redshift bin, 2.14 < z < 3.10, agreement between the old and new model is good below 0.1 TeV, although above this energy the old model predicts ≈50% higher τ_{γ γ} . Measurements of τ_{γ γ} are not found in this redshift range for energies > 0.1 TeV, since the values of τ_{γ γ} are likely too high here to be measured (see Section 3.1).

Zoom In Zoom Out Reset image size

IACTs are sensitive to higher-energy γ-rays than the LAT, which makes them sensitive to γ-ray absorption at lower redshifts and higher energies. The absorption measurements with IACTs are shown in the bottom two panels of Figure 6, along with the comparison to our updated model and the Finke et al. (2010) model. Agreement between data and both the models is very good in most places. The data and the older model are ≈30% higher in the 1–6 TeV range than the updated model for both redshift bins. In this energy range, the γ-rays interact with EBL photons in the range ≈4–24 μm (see Equation (1)). In this range, the EBL photons are due to the red end of the stellar photons, and the PAH component, where the old model from Finke et al. (2010) is clearly below the updated model (Figure 5).

The Model A fit is mostly constrained by the luminosity density measurements, since these are more plentiful and have smaller errors. So the discrepancy in the model and observed γ-ray absorption here could be interpreted as a discrepancy between the luminosity density and γ-ray absorption measurements.

The SFRD for the result of Model A is shown in Figure 7. The 68% SFR shows little uncertainty, rarely more than 10%. The curve from Madau & Fragos (2017) is within our model uncertainty at z ≲ 3. But compared with Madau & Dickinson (2014), our model predicts a lower SFRD by a factor ≈2 at z ≲ 3, and by a factor ≈5 at z ≳ 3. The likely discrepancy between our model and that of Madau & Dickinson (2014) is that the latter used a Salpeter IMF, while we used the Baldry–Glazebrook IMF. Madau & Fragos (2017) used the Kroupa IMF (Kroupa 2001) which is more similar to the IMF we used, in that it deviates from a power-law at low masses, with fewer stars. The discrepancy at z ≲ 3 is likely due to different assumptions about the evolution of metallicity and dust extinction with redshift.

Zoom In Zoom Out Reset image size

The resulting stellar mass density and gas-phase metallicity as a function of z are shown in Figure 8. The model mass density reproduces the data quite well. Metallicity measurements from damped Lyα systems (Péroux & Howk 2020 and references therein) corrected for dust depletion following De Cia et al. (2018) are also shown on the figure. Our model is clearly below these data points. We have also performed a model fit identical to Model A, but also including a fit to these points. The result is nearly identical; the metallicity data do not have enough constraining power to impact the model fit. The metallicities from Lyα systems compiled by Péroux & Howk (2020; and also Biffi & Maio 2013) have large dispersions around the mean. This can also be seen in hydrodynamic chemical evolution simulations of early galaxy star formation; these simulations indicate galaxies at z ≈ 9 have 10⁻⁹ < Z < 10⁻⁴ (Biffi & Maio 2013). The model mean metallicity is compared with the mean metallicity curve from Madau & Fragos (2017). For all metallicity conversions in this figure we assumed Z_⊙ = 0.02. The Madau & Fragos (2017) curve was found by convolving the galaxy mass–metallicity relation for star-forming galaxies with various galaxy mass functions at different redshifts. Their results are clearly discrepant with our model, where they get much higher metallicities by ≈1.5–2.5 orders of magnitude. It is not clear if these two methods are really measuring the same thing. The Madau & Fragos (2017) curve is measuring the galaxy mass-weighted mean gas-phase metallicity, while the interpretation of our mean metallicity from a one-zone closed box model is not clear. The mass–metallicity relation used by Madau & Fragos (2017) has been derived from only star-forming galaxies at z ≲ 1.6, so applying it to all galaxies may not be appropriate, nor using it at z ≳ 1.6. The Madau & Fragos (2017) curve is also above almost all of the damped Lyα data points, even when corrected for dust depletion. Madau & Fragos (2017) also note the metallicity normalization of the relation is quite uncertain. It may also be that the closed box model we use is not a good approximation for the metallicity evolution of the universe; indeed it is almost certainly an over-simplification.

Zoom In Zoom Out Reset image size

The 68% contour of our Model A FUV (1500 Å) dust extinction as a function of redshift (Equation (12)) is shown in Figure 9. It clearly agrees with the data. The extinction increases with increasing redshift until about z = 2, after which it decreases. As z → 10, A_FUV(z) → 0, as expected, since dust is made of metals, and at high redshift the metallicity Z → 0.

Zoom In Zoom Out Reset image size

Not only does A(z), i.e., the normalization of the dust extinction curve, evolve, the shape of the curve evolves as well. The dust escape fraction, f_esc,d (λ, z) is shown in Figure 10 as a function of λ for four different redshifts. The general pattern is that as the overall absorption decreases (i.e., more photons escape), it decreases at longer wavelengths much more than shorter wavelengths. At the highest redshift, z = 6, above ≈0.13 μm, the dust escape fraction goes to unity, indicating all the photons escape and none are absorbed. One can also see in this figure that the dust escape fraction decreases (absorption increases) until a peak at about z = 2, and the dust escape fraction increases (absorption decreases) at increasing redshift at z > 2. This is consistent with what is seen in Figure 9.

Zoom In Zoom Out Reset image size

We can compute the contribution to the local (z = 0) EBL intensity from sources at redshifts <z for our model as (Finke et al. 2010)

$$\begin{eqnarray}&&\epsilon I(\epsilon ;\lt z)=\displaystyle \frac{c}{4\pi }{\int }_{0}^{z}{{dz}}_{1}\displaystyle \frac{{\epsilon }^{{\prime} }j({\epsilon }^{{\prime} };{z}_{1})}{1+{z}_{1}}\left|\displaystyle \frac{{{dt}}_{* }}{{{dz}}_{1}}\right|,\end{eqnarray} \tag{ 32 }$$

where ${\epsilon }^{{\prime} }=\epsilon (1+{z}_{1})$. The results of this calculation for Model A for various wavelengths are shown Figure 11. Also shown are the resolved contributions to the build-up of the EBL, as measured by BLAST (Marsden et al. 2009) and Herschel-SPIRE (Béthermin et al. 2012). These are lower limits on the contribution to the local EBL. Based on our modeling, most of the EBL has been resolved by Herschel-SPIRE at 250 μm and 350 μm, while less has been resolved at 500 μm. It is also apparent that at shorter wavelengths for the stellar component of the EBL, the contribution to the local EBL is much more "local" than for the sub-millimeter wavelengths. Most of the local EBL comes from z ≲ 1 for λ < 2.2 μm, while for the dust component in the sub-millimeter, one has to go out to z ≈ 2.5 to have all of the local EBL. This is consistent with the result of Saldana-Lopez et al. (2021).

Zoom In Zoom Out Reset image size

Figure 7 also shows the resulting fit to the LAT γ-ray absorption optical depth data only (Model B). The addition of luminosity and stellar mass density data significantly increase the constraint on the SFRD. Indeed, it is not clear if the γ-ray absorption data are constraining the SFRD much more than the luminosity and stellar mass density data. To test this, we did a run similar to Model A, fitting the luminosity, stellar mass density, and dust extinction data, but without the γ-ray absorption data. This we called model C (Table 3). The results are nearly identical, as seen in Figure 7, with the Model C result having slightly larger uncertainty (≈5%) which is nearly imperceptible in the figure. It does not seem that the γ-ray absorption data are giving a much more tightly constrained SFRD.

In Figure 12, we plot the luminosity density model results for the fit to the γ-ray absorption data only. It demonstrates that the γ-ray data alone can constrain the luminosity density consistent with the observed luminosity density data. Thus, although the γ-ray absorption data do not contribute strongly to the constraint due to their low signal-to-noise ratio, they are valuable as an alternative, independent measurement.

Zoom In Zoom Out Reset image size

The fractions of absorbed light being re-emitted by large warm grains (f₁), small hot grains (f₂), and PAHs are constrained by our models that include luminosity density (Table 3), although uncertainties are large. We find results that are consistent with Finke et al. (2010). For the large warm grains, we find ${f}_{1}={0.56}_{-0.18}^{+0.17}$, compared with f₁ = 0.6 from the previous model; and for the small hot grains we find ${f}_{2}={0.26}_{-0.17}^{+0.18}$, a little higher than the old model (f₂ = 0.05), but within the uncertainty.

The results constrained by only the γ-rays (Model B) demonstrate that the existing γ-ray data do not strongly constrain the dust parameters, although the results of this model are completely consistent with Model A. A similar analysis was performed by Pimentel & Moura-Santos (2019), who used the VHE γ-ray spectrum of Mrk 501 from HEGRA to attempt to constrain these dust parameters as well. Our Model B results are consistent with theirs. Since the HEGRA spectrum for Mrk 501 extended to 20 TeV, they were able to a rule out a single PAH dust component; at least two dust components were required to explain the γ γ absorption seen in this VHE spectrum. As Pimentel & Moura-Santos (2019) discuss, the upcoming CTA will reach yet higher γ-ray energies, and may be able to provide stronger constraints on the dust parameters.

The temperature we find for the warm, large grains, T₁ ≈ 60 K is consistent across all models where this was a free parameter, except Model A.c. It is quiet a bit higher than the value used by Finke et al. (2010), which was T₁ = 40 K. Model A.c gives a value ${T}_{1}={40.5}_{-4.2}^{+11.7}\ {\rm{K}}$, consistent with Finke et al. (2010). This model has H₀ and Ω_m as free parameters, and returns greater uncertainty for this parameter.

We wished to test our model for sensitivity to the IMF and the SFRD parameterization. The results of our model with the Salpeter IMF (Model D) can be seen in Figure 7. The results for Model D agree with Madau & Dickinson (2014) at z ≲ 1, while they are increasingly below that model at higher z. Model D agrees reasonably well with the fiducial Model A (with a Baldry–Glazebrook IMF) at z ≲ 4, but the Salpeter model gives a higher SFRD, more consistent with the one form Madau & Fragos (2017). The difference seems to be that the dust extinction is a free parameter in our model, which for Madau & Dickinson (2014) and Madau & Fragos (2017) it was not in their simpler calculation. We plan to explore different IMFs in more detail in a future publication. The EBL intensity for Model D is still quite low as with Model A, as is the absorption optical depth in the 0.03 < z < 0.23 redshift bin (see Section 4.1).

We also want to test that the results do not depend strongly on the SFRD parameterization. In the right panel of Figure 7 we present our results with our standard MD14 SFRD parameterization and our piecewise SFRD parameterization (Equation (15)). Our results are again very similar for z ≲ 4, and the piecewise SFRD model at z ≳ 4 is higher than ours. The model results do have some sensitivity to the SFRD parameterization as well. Indeed, their seems to be some degeneracy between the SFRD and dust parameters. This is demonstrated in Figure 9, where clearly the dust extinction for Models A and E differ significantly. The greater the dust extinction, the greater the SFRD must be for there to be more light to make up for that extra absorption and fit the luminosity density data. The EBL intensity and absorption optical depths for Model E are very similar to those for the fiducial Model A.

Since the Model A z = 0 EBL intensity is below the lower limits from galaxy counts at many wavelengths (Figure 5), and γ-ray absorption optical depth is consistently below the observed data (Figure 6), we ran model fits with the cosmological parameters H₀ and Ω_m free in the fit, to see if this could resolve this discrepancy. We did model fits with both the BG03 (Model A.c) and Salpeter IMF (Model D.c). These models did not resolve the discrepancy.

Nevertheless, it is interesting to compare the results of these cosmological parameters to values determined with other methods. For Model A.c, our value of ${H}_{0}\,={69.8}_{-3.2}^{+3.6}\ \mathrm{km}\,{{\rm{s}}}^{-1}\ {\mathrm{Mpc}}^{-1}$ is consistent with measurements in the near universe from Type Ia supernovae (73.2 ±1.3 km s⁻¹Mpc⁻¹; Riess et al. 2022) and in the far universe from the CMB (67.4 ± 0.5 km s⁻¹ Mpc⁻¹; Aghanim et al. 2020). The values of ${{\rm{\Omega }}}_{m}={0.23}_{-0.4}^{+0.3}$ is marginally inconsistent (between 2σ and 3σ) with the values measured from Type Ia supernovae (0.327 ± 0.016) and the CMB (0.315 ±0.007). For Model D.c, the value of ${H}_{0}={79.4}_{-4.3}^{+8.1}\ \mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ is 1.4σ of from the Type Ia supernova value, and 2.8σ from the CMB value; the value of Ω_m = 0.29 ± 0.04 is consistent with both the Type Ia supernova and CMB values. It is interesting to note that for our two models with varying cosmological parameters, one is completely consistent with H₀ but slightly discrepant with Ω_m, and vice versa for the other. The main difference between Models A.c and D.c, with different IMFs, is the number of low-mass stars produced. It seems that light produced by low-mass stars in this model has some sensitivity to cosmological parameters.

Domínguez et al. (2019) used the same EBL γ-ray absorption data we use (Section 3.1) to measure cosmological parameters, and their results were ${H}_{0}={67.4}_{-6.2}^{+6.0}\ \mathrm{km}\,{{\rm{s}}}^{-1}\ {\mathrm{Mpc}}^{-1}$ and ${{\rm{\Omega }}}_{m}\,={0.14}_{-0.07}^{+0.06}$. Our errors are much smaller than theirs, likely due to (1) their inclusion of systematic uncertainty, which we do not include, and (2) further constraints imposed on our EBL model from stellar models and other physical constraints. Our values from Model A.c are within 1.3σ of those from Domínguez et al. (2019), while our values for Model D.c are within 2.1σ of theirs. Our values from Model A.c are clearly more consistent with this previous work, including the low value for Ω_m .

We emphasize that the cosmological parameters H₀ and Ω_m from our model fits should not be considered measurements of these cosmological parameters. Unlike the work by Domínguez et al. (2019), we do not include systematic uncertainties. However, it is worth noting the cosmological parameter results depend on the chosen IMF, and light produced by low-mass stars. This can inform future efforts to constrain cosmological parameters with the EBL.