A $\gamma$-ray determination of the Universe's star-formation history

TL;DR

This paper uses gamma-ray observations of distant sources to measure the extragalactic background light, enabling reconstruction of the Universe's star-formation history over most of cosmic time, consistent with galaxy survey data.

Contribution

It introduces a novel gamma-ray based method to trace the Universe's star-formation history, extending measurements to earlier epochs.

Findings

Star-formation rate peaks at redshift z~2.

EBL measurements are consistent with galaxy surveys.

Hints of a decline in faint galaxy abundance at z~6.

Abstract

The light emitted by all galaxies over the history of the Universe produces the extragalactic background light (EBL) at ultraviolet, optical, and infrared wavelengths. The EBL is a source of opacity for $\gamma$ rays via photon-photon interactions, leaving an imprint in the spectra of distant $\gamma$-ray sources. We measure this attenuation using {739} active galaxies and one gamma-ray burst detected by the {\it Fermi} Large Area Telescope. This allows us to reconstruct the evolution of the EBL and determine the star-formation history of the Universe over 90\% of cosmic time. Our star-formation history is consistent with independent measurements from galaxy surveys, peaking at redshift $z\sim2$. Upper limits of the EBL at the epoch of re-ionization suggest a turnover in the abundance of faint galaxies at $z\sim 6$.

Figures

Tables

Table 1. Table S1: Criteria, optimized on simulations, adopted to choose a spectral model.

TSc,1	TSc,2	Model Chosen
		Log-parabola
$>1$		Power law with exponential cut-off with ${\gamma }_1$=0.5
	$>3$	Power law with exponential cut-off with ${\gamma }_1$ free

Table 2. Table S2: Joint-likelihood results for different EBL models (first and second columns) ordered by decreasing value of the last column. The third column refers to the significance, in units of σ 𝜎 \sigma , of the attenuation in the spectra of blazars when a given EBL model is scaled by the factor b 𝑏 b . In this case b 𝑏 b =0 (i.e., no EBL absorption) constitutes the null hypothesis. The fourth column lists the maximum likelihood values and 1  σ 𝜎 \sigma confidence ranges for the opacity scaling factor. In the last column, the b 𝑏 b =1 case (i.e., EBL absorption as predicted by a given EBL model) constitutes the null hypothesis. This column shows the compatibility (expressed in units of σ 𝜎 \sigma ) of the predictions of EBL models with the Fermi -LAT observations. Large values mean less likely to be compatible.

Model	Ref.	Significance of $b$=0	$b$	Significance of $b$=1
		Rejection		Rejection
Scully et al. (2014) – high	(?)	16.0	0.42$\pm 0.03$	17.4
Kneiske et al. (2004) – best -fit	(?)	16.9	0.68$\pm 0.05$	6.0
Gilmore et al. (2012) – fixed	(?)	16.7	1.30$\pm 0.10$	3.0
Gilmore et al. (2012) – fiducial	(?)	16.6	0.81$\pm 0.06$	2.9
Dominguez et al. (2011)	(?)	16.6	1.31$\pm 0.10$	2.9
Franceschini et al. (2017)	(?)	16.4	1.25$\pm 0.10$	2.5
Gilmore et al. (2009)	(?)	16.7	1.03$\pm 0.08$	2.4
Inoue et al. (2013)	(?)	16.2	0.87$\pm 0.06$	2.1
Kneiske & Dole (2010)	(?)	16.8	0.94$\pm 0.08$	1.7
Helgason et al. (2012)	(?)	16.5	1.10$\pm 0.08$	1.3
Finke et al. (2010) – model C	(?)	17.1	1.03$\pm 0.08$	0.4
Scully et al. (2014) – low	(?)	16.0	1.00$\pm 0.07$	0.1

Table 3. Table S3: EBL Optical Depths τ γ ​ γ subscript 𝜏 𝛾 𝛾 \tau_{\gamma\gamma} in bins of redshift and energy as reported in Figure  S3 .

$\overline{z}$	$z_{min}$	$z_{max}$	[10.0-21.4]	[21.4-46.4]	[46.4-100.0]	[100.0-215.4]	[215.4-464.1]	[464.1-1000]
			(GeV)	(GeV)	(GeV)	(GeV)	(GeV)	(GeV)
0.14	0.03	0.23				${0.14}_{-0.04}^{+0.04}$	${0.44}_{-0.09}^{+0.09}$	${0.96}_{-0.21}^{+0.23}$
0.27	0.23	0.34				${0.44}_{-0.14}^{+0.15}$	${0.90}_{-0.28}^{+0.31}$	$>4.32$
0.39	0.34	0.44				${0.34}_{-0.14}^{+0.16}$	${1.59}_{-0.44}^{+0.56}$	$>3.30$
0.52	0.44	0.60				${0.69}_{-0.22}^{+0.24}$	${2.55}_{-0.73}^{+0.88}$	$>3.72$
0.65	0.60	0.68			${0.32}_{-0.14}^{+0.16}$	${1.24}_{-0.42}^{+0.52}$	$>1.61$	$>2.46$
0.79	0.69	0.89			${0.36}_{-0.15}^{+0.16}$	${1.54}_{-0.49}^{+0.61}$	$>5.57$	$>1.73$
0.92	0.89	0.94		${0.05}_{-0.04}^{+0.04}$	${0.51}_{-0.17}^{+0.18}$	${1.92}_{-0.63}^{+0.78}$	$>4.56$	$>1.25$
1.01	0.94	1.09			${0.60}_{-0.16}^{+0.18}$	${1.67}_{-0.50}^{+0.60}$	$>4.38$	$>2.16$
1.24	1.10	1.37	${0.01}_{-0.01}^{+0.01}$	${0.09}_{-0.04}^{+0.04}$	${0.69}_{-0.16}^{+0.17}$	${3.43}_{-0.91}^{+1.32}$	$>5.59$	$>3.16$
1.47	1.38	1.59		${0.21}_{-0.06}^{+0.06}$	${1.48}_{-0.33}^{+0.38}$	$>4.38$	$>2.15$
1.82	1.60	2.14	${0.03}_{-0.02}^{+0.02}$	${0.28}_{-0.08}^{+0.08}$	${1.27}_{-0.30}^{+0.32}$	${4.21}_{-1.59}^{+2.90}$	$>3.03$
2.40	2.15	3.10	${0.05}_{-0.04}^{+0.05}$	${0.60}_{-0.17}^{+0.19}$	${2.97}_{-1.02}^{+1.45}$	$>2.13$

Table 4. Table S4: Spectral intensity a a footnotemark:   ( λ ​ I λ 𝜆 subscript 𝐼 𝜆 \lambda I_{\lambda} ) of the EBL as a function of redshift and wavelengths as reported in Figure  2 . The intensities reported in the this table are in comoving coordinates. To reproduce the results of Figure  2 (shown in physical coordinates) they need to be multiplied by a ( 1 + z ) 3 superscript 1 𝑧 3 (1+z)^{3} factor. A machine-readable version of this table is provided in Data S1 and \href https://figshare.com/s/14f943002230d69a4afdonline.

$\lambda$	$z=0$	$z=1$	$z=2$	$z=3$
($\mu$m)	(nW m-2 sr-1)	(nW m-2 sr-1)	(nW m-2 sr-1)	(nW m-2 sr-1)
0.102	${0.12}_{-0.07}^{+0.12}$	${0.41}_{-0.06}^{+0.08}$	${0.36}_{-0.07}^{+0.11}$	${0.21}_{-0.08}^{+0.07}$
0.111	${0.33}_{-0.20}^{+0.31}$	${0.90}_{-0.13}^{+0.18}$	${0.72}_{-0.14}^{+0.20}$	${0.39}_{-0.15}^{+0.14}$
0.120	${0.59}_{-0.34}^{+0.51}$	${1.42}_{-0.21}^{+0.29}$	${1.05}_{-0.21}^{+0.28}$	${0.55}_{-0.23}^{+0.21}$
0.130	${0.90}_{-0.49}^{+0.68}$	${1.92}_{-0.29}^{+0.38}$	${1.33}_{-0.28}^{+0.36}$	${0.68}_{-0.30}^{+0.28}$
0.141	${1.22}_{-0.63}^{+0.83}$	${2.37}_{-0.37}^{+0.44}$	${1.55}_{-0.34}^{+0.40}$	${0.78}_{-0.35}^{+0.34}$
0.153	${1.55}_{-0.74}^{+0.92}$	${2.74}_{-0.43}^{+0.50}$	${1.71}_{-0.38}^{+0.45}$	${0.84}_{-0.39}^{+0.39}$
0.166	${1.87}_{-0.84}^{+0.98}$	${2.99}_{-0.46}^{+0.56}$	${1.80}_{-0.42}^{+0.47}$	${0.87}_{-0.42}^{+0.43}$
0.180	${2.16}_{-0.87}^{+0.99}$	${3.15}_{-0.48}^{+0.59}$	${1.82}_{-0.44}^{+0.48}$	${0.87}_{-0.42}^{+0.45}$
0.195	${2.44}_{-0.89}^{+0.93}$	${3.21}_{-0.48}^{+0.58}$	${1.79}_{-0.45}^{+0.47}$	${0.84}_{-0.42}^{+0.46}$
0.212	${2.68}_{-0.86}^{+0.87}$	${3.19}_{-0.47}^{+0.55}$	${1.72}_{-0.44}^{+0.46}$	${0.80}_{-0.41}^{+0.45}$
0.230	${2.86}_{-0.79}^{+0.79}$	${3.10}_{-0.44}^{+0.50}$	${1.62}_{-0.42}^{+0.45}$	${0.75}_{-0.38}^{+0.45}$
0.249	${3.01}_{-0.69}^{+0.70}$	${2.98}_{-0.40}^{+0.45}$	${1.52}_{-0.40}^{+0.43}$	${0.71}_{-0.36}^{+0.42}$
0.270	${3.12}_{-0.56}^{+0.60}$	${2.85}_{-0.36}^{+0.41}$	${1.44}_{-0.39}^{+0.40}$	${0.68}_{-0.35}^{+0.39}$
0.293	${3.23}_{-0.48}^{+0.50}$	${2.75}_{-0.33}^{+0.38}$	${1.40}_{-0.39}^{+0.38}$	${0.66}_{-0.36}^{+0.41}$
0.318	${3.33}_{-0.41}^{+0.44}$	${2.72}_{-0.32}^{+0.34}$	${1.39}_{-0.43}^{+0.42}$	${0.68}_{-0.39}^{+0.46}$
0.345	${3.46}_{-0.42}^{+0.41}$	${2.77}_{-0.35}^{+0.38}$	${1.45}_{-0.53}^{+0.52}$	${0.70}_{-0.43}^{+0.56}$
0.374	${3.63}_{-0.48}^{+0.46}$	${2.96}_{-0.43}^{+0.44}$	${1.57}_{-0.65}^{+0.71}$	${0.74}_{-0.48}^{+0.75}$
0.405	${3.87}_{-0.62}^{+0.59}$	${3.27}_{-0.56}^{+0.57}$	${1.77}_{-0.83}^{+0.96}$	${0.84}_{-0.58}^{+0.98}$
0.440	${4.21}_{-0.80}^{+0.75}$	${3.73}_{-0.77}^{+0.72}$	${2.04}_{-1.06}^{+1.28}$	${0.99}_{-0.72}^{+1.26}$
0.477	${4.64}_{-0.98}^{+0.92}$	${4.34}_{-1.01}^{+0.94}$	${2.38}_{-1.33}^{+1.69}$	${1.20}_{-0.92}^{+1.63}$
0.517	${5.19}_{-1.20}^{+1.12}$	${5.06}_{-1.27}^{+1.20}$	${2.80}_{-1.63}^{+2.12}$	${1.47}_{-1.18}^{+2.09}$
0.561	${5.84}_{-1.42}^{+1.32}$	${5.91}_{-1.57}^{+1.48}$	${3.27}_{-1.95}^{+2.58}$	${1.78}_{-1.45}^{+2.66}$
0.608	${6.59}_{-1.65}^{+1.50}$	${6.81}_{-1.93}^{+1.77}$	${3.79}_{-2.28}^{+3.07}$	${2.13}_{-1.75}^{+3.32}$
0.660	${7.41}_{-1.87}^{+1.67}$	${7.71}_{-2.23}^{+2.14}$	${4.32}_{-2.60}^{+3.55}$	${2.47}_{-2.06}^{+4.35}$
0.716	${8.25}_{-2.02}^{+1.82}$	${8.56}_{-2.55}^{+2.50}$	${4.88}_{-2.95}^{+4.01}$	${2.88}_{-2.41}^{+5.38}$
0.776	${9.09}_{-2.13}^{+1.90}$	${9.30}_{-2.84}^{+2.87}$	${5.48}_{-3.35}^{+4.56}$	${3.27}_{-2.72}^{+6.58}$
0.842	${9.88}_{-2.16}^{+1.95}$	${9.90}_{-3.09}^{+3.23}$	${5.98}_{-3.72}^{+5.15}$	${3.64}_{-3.02}^{+7.78}$
0.913	${10.60}_{-2.18}^{+1.96}$	${10.32}_{-3.18}^{+3.56}$	${6.52}_{-4.14}^{+5.88}$	${3.94}_{-3.25}^{+9.86}$
0.990	${11.15}_{-2.06}^{+1.88}$	${10.60}_{-3.36}^{+3.84}$	${6.83}_{-4.40}^{+6.60}$	${4.26}_{-3.52}^{+12.24}$
1.074	${11.54}_{-1.90}^{+1.79}$	${10.73}_{-3.48}^{+4.04}$	${6.97}_{-4.51}^{+7.81}$	${4.51}_{-3.73}^{+15.08}$
1.164	${11.79}_{-1.76}^{+1.65}$	${10.64}_{-3.47}^{+4.23}$	${7.04}_{-4.55}^{+9.00}$	${4.65}_{-3.87}^{+19.68}$
1.263	${11.86}_{-1.58}^{+1.50}$	${10.46}_{-3.65}^{+4.47}$	${6.91}_{-4.47}^{+11.27}$	${4.72}_{-3.94}^{+28.24}$
1.370	${11.73}_{-1.41}^{+1.40}$	${10.05}_{-3.72}^{+4.94}$	${6.87}_{-4.59}^{+13.94}$	${4.79}_{-4.00}^{+40.34}$
1.485	${11.50}_{-1.38}^{+1.28}$	${9.56}_{-3.78}^{+5.80}$	${6.81}_{-4.68}^{+18.49}$	${4.81}_{-4.05}^{+58.87}$
1.611	${11.07}_{-1.35}^{+1.35}$	${9.09}_{-3.86}^{+6.77}$	${6.52}_{-4.61}^{+25.41}$	${4.67}_{-3.96}^{+84.71}$
1.747	${10.53}_{-1.42}^{+1.46}$	${8.54}_{-3.96}^{+8.48}$	${6.20}_{-4.52}^{+37.72}$	${4.51}_{-3.86}^{+117.34}$
1.895	${9.94}_{-1.53}^{+1.63}$	${7.97}_{-4.01}^{+10.94}$	${5.77}_{-4.28}^{+55.07}$	${4.38}_{-3.84}^{+156.29}$
2.055	${9.34}_{-1.70}^{+1.87}$	${7.57}_{-4.07}^{+14.47}$	${5.51}_{-4.22}^{+76.16}$	${4.40}_{-3.90}^{+213.73}$
2.229	${8.72}_{-1.86}^{+2.14}$	${7.01}_{-4.10}^{+20.41}$	${5.38}_{-4.31}^{+101.00}$	${4.47}_{-4.03}^{+271.37}$
2.417	${8.16}_{-2.03}^{+2.55}$	${6.70}_{-4.28}^{+28.53}$	${5.24}_{-4.33}^{+135.49}$	${4.51}_{-4.11}^{+341.38}$
2.621	${7.57}_{-2.16}^{+3.20}$	${6.31}_{-4.30}^{+39.23}$	${5.28}_{-4.50}^{+174.14}$	${4.72}_{-4.37}^{+415.24}$
2.843	${7.05}_{-2.33}^{+4.08}$	${5.95}_{-4.31}^{+52.73}$	${5.38}_{-4.71}^{+223.92}$	${4.86}_{-4.55}^{+494.98}$
3.083	${6.54}_{-2.48}^{+5.51}$
3.344	${6.15}_{-2.69}^{+7.44}$
3.626	${5.85}_{-2.91}^{+9.76}$
3.933	${5.46}_{-3.03}^{+13.33}$
4.265	${5.23}_{-3.21}^{+17.70}$
4.626	${5.11}_{-3.44}^{+23.25}$
5.017	${5.07}_{-3.70}^{+30.88}$

Table 5. Table S5: The Cosmic star-formation history as reported in Figure  3 , also available \href https://figshare.com/s/14f943002230d69a4afdonline.

$z$	Physical EBL model	EBL Reconstruction
	(${10}^{-2}$M⊙ yr-1 Mpc-3)	(${10}^{-2}$M⊙ yr-1 Mpc-3)
0.0	${0.8}_{-0.3}^{+0.7}$	${1.2}_{-0.7}^{+1.3}$
0.1	${1.1}_{-0.4}^{+0.8}$	${1.6}_{-0.9}^{+1.4}$
0.2	${1.6}_{-0.5}^{+0.9}$	${2.1}_{-1.1}^{+1.4}$
0.3	${2.0}_{-0.5}^{+0.9}$	${2.7}_{-1.2}^{+1.4}$
0.4	${2.6}_{-0.6}^{+0.8}$	${3.4}_{-1.4}^{+1.4}$
0.5	${3.3}_{-0.6}^{+0.8}$	${4.2}_{-1.5}^{+1.4}$
0.6	${4.1}_{-0.7}^{+0.8}$	${5.0}_{-1.5}^{+1.4}$
0.8	${6.0}_{-0.9}^{+1.0}$	${6.8}_{-1.4}^{+1.5}$
1.0	${8.2}_{-1.4}^{+1.6}$	${8.5}_{-1.3}^{+1.6}$
1.2	${10.7}_{-2.2}^{+2.0}$	${10.0}_{-1.5}^{+2.1}$
1.4	${12.9}_{-2.8}^{+2.5}$	${11.0}_{-1.8}^{+2.3}$
1.6	${14.5}_{-2.9}^{+2.3}$	${11.3}_{-2.0}^{+2.9}$
1.8	${15.0}_{-2.5}^{+2.3}$	${11.0}_{-2.0}^{+3.1}$
2.0	${14.7}_{-2.4}^{+2.4}$	${10.4}_{-1.9}^{+3.0}$
2.5	${11.6}_{-2.9}^{+3.2}$	${7.9}_{-1.7}^{+2.2}$
3.0	${7.6}_{-2.4}^{+3.3}$	${5.6}_{-1.8}^{+1.8}$
3.5	${4.8}_{-1.8}^{+2.6}$	${4.0}_{-1.7}^{+1.5}$
4.0	${3.1}_{-1.3}^{+2.0}$	${2.8}_{-1.5}^{+1.5}$
4.5	${2.0}_{-0.9}^{+1.6}$	${2.1}_{-1.3}^{+1.4}$
5.0	${1.3}_{-0.7}^{+1.3}$	${1.6}_{-1.1}^{+1.4}$
5.5	${0.9}_{-0.5}^{+1.1}$	${1.2}_{-0.9}^{+1.3}$
6.0	${0.6}_{-0.3}^{+1.0}$	${1.0}_{-0.8}^{+1.3}$

Equations

$${\left(\frac{\mathrm{d}N}{\mathrm{d}E}\right)}_{\mathrm{obs}}={\left(\frac{\mathrm{d}N}{\mathrm{d}E}\right)}_{\mathrm{int}}\times {\mathrm{e}}^{-b\cdot {\tau }_{\gamma \gamma }(E,z)}$$

$$\frac{\mathrm{d}N}{\mathrm{d}E}=N_0{\left(\frac{E}{E_b}\right)}^{-\alpha +\beta \log (E/E_b)}$$

$$\frac{\mathrm{d}N}{\mathrm{d}E}=N_0{\left(\frac{E}{E_c}\right)}^{\alpha }{e}^{-(E/E_b{)}^{{\gamma }_1}}$$

$$T{S}_{\mathrm{c},1}$$

$$T{S}_{\mathrm{c},2}$$

$${\tau }_{\gamma \gamma }(E_{\gamma },z_s)=c{\int }_0^{z_s}\left|\frac{dt}{dz}\right|dz{\int }_{-1}^1(1-\mu )\frac{d\mu }{2}{\int }_{2m_e^2{c}^4/{ϵ}_{\gamma }(1-\mu )}^{\infty }\sigma ({ϵ}_{\mathrm{EBL}},{ϵ}_{\gamma },\mu )n_{\mathrm{EBL}}(ϵ,z)d{ϵ}_{\mathrm{EBL}}$$

$$\sigma ({ϵ}_{\mathrm{EBL}},{ϵ}_{\gamma },\mu )=\frac{3{\sigma }_T}{16}(1-{\beta }^2)\left[2\beta ({\beta }^2-2)+(3-{\beta }^4)\ln \left(\frac{1+\beta }{1-\beta }\right)\right],$$

$$\beta =\sqrt{1-\frac{2m_e^2{c}^4}{{ϵ}_{\mathrm{EBL}}{ϵ}_{\gamma }(1-\mu )}}.$$

$$n_{\mathrm{EBL}}(ϵ,z)=(1+z{)}^3{\int }_z^{\infty }\frac{j({ϵ}^{\prime },\bar{z})}{{ϵ}^{\prime }}\left|\frac{dt}{d\bar{z}}\right|d\bar{z},$$

$$j(\lambda )=\sum _i{a}_i\cdot \exp \left[-\frac{{\left(\log \lambda -\log {\lambda }_i\right)}^2}{2{\sigma }^2}\right][\mathrm{erg}\cdot {\mathrm{s}}^{-1}\mathrm{c}{\mathrm{m}}^{-3}\mathrm{Hz}{}^{-1}]$$

$$j({\lambda }_i,z)=\sum _i{a}_i\cdot \exp \left[-\frac{{\left(\log \lambda -\log {\lambda }_i\right)}^2}{2{\sigma }^2}\right]\times \frac{(1+z{)}^{b_i}}{1+{\left(\frac{1+z}{c_i}\right)}^{d_i}}.$$

$$j({\lambda }_i,z)\propto \frac{a_i+b_{i}z}{1+(z/c_i{)}^{d_i}},$$

$${\chi }^2=\sum _{i=1}^N\sum _{j=1}^M\frac{[{\tau }_{\mathrm{data}}(E_i,z_j)-{\tau }_{\mathrm{model}}(\vec{\theta }|E_i,z_j){]}^2}{{\sigma }_{i,j}^2}$$

$$j_{\mathrm{high}-\mathrm{z}}(\lambda ,z)=a_{\mathrm{high}-\mathrm{z}}\exp \left[-\frac{(z-z_0{)}^2}{2{\sigma }^2}\right]{\left(\frac{\lambda }{0.17\text{upmu}\mathrm{m}}\right)}^{-0.5}$$

$$j_{\mathrm{UV}}=2.5\times 10^{26}{ϵ}_{53}^{-1}{\left(\frac{1+z}{6}\right)}^3{\left(\frac{{\Omega }_b{h}_{70}^2}{0.0461}\right)}^2{\left(\frac{C/f_{\mathrm{esc}}}{30}\right)}^2\mathrm{erg}{\mathrm{s}}^{-1}\mathrm{Mp}{\mathrm{c}}^{-3}\mathrm{H}{\mathrm{z}}^{-1}.$$

$$A_V(z)=(1.49\pm 0.07)\frac{(1+z{)}^{(0.64\pm 0.19)}}{1+[(1+z)/(3.40\pm 0.44){]}^{(3.54\pm 0.47)}}.$$

$${\rho }_{\star }(z)=j_{\mathrm{UV}}(z)\cdot {\mathcal{K}}_{\mathrm{UV}}\cdot 10^{0.4A_V}$$

Full text

A $\gamma$-ray determination of the Universe’s star-formation history

The Fermi-LAT Collaboration*∗,†*

$*$ The list of authors and affiliations can be found at the beginning of the supplementary material

$\dagger$ Corresponding Authors: majello@g.clemson.edu, helgason@hi.is, vpaliya@g.clemson.edu,

justin.finke@nrl.navy.mil, abhishd@g.clemson.edu, alberto@gae.ucm.es

The light emitted by all galaxies over the history of the Universe produces the extragalactic background light (EBL) at ultraviolet, optical, and infrared wavelengths. The EBL is a source of opacity for $\gamma$ rays via photon-photon interactions, leaving an imprint in the spectra of distant $\gamma$-ray sources. We measure this attenuation using 739 active galaxies and one gamma-ray burst detected by the Fermi Large Area Telescope. This allows us to reconstruct the evolution of the EBL and determine the star-formation history of the Universe over 90% of cosmic time. Our star-formation history is consistent with independent measurements from galaxy surveys, peaking at redshift $z\sim 2$. Upper limits of the EBL at the epoch of re-ionization suggest a turnover in the abundance of faint galaxies at $z\sim 6$.

Stars produce the bulk of the optical light in the Universe and synthesize most of the elements found in galaxies. The cosmic star-formation history (SFH), i.e. the stellar birth rate as a function of the Universe’s age, summarizes the history of stellar formation since the Big Bang  (?). The rate of star formation is commonly estimated by measuring direct emission of light from massive short-lived stars, typically in the ultraviolet (UV) and (or) by detecting the reprocessed radiation from dusty star-forming regions in the infrared (IR). The conversion from the UV light emitted by a minority of stars to the stellar mass formed per year relies on assumptions about the mass distribution of the newly formed stellar population (the initial mass function, IMF), the element enrichment history of the interstellar medium, and obscuration by dust. Such estimates of the SFH rely on the detection of many individual galaxies in deep surveys  (?, ?, ?). Because not even the most powerful telescope can detect all the galaxies in a representative field, one of the largest sources of uncertainty in the SFH is estimating the amount of light from undetected galaxies, and the star formation associated with them. This difficulty becomes particularly relevant within the first billion years after the Big Bang when a large population of faint, still undetected, galaxies existed  (?). These galaxies are expected to drive the re-ionization of the Universe: the period when energetic UV photons from young stars escaped into intergalactic space and ionized the neutral hydrogen of the intergalactic medium. Similarly, recent (i.e. within one billion years from the present age) star formation measured using space-borne UV observatories is based on surveys extending over small solid angles  (?), and are therefore subject to density fluctuations in the large-scale structure, an effect known as cosmic variance.

Observational estimates of the SFH are sufficiently uncertain that measurements with multiple independent methodologies are desirable. Starlight that escapes galaxies is almost never destroyed and becomes part of the extragalactic background light (EBL), the total light accumulated by all sources over the lifetime of the Universe  (?, ?, ?). While extremely important, accurate measurements of this diffuse all-sky background at UV to IR wavelengths, and particularly its build-up over time, have only just become possible  (?).

We present an alternative approach to measure the SFH based on the attenuation that the EBL produces in the $\gamma$-ray spectra of distant sources. $\gamma$ rays with sufficient energy can annihilate when they collide with EBL photons and produce electron-positron pairs (i.e. the reaction $\gamma\gamma\rightarrow e^{+}e^{-}$), effectively being absorbed as a result of the interaction  (?). Above a given threshold energy, the attenuation experienced by every $\gamma$-ray source at a similar distance depends on the number density of the EBL target photons integrated along the line of sight; observations of $\gamma$-ray sources at different distances (as measured by the sources’ redshifts) can be used to measure the density of EBL photons at different cosmic times.

We analyze $\gamma$-ray photons detected by the Large Area Telescope (LAT) instrument on the Fermi Gamma-ray Space Telescope, over 9 years of operations. Our sample of suitable objects for this analysis consists of 739 blazars, galaxies hosting a super-massive black hole with a relativistic jet pointed at a small angle to the line of sight. The distances of these blazars correspond to lookback times of 0.2-11.6 billion years according to the standard cosmological model  (?). We perform a likelihood analysis to find the EBL attenuation experienced by all blazars whilst simultaneously optimizing the spectral parameters independently for each blazar  (?). This is accomplished individually for each source, by defining a region of interest that comprises all $\gamma$ rays detected within 15*∘* of the source position and creating a sky model that includes all sources of $\gamma$ rays in the field. The parameters of the sky model are then optimized by a maximum likelihood method. For every blazar, the fitting is performed below an energy at which the EBL attenuation is negligible and thus yields a measurement of the intrinsic (i.e., unabsorbed) blazar spectrum. The intrinsic spectra are described using simple empirical functions  (?) and extrapolated to higher energy, where the $\gamma$ rays are expected to be attenuated by the EBL.

Potential EBL absorption is added to the fitted spectra as follows:

[TABLE]

where $\left(\frac{{\rm d}N}{{\rm d}E}\right)_{\rm obs}$ and $\left(\frac{{\rm d}N}{{\rm d}E}\right)_{\rm int}$ are the observed and intrinsic blazar spectra respectively, $\tau_{\gamma\gamma}(E,z)$ is the EBL optical depth as estimated from models (at a given energy $E$ and redshift $z$) and $b$ is a free parameter. The data from all blazars are combined to yield the best-fitting value of $b$ for each model. A value of $b=0$ implies no EBL attenuation is present in the spectra of blazars, while $b\approx 1$ implies an attenuation compatible with the model prediction. Twelve of the most recent models that predict the EBL attenuation up to a redshift of $z=3.1$ have been tested in this work. We detect the attenuation due to the EBL in the spectra of blazars at $\gtrsim 16$ standard deviations ($\sigma$) for all models tested (see Table S2).

Our analysis leads to detections of the EBL attenuation across the entire $0.03<z<3.1$ redshift range of the blazars. From this, we identify the redshift at which, for a given energy, the Universe becomes opaque to $\gamma$ rays, known as cosmic $\gamma$-ray horizon (Figure 1). With the optical depths measured in six energy bins ($10-1000$ GeV) across twelve redshift bins  (?) we are able to reconstruct the intensity of the EBL at different epochs (Figure 2). We model the cosmic emissivity (luminosity density) of sources as several simple spectral components at UV, optical, and near-IR (NIR) wavelengths. These components are allowed to vary in amplitude and evolve with redshift independently of each other to reproduce, through a Markov Chain Monte Carlo (MCMC) analysis, the optical depth data. The emissivities as a function of wavelength and redshift allow us to reconstruct the history of the EBL over $\sim 90$ % of cosmic time.

At $z=0$ the energy spectrum of the EBL is close to the one inferred by resolving individual galaxies in deep fields  (?). At all other epochs, Fermi LAT is most sensitive to the UV-optical component of the EBL, and is only able to constrain the NIR component at more recent times (see Figure 2). The intensity of the UV background in the local Universe remains uncertain, with independent studies reporting differing values  (?, ?, ?). Our determination of 2.56${}^{+0.92(2.23)}_{-0.87(1.49)}$ nW m*-2* sr*-1*, 1$\sigma$(2$\sigma$), at 0.2 $\upmu$m favors an intermediate UV intensity in agreement with  (?). In the NIR our measurement of $11.6^{+1.3(2.6)}_{-1.4(3.1)}~{}{\rm nW~{}m^{-2}sr^{-1}}$, 1$\sigma$ (2$\sigma$), at 1.4 $\upmu$m is consistent with integrated galaxy counts  (?, ?), leaving little room for additional components, contrary to some suggestions  (?, ?). This notably includes contributions from stars that have been stripped from galaxies as the technique presented here is sensitive to all photons  (?, ?).

At any epoch, the EBL is composed of the emission of all stars  (?) that existed up to that point in time and can therefore be used to infer properties related to the evolution of galaxy populations. We focus on the cosmic SFH, which we determine using two independent methods. First, we use the reconstructed UV emissivity across cosmic time to derive the SFH from established relations between the UV luminosity and star-formation rate  (?), taking into account the mean dust extinction within galaxies  (?, ?, ?). The second approach uses a physical EBL model  (?) to calculate the optical depth due to the EBL directly from the SFH. The SFH is then optimized using a MCMC to reproduce the Fermi-LAT optical depth data  (?). The two approaches yield consistent results for the SFH, which is well constrained out to a redshift of $z\approx 5$, i.e., to the epoch 1.5 billion years after the Big Bang (Figure 3).

Because the optical depth increases with the distance traveled by the $\gamma$ rays, we obtain the tightest constraints in the redshift range $0.1<z<1.5$, beyond which our sensitivity decreases due to the lower number of observed blazars. To improve the constraint of the SFH beyond $z=3$, we have complemented the blazar sample with a gamma-ray burst (GRB 080916C) at $z=4.35$  (?). This allows us to place upper limits on the SFH at $z\gtrsim 5$, because photons generated at redshifts higher than the $z=4.35$ limit of our sample remain in the EBL, become redshifted, and start interacting with the $\gamma$ rays from the blazars and the GRB used here at $z<4.35$.

At $z\gtrsim 6$ the far-UV background (photon energy $>13.6$ eV) is responsible for the re-ionization of the neutral hydrogen in the Universe, but the nature of ionizing sources has not been conclusively identified. One possibility is that ultra-faint galaxies existing in large numbers can provide the required ionizing photons  (?, ?). In this case, the galaxy UV luminosity function must be steep at the faint end. Recent measurements of the luminosity function in the deepest Hubble fields remain inconclusive at the faintest levels (absolute AB magnitude $M_{\rm AB}\gtrsim-15$) with some suggesting a continued steep faint-end slope  (?, ?) and others claiming a turnover  (?, ?). Our upper limits at $z=5-6$ on the UV emissivity $\rho_{\rm UV}<3.2(5.3)\times 10^{26}$ erg s*-1* Mpc*-3* Hz*-1* 1 $\sigma$ (2 $\sigma$), see Figure 4, suggest a turnover of the luminosity function at $M_{\rm AB}\sim-14$ in agreement with  (?) and  (?). This still allows for abundant photons to drive the re-ionization.

References and Notes

P. Madau, M. Dickinson, ARA&A** 52, 415 (2014).

• 2.

N. A. Grogin, et al., ApJS** 197, 35 (2011).

• 3.

G. D. Illingworth, et al., ApJS** 209, 6 (2013).

• 4.

J. M. Lotz, et al., ApJ** 837, 97 (2017).

• 5.

R. J. McLure, et al., MNRAS** 432, 2696 (2013).

• 6.

D. Schiminovich, et al., ApJ** 619, L47 (2005).

• 7.

M. G. Hauser, E. Dwek, ARA&A** 39, 249 (2001).

• 8.

A. Kashlinsky, Phys. Rep.** 409, 361 (2005).

• 9.

E. Dwek, F. Krennrich, Astroparticle Physics 43, 112 (2013).

• 10.

S. K. Andrews, et al., MNRAS** 470, 1342 (2017).

• 11.

A. Nikishov, Journal of Experimental and Theoretical Physics 14, 393 (1962).

• 12.

We adopted the following values for the Hubble and cosmological parameters: $H_{0}$=70 km s*-1* Mpc*-1*, $\Omega_{M}$=0.3, and $\Omega_{\Lambda}$=0.7.

• 13.

M. Ackermann, et al., Science 338, 1190 (2012).

• 14.

Materials and methods are available as supplementary materials on Science Online.

• 15.

S. P. Driver, et al., ApJ** 827, 108 (2016).

• 16.

J. P. Gardner, T. M. Brown, H. C. Ferguson, ApJ** 542, L79 (2000).

• 17.

C. K. Xu, et al., ApJ** 619, L11 (2005).

• 18.

E. N. Voyer, J. P. Gardner, H. I. Teplitz, B. D. Siana, D. F. de Mello, ApJ** 736, 80 (2011).

• 19.

R. C. Keenan, A. J. Barger, L. L. Cowie, W.-H. Wang, ApJ** 723, 40 (2010).

• 20.

M. L. N. Ashby, et al., ApJS** 209, 22 (2013).

• 21.

R. A. Bernstein, ApJ** 666, 663 (2007).

• 22.

S. Matsuura, et al., ApJ** 839, 7 (2017).

• 23.

M. Zemcov, et al., Science 346, 732 (2014).

• 24.

C. Burke, M. Hilton, C. Collins, MNRAS** 449, 2353 (2015).

• 25.

The contribution of active galactic nuclei is small in comparison, see Supplementary Materials.

• 26.

R. C. Kennicutt, Jr., ARA&A** 36, 189 (1998).

• 27.

R. J. Bouwens, et al., ApJ** 754, 83 (2012).

• 28.

D. Burgarella, et al., A&A** 554, A70 (2013).

• 29.

J. D. Finke, S. Razzaque, C. D. Dermer, ApJ** 712, 238 (2010).

• 30.

A. A. Abdo, et al., Science 323, 1688 (2009).

• 31.

S. L. Finkelstein, et al., ApJ** 810, 71 (2015).

• 32.

B. E. Robertson, R. S. Ellis, S. R. Furlanetto, J. S. Dunlop, ApJ** 802, L19 (2015).

• 33.

R. C. Livermore, S. L. Finkelstein, J. M. Lotz, ApJ** 835, 113 (2017).

• 34.

M. Ishigaki, et al., ApJ** 854, 73 (2018).

• 35.

R. J. Bouwens, P. A. Oesch, G. D. Illingworth, R. S. Ellis, M. Stefanon, ApJ** 843, 129 (2017).

• 36.

H. Atek, J. Richard, J.-P. Kneib, D. Schaerer, arXiv:1803.09747 (2018).

• 37.

A. Domínguez, et al., MNRAS** 410, 2556 (2011).

• 38.

K. Helgason, A. Kashlinsky, ApJ** 758, L13 (2012).

• 39.

A. Abramowski, et al., A&A** 550, A4 (2013).

• 40.

J. Biteau, D. A. Williams, ApJ** 812, 60 (2015).

• 41.

H. Abdalla, et al., A&A** 606, A59 (2017).

• 42.

M. Meyer, M. Raue, D. Mazin, D. Horns, A&A** 542, A59 (2012).

• 43.

M. Ackermann, et al., ApJ** 810, 14 (2015).

• 44.

M. S. Shaw, et al., ApJ** 764, 135 (2013).

• 45.

F. Acero, et al., ApJS** 223, 26 (2016).

• 46.

The templates used are gll$\_$iem$\_$v06.fits and iso$\_$P8R2$\_$SOURCE$\_$V6$\_$v06.txt. These are available at \urlhttps://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.

• 47.

F. Acero, et al., ApJS** 218, 23 (2015).

• 48.

M. Ajello, et al., ApJ** 751, 108 (2012).

• 49.

S. T. Scully, M. A. Malkan, F. W. Stecker, ApJ** 784, 138 (2014).

• 50.

T. M. Kneiske, T. Bretz, K. Mannheim, D. H. Hartmann, A&A** 413, 807 (2004).

• 51.

R. C. Gilmore, R. S. Somerville, J. R. Primack, A. Domínguez, MNRAS** 422, 3189 (2012).

• 52.

A. Franceschini, G. Rodighiero, A&A** 603, A34 (2017).

• 53.

R. C. Gilmore, P. Madau, J. R. Primack, R. S. Somerville, F. Haardt, MNRAS** 399, 1694 (2009).

• 54.

Y. Inoue, et al., ApJ** 768, 197 (2013).

• 55.

T. M. Kneiske, H. Dole, A&A** 515, A19 (2010).

• 56.

G. G. Fazio, F. W. Stecker, Nature** 226, 135 (1970).

• 57.

A. Domínguez, et al., ApJ** 770, 77 (2013).

• 58.

J. Greiner, et al., A&A** 498, 89 (2009).

• 59.

R. C. Gilmore, MNRAS** 420, 800 (2012).

• 60.

Y. Inoue, Y. T. Tanaka, G. M. Madejski, A. Domínguez, ApJ** 781, L35 (2014).

• 61.

W. B. Atwood, et al., ApJ** 774, 76 (2013).

• 62.

A. Desai, et al., ApJ** 850, 73 (2017).

• 63.

We used the following templates: gll_iem_v06.fits and iso_P8R2_TRANSIENT020_V6_v06.txt. These are available at \urlhttps://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.

• 64.

A. Tramacere, P. Giommi, M. Perri, F. Verrecchia, G. Tosti, A&A** 501, 879 (2009).

• 65.

A. Tramacere, E. Massaro, A. M. Taylor, ApJ** 739, 66 (2011).

• 66.

M. Ajello, et al., ApJS** 232, 18 (2017).

• 67.

M. Ackermann, et al., ApJ** 771, 57 (2013).

• 68.

S. Abdollahi, et al., ApJ** 846, 34 (2017).

• 69.

M. Ackermann, et al., ApJS** 203, 4 (2012).

• 70.

R. J. Gould, G. P. Schréder, Physical Review 155, 1404 (1967).

• 71.

F. W. Stecker, Cosmic gamma rays, vol. 249 (NASA, 1971).

• 72.

R. W. Brown, K. O. Mikaelian, R. J. Gould, Astrophys. Lett.** 14, 203 (1973).

• 73.

S. Cole, et al., MNRAS** 326, 255 (2001).

• 74.

D. Foreman-Mackey, D. W. Hogg, D. Lang, J. Goodman, PASP** 125, 306 (2013).

• 75.

J. Goodman, J. Weare, Communications in Applied Mathematics and Computational Science 5, 65 (2010).

• 76.

F. W. Stecker, S. T. Scully, M. A. Malkan, ApJ** 827, 6 (2016).

• 77.

S. K. Andrews, S. P. Driver, L. J. M. Davies, C. d. P. Lagos, A. S. G. Robotham, MNRAS** 474, 898 (2018).

• 78.

V. Khaire, R. Srianand, ApJ** 805, 33 (2015).

• 79.

E. L. Wright, E. D. Reese, ApJ** 545, 43 (2000).

• 80.

T. Matsumoto, et al., ApJ** 626, 31 (2005).

• 81.

F. Aharonian, et al., Nature** 440, 1018 (2006).

• 82.

Albert, et al., Science 320, 1752 (2008).

• 83.

A. Domínguez, M. Ajello, ApJ** 813, L34 (2015).

• 84.

T. Armstrong, A. M. Brown, P. M. Chadwick, MNRAS** 470, 4089 (2017).

• 85.

M. L. Ahnen, et al., ApJ** 815, 23L (2015).

• 86.

M. L. Ahnen, et al., A&A 595, 98 (2016).

• 87.

P. Madau, After the Dark Ages: When Galaxies were Young (the Universe at 2 $<$ Z $<$ 5), S. Holt, E. Smith, eds. (1999), vol. 470 of American Institute of Physics Conference Series, pp. 299–311.

• 88.

S. L. Finkelstein, et al., ApJ** 758, 93 (2012).

• 89.

A. Kashlinsky, ApJ** 633, L5 (2005).

• 90.

R. J. Bouwens, et al., ApJ** 793, 115 (2014).

• 91.

R. S. Ellis, et al., ApJ** 763, L7 (2013).

• 92.

G. Chabrier, Publ. Astron. Soc. Pac., 115, 763, (2003).

• 93.

E. .E. Salpeter, ApJ**, 121, 161 (1995).

• 94.

R. J. Bouwens, et al., ApJ** 833, 72 (2016).

• 95.

O. Cucciati, et al., A&A** 539, A31 (2012).

• 96.

T. T. Takeuchi, V. Buat, D. Burgarella, A&A** 440, L17 (2005).

• 97.

N. A. Reddy, C. C. Steidel, ApJ** 692, 778 (2009).

• 98.

T. K. Wyder, et al., ApJ** 619, L15 (2005).

• 99.

A. S. G. Robotham, S. P. Driver, MNRAS** 413, 2570 (2011).

• 100.

T. Dahlen, et al., ApJ** 654, 172 (2007).

• 101.

M. Yoshida, et al., ApJ** 653, 988 (2006).

• 102.

L. Tresse, et al., A&A** 472, 403 (2007).

• 103.

M. Raue, M. Meyer, apj 426, 1097 (2012).

• 104.

N. Bourne, et al., MNRAS** 467, 1360 (2017).

• 105.

M. D. Kistler, H. Yüksel, J. F. Beacom, A. M. Hopkins, J. S. B. Wyithe, ApJ** 705, L104 (2009).

• 106.

S. P. Driver, et al., MNRAS** 475, 2891 (2018).

• 107.

R. J. Bouwens, et al., ApJ** 803, 34 (2015).

• 108.

S. Razzaque, C. D. Dermer, J. D. Finke, ApJ** 697, 483 (2009).

• 109.

P. P. Eggleton, M. J. Fitchett, C. A. Tout, ApJ** 347, 998 (1989).

• 110.

S. P. Driver, et al., ApJ** 678, L101 (2008).

• 111.

Y. Gong, A. Cooray, ApJ** 772, L12 (2013).

• 112.

A. M. Hopkins, J. F. Beacom, ApJ** 651, 142 (2006).

• 113.

P. Madau, F. Haardt, ApJ** 813, L8 (2015).

• 114.

E. Lusso, et al., MNRAS** 449, 4204 (2015).

• 115.

V. Khaire, R. Srianand, arXiv:1801.09693 (2018).

Acknowledgments

Funding

M. Ajello, V. Paliya, and A. Desai acknowledge support from NSF and NASA through grants AST-1715256, NNX16AR72G, and 80NSSC17K0506. K. Helgason acknowledges support from the Icelandic Research Fund, grant number 173728-051. J. Finke was supported by the Chief of Naval Research and by a grant of computer time from the Department of Defense High Performance Computing Modernization Program at the Naval Research Laboratory. A. Domínguez thanks the support of the Juan de la Cierva program from the Spanish MEC.

Fermi-LAT collaboration: The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515. J. Conrad acknowledge funding from Wallenberg Academy. M. Razzano is funded by contract FIRB-2012-RBFR12PM1F from the Italian Ministry of Education, University and Research (MIUR)

Author contribution:

All authors meet the journal authorship criteria. M. Ajello designed the project and wrote most of the paper. V. Paliya performed the analysis of the $\gamma$-ray data, while K. Helgason derived the constrains on the EBL and the star-formation history and wrote all the corresponding text. J. Finke derived the results of the stellar population method and wrote the corresponding text. A. Desai tested all the EBL models while A. Domínguez provided all the EBL-related data reported in the figures and wrote the corresponding text. All co-authors have read, provided comments and approved the manuscript.

Competing interests:

All co-authors declare that there are no competing interesting.

Data and materials availability:

The data used to derive the results presented in this paper are provided in tabular form in the supplementary materials. The Fermi-LAT data and software needed to analyze those are available from the Fermi Science Support Center \urlhttp://fermi.gsfc.nasa.gov/ssc. The reconstructed optical depth templates, EBL, and SFH are also available at \urlhttps://figshare.com/s/14f943002230d69a4afd. The tool to produce physical models of blazars’ SEDs is available at \urlwww.isdc.unige.ch/sedtool.

Supplementary Materials

\url

www.sciencemag.org

Authors and Affiliations

Materials and Methods

Figs. S1 to S12

Table S1 to S5

References (43-115)

Authors and Affiliations

S. Abdollahi1, M. Ackermann2, M. Ajello${}^{3}\dagger$, W. B. Atwood4, L. Baldini5, J. Ballet6, G. Barbiellini7,8, D. Bastieri9,10, J. Becerra Gonzalez11,12, R. Bellazzini13, E. Bissaldi14,15, R. D. Blandford16, E. D. Bloom16, R. Bonino17,18, E. Bottacini16,19, S. Buson11, J. Bregeon20, P. Bruel21, R. Buehler2, R. A. Cameron16, R. Caputo22, P. A. Caraveo23, E. Cavazzuti24, E. Charles16, S. Chen9,19, C. C. Cheung25, G. Chiaro23, S. Ciprini26,27, J. Cohen-Tanugi20, L. R. Cominsky28, J. Conrad29,30, D. Costantin10, S. Cutini26,27, F. D’Ammando31,32, F. de Palma15,33, A. Desai${}^{3}\dagger$, S. W. Digel16, N. Di Lalla5, M. Di Mauro16, L. Di Venere14,15, A. Domínguez${}^{34}\dagger$, C. Favuzzi14,15, S. J. Fegan21, J. Finke${}^{25}\dagger$, A. Franckowiak2, Y. Fukazawa1, S. Funk35, P. Fusco14,15, G. Gallardo Romero2,36, F. Gargano15, D. Gasparrini26,27, N. Giglietto14,15, F. Giordano14,15, M. Giroletti31, D. Green12,11, I. A. Grenier6, L. Guillemot37,38, S. Guiriec39,11, D. H. Hartmann3, E. Hays11, K. Helgason${}^{40,42}\dagger$, D. Horan21, G. Jóhannesson41,42, D. Kocevski11, M. Kuss13, S. Larsson43,30, L. Latronico17, J. Li2, F. Longo7,8, F. Loparco14,15, B. Lott44, M. N. Lovellette25, P. Lubrano27, G. M. Madejski16, J. D. Magill12, S. Maldera17, A. Manfreda5, L. Marcotulli3, M. N. Mazziotta15, J. E. McEnery11,12, M. Meyer16, P. F. Michelson16, T. Mizuno45, M. E. Monzani16, A. Morselli46, I. V. Moskalenko16, M. Negro17,18, E. Nuss20, R. Ojha11, N. Omodei16, M. Orienti31, E. Orlando16, J. F. Ormes47, M. Palatiello7,8, V. S. Paliya${}^{3}\dagger$, D. Paneque48, J. S. Perkins11, M. Persic7,49, M. Pesce-Rollins13, V. Petrosian16, F. Piron20, T. A. Porter16, J. R. Primack4, G. Principe35, S. Rainò14,15, R. Rando9,10, M. Razzano13,51, S. Razzaque50, A. Reimer51,16, O. Reimer51,16, P. M. Saz Parkinson4,52,53, C. Sgrò13, E. J. Siskind54, G. Spandre13, P. Spinelli14,15, D. J. Suson55, H. Tajima56,16, M. Takahashi48, J. B. Thayer16, L. Tibaldo57,58, D. F. Torres59,60, E. Torresi61, G. Tosti27,62, A. Tramacere63, E. Troja11,12, J. Valverde21, G. Vianello16, M. Vogel64, K. Wood65, G. Zaharijas66,67

Department of Physical Sciences, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan 2. 2.

Deutsches Elektronen Synchrotron DESY, D-15738 Zeuthen, Germany 3. 3.

Department of Physics and Astronomy, Clemson University, Kinard Lab of Physics, Clemson, SC 29634-0978, USA 4. 4.

Santa Cruz Institute for Particle Physics, Department of Physics and Department of Astronomy and Astrophysics, University of California at Santa Cruz, Santa Cruz, CA 95064, USA 5. 5.

Università di Pisa and Istituto Nazionale di Fisica Nucleare, Sezione di Pisa I-56127 Pisa, Italy 6. 6.

Laboratoire Astrophysique Interactions Multi-échelles, Commissariat á l’énergie atomique-Institute of Research into the Fundamental Laws of the Universe/Centre national de la recherche scientifique/Université Paris Diderot, Service d’Astrophysique, Commissariat á l’énergie atomique Saclay, F-91191 Gif sur Yvette, France 7. 7.

Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, I-34127 Trieste, Italy 8. 8.

Dipartimento di Fisica, Università di Trieste, I-34127 Trieste, Italy 9. 9.

Istituto Nazionale di Fisica Nucleare, Sezione di Padova, I-35131 Padova, Italy 10. 10.

Dipartimento di Fisica e Astronomia “G. Galilei”, Università di Padova, I-35131 Padova, Italy 11. 11.

National Aeronautics and Space Administration Goddard Space Flight Center, Greenbelt, MD 20771, USA 12. 12.

Department of Astronomy, University of Maryland, College Park, MD 20742, USA 13. 13.

Istituto Nazionale di Fisica Nucleare, Sezione di Pisa, I-56127 Pisa, Italy 14. 14.

Dipartimento di Fisica “M. Merlin” dell’Università e del Politecnico di Bari, I-70126 Bari, Italy 15. 15.

Istituto Nazionale di Fisica Nucleare, Sezione di Bari, I-70126 Bari, Italy 16. 16.

W. W. Hansen Experimental Physics Laboratory, Kavli Institute for Particle Astrophysics and Cosmology, Department of Physics and SLAC National Accelerator Laboratory, Stanford University, Stanford, CA 94305, USA 17. 17.

Istituto Nazionale di Fisica Nucleare, Sezione di Torino, I-10125 Torino, Italy 18. 18.

Dipartimento di Fisica, Università degli Studi di Torino, I-10125 Torino, Italy 19. 19.

Department of Physics and Astronomy, University of Padova, Vicolo Osservatorio 3, I-35122 Padova, Italy 20. 20.

Laboratoire Univers et Particules de Montpellier, Université Montpellier, Centre national de la recherche scientifique/Institut national de physique nucléaire de physique des particules, F-34095 Montpellier, France 21. 21.

Laboratoire Leprince-Ringuet, École polytechnique, Centre national de la recherche scientifique/Institut national de physique nucléaire de physique des particules, F-91128 Palaiseau, France 22. 22.

Center for Research and Exploration in Space Science and Technology (CRESST) and National Aeronautics and Space Administration Goddard Space Flight Center, Greenbelt, MD 20771, USA 23. 23.

Istituto Nazionale di Astrofisica-Istituto di Astrofisica Spaziale e Fisica Cosmica Milano, via E. Bassini 15, I-20133 Milano, Italy 24. 24.

Italian Space Agency, Via del Politecnico snc, 00133 Roma, Italy 25. 25.

Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA 26. 26.

Space Science Data Center - Agenzia Spaziale Italiana, Via del Politecnico, snc, I-00133, Roma, Italy 27. 27.

Istituto Nazionale di Fisica Nucleare, Sezione di Perugia, I-06123 Perugia, Italy 28. 28.

Department of Physics and Astronomy, Sonoma State University, Rohnert Park, CA 94928-3609, USA 29. 29.

Department of Physics, Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden 30. 30.

The Oskar Klein Centre for Cosmoparticle Physics, AlbaNova, SE-106 91 Stockholm, Sweden 31. 31.

Istituto Nazionale di Astrofisica Istituto di Radioastronomia, I-40129 Bologna, Italy 32. 32.

Dipartimento di Astronomia, Università di Bologna, I-40127 Bologna, Italy 33. 33.

Università Telematica Pegaso, Piazza Trieste e Trento, 48, I-80132 Napoli, Italy 34. 34.

Grupo de Altas Energías, Universidad Complutense de Madrid, E-28040 Madrid, Spain 35. 35.

Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, 91058 Erlangen, Germany 36. 36.

Circolo Astrofili Talmassons, I-33030 Campoformido (UD), Italy 37. 37.

Laboratoire de Physique et Chimie de l’Environnement et de l’Espace – Université d’Orléans / Centre national de la recherche scientifique, F-45071 Orléans Cedex 02, France 38. 38.

Station de radioastronomie de Nançay, Observatoire de Paris, Centre national de la recherche scientifique/Institut national des sciences de l’univers, F-18330 Nançay, France 39. 39.

The George Washington University, Department of Physics, 725 21st St, NW, Washington, DC 20052, USA 40. 40.

Max-Planck-Institut für Astrophysik, Postfach 1317, D-85741 Garching, Germany 41. 41.

Science Institute, University of Iceland, IS-107 Reykjavik, Iceland 42. 42.

Kungliga Tekniska högskolan Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden 43. 43.

Department of Physics, Kungliga Tekniska högskolan Royal Institute of Technology, AlbaNova, SE-106 91 Stockholm, Sweden 44. 44.

Centre d’Études Nucléaires de Bordeaux Gradignan, Institut national de physique nucléaire de physique des particules/Centre national de la recherche scientifique, Université Bordeaux 1, BP120, F-33175 Gradignan Cedex, France 45. 45.

Hiroshima Astrophysical Science Center, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan 46. 46.

Istituto Nazionale di Fisica Nucleare, Sezione di Roma “Tor Vergata”, I-00133 Roma, Italy 47. 47.

Department of Physics and Astronomy, University of Denver, Denver, CO 80208, USA 48. 48.

Max-Planck-Institut für Physik, D-80805 München, Germany 49. 49.

Osservatorio Astronomico di Trieste, Istituto Nazionale di Astrofisica, I-34143 Trieste, Italy 50. 50.

Department of Physics, University of Johannesburg, PO Box 524, Auckland Park 2006, South Africa 51. 51.

Institut für Astro- und Teilchenphysik and Institut für Theoretische Physik, Leopold-Franzens-Universität Innsbruck, A-6020 Innsbruck, Austria 52. 52.

Department of Physics, The University of Hong Kong, Pokfulam Road, Hong Kong, China 53. 53.

Laboratory for Space Research, The University of Hong Kong, Hong Kong, China 54. 54.

Nycb Real-Time Computing Inc., Lattingtown, NY 11560-1025, USA 55. 55.

Purdue University Northwest, Hammond, IN 46323, USA 56. 56.

Solar-Terrestrial Environment Laboratory, Nagoya University, Nagoya 464-8601, Japan 57. 57.

Centre national de la recherche scientifique, Institut de Recherche en Astrophysique et Planétologie, F-31028 Toulouse cedex 4, France 58. 58.

Galaxies Astrophysique des Hautes Énergies Cosmologie, Universit de Toulouse, Université Paul Sabatier-Observatoire midi-pyrénées, Institut de Recherche en astrophysique et planétologie, F-31400 Toulouse, France 59. 59.

Institute of Space Sciences (Consejo Superior de Investigaciones Científicas-Institut d’Estudis Espacials de Catalunya), Campus Universitat Autónoma de Barcelona, Carrer de Magrans s/n, E-08193 Barcelona, Spain 60. 60.

Institució Catalana de Recerca i Estudis Avançats (ICREA), E-08010 Barcelona, Spain 61. 61.

Istituto Nazionale di Fisica Nucleare-Istituto di Astrofisica Spaziale e Fisica Cosmica Bologna, via P. Gobetti 101, I-40129 Bologna, Italy 62. 62.

Dipartimento di Fisica, Università degli Studi di Perugia, I-06123 Perugia, Italy 63. 63.

INTErnational Gamma-Ray Astrophysics Laboratory Science Data Centre, CH-1290 Versoix, Switzerland 64. 64.

California State University, Los Angeles, Department of Physics and Astronomy, Los Angeles, CA 90032, USA 65. 65.

Praxis Inc., Alexandria, VA 22303, resident at Naval Research Laboratory, Washington, DC 20375, USA 66. 66.

Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, and Università di Trieste, I-34127 Trieste, Italy 67. 67.

Center for Astrophysics and Cosmology, University of Nova Gorica, Nova Gorica, Slovenia 68. $\dagger$

majello@g.clemson.edu, helgason@hi.is, vpaliya@g.clemson.edu,

justin.finke@nrl.navy.mil, abhishd@g.clemson.edu, alberto@gae.ucm.es

Materials and Methods

Sample Selection and Data Analysis

Our sample is selected starting from the objects reported in the third catalog of active galactic nuclei detected by the LAT, 3LAC,  (?). We exclude all the blazars reported there with a double association and those lacking a redshift measurement. Most redshift measurements for BL Lacs reported in 3LAC come from  (?). For each source we assess the significance of the detection (between 1 GeV and 1 TeV) defining a test statistics ($TS$) as $TS=2\Delta\log\mathscr{L}$, where $\mathscr{L}$ represents the likelihood function between models with and without the source of interest. We use this to exclude all sources that have a $TS<25$ in this analysis. Our final sample comprises 419 FSRQs and 320 BL Lacs distributed (see Figure S1) between a redshift of 0.03 and 3.1.

The analysis relies on 101 months (Aug. 2008 to Jan. 2017) of Pass 8 (P8) class ‘SOURCE’ photons detected by the LAT between 1 GeV and 1 TeV. This dataset was filtered to eliminate times when the spacecraft was over the South Atlantic Anomaly and to remove photons detected at angles larger than 100*∘* from the zenith. For the analysis of each source we use photons within 15*∘* of the source position (region of interest, ROI). For each ROI we define a sky model that comprises the diffuse Galactic  (?) and extragalactic emission  (?) as well as the emission from background sources in the ROI. The latter includes sources detected in the third Fermi-LAT catalog, 3FGL,  (?) as well as any new source that is detected because of the additional exposure (with respect to the 3FGL) used here. These sources are found generating a $TS$ map and identified as excesses above a $TS=25$ threshold and added to the sky model with a power-law spectrum. The LAT ‘P8R2$\_$SOURCE$\_$V6’ instrumental response function (IRF) and a binned likelihood method are used to fit the sky model to the data.

Intrinsic Spectra of Blazars

To capture the intrinsic curvature in the spectra of blazars we adopt the following strategy that has been optimized using simulations prior to the data analysis (see below). The data are fitted only to a maximum energy up to which the attenuation of the EBL is negligible. This is defined as the energy at which the optical depth \tau_{\gamma\gamma}$$<0.1 for the model of  (?). However, we tested that our analysis is robust against changes of EBL model used to define this maximum energy and changes to the threshold (from e.g. \tau_{\gamma\gamma}$$<0.1 to \tau_{\gamma\gamma}$$<0.05). The optical depth decreases sharply in this regime where not many EBL photons are expected due to a characteristic drop-off at the Lyman-limit (13.6 eV). Our baseline model for the intrinsic blazar spectrum is a log-parabola:

[TABLE]

where $N_{0}$ (the normalization), $\alpha$ (photon index), $\beta$ (curvature) are all free parameters and $E_{b}$ is a scaling energy. We also test whether an exponential power law could be a better representation of the blazar spectrum and this is defined as:

[TABLE]

where $E_{c}$ (cut off energy) and $\gamma_{1}$ (the exponential index) are all free parameters. Smoothed broken power laws and broken power laws were also tested, but they were never found to describe the blazar intrinsic spectrum better than the two models reported above in the energy range used in this work.

When testing the exponential cut-off model we perform a first fit with $\gamma_{1}$ fixed at 0.5 (justified from the observations of hundreds of FSRQs, see  (?)) and then another fit leaving $\gamma_{1}$ free to vary. We define two TS of curvature $TS_{\rm c,1}$ and $TS_{\rm c,2}$ as follows:

[TABLE]

where $\log L_{{\rm exp},\gamma_{1}=0.5}$ and $\log L_{{\rm exp},\gamma_{1}={\rm free}}$ are the log-likelihoods derived using the exponential cut-off model with $\gamma_{1}=0.5$ and $\gamma_{1}$ free to vary respectively and $\log L_{\log-{\rm parabola}}$ is the log-likelihood of the log-parabola model.

We adopt the criteria reported in Table S1 to choose the model used to describe each blazar’s intrinsic spectrum. In order to avoid convergence problems, in the analysis presented above, the exponential index $\gamma_{1}$ remains fixed at either 0.5 or the best-fitting value found at this step. The median of the distribution of fitted $\gamma_{1}$ values is $\approx$0.5.

Analysis

Results for Blazars

Once the choice of the intrinsic spectrum for the sources has been made, the analysis reverts to using the full, 1 GeV–1 TeV, energy band and the modeled spectra of all sources include the EBL attenuation as reported in Equation 1, where $b$ is a parameter, common to all sources, that is varied to fit the EBL model prediction to the data. A $b=1$ would mean that the EBL model predictions are in agreement with the LAT data, while a $b=0$ would imply that there is no evidence for attenuation due to absorption by EBL photons in the spectra of blazars.

Because of the complexity of the problem, the $b$ parameter is not optimized in one stage. Instead, for each source we scan the likelihood function in very small steps of $b$ creating a profile likelihood. In this process, the parameters of the diffuse emission, those of the brightest sources, and those of the source of interest (except $\gamma_{1}$) are all left free to vary. For each source, the best-fitting $b$ value is the one that maximizes the log-likelihood. A $TS$ of the detection of the EBL can be built comparing the log-likelihood values at the best-fitting $b=b_{0}$ and at $b=0$ as $TS_{\rm EBL}=2[\log L(b_{0})-\log L(b=0)]$. Because log-likelihoods (and thus $TS$) are additive, we can determine the $b$ value that maximizes the global (for all sources) likelihood and produces the largest $TS_{\rm EBL}$. In Figure S2, we plot the $TS$ profile, as a function of $b$, for all sources (and separately for BL Lacs and FSRQs) for the model of  (?). A $b$=1.03 improves the fit by a $TS$ of $\sim$300, which corresponds to $\sim$17 $\sigma$ for one degree of freedom. We note that the spectral evolution of the blazar class with redshift has a negligible effect on this analysis, as apparent from Figure S2, which shows that the level of EBL measured using (mostly) hard-spectrum BL Lacs is in very good agreement with that found using soft-spectrum FSRQs. As an additional test we report the values of the $b$ parameter for the model of  (?) for BL Lacs with a synchrotron peak frequency $>10^{16}$ Hz (called HSPs) and for the remaining BL Lacs; these are respectively $b_{\rm HSPs}=0.98^{+0.09}_{-0.13}$ ($TS_{\rm EBL}$=125.8) and $b_{\rm rest}=0.86^{+0.16}_{-0.10}$ ($TS_{\rm EBL}=45.1$). These highlight once more that there is no bias in the level of the EBL due to the spectral evolution of the blazar class.

One can also measure the compatibility of a model prediction with the Fermi-LAT data defining a $TS$ as $TS_{b=1}=2[\log L(b_{0})-\log L(b=1)]$. By definition a large $TS_{b=1}$ implies that the model predictions are in tension with the Fermi-LAT data; this typically happens when the model predicts a larger-than-observed attenuation. Table S2 shows the results of our analysis for some of the models available in the literature that have not been found in tension with previous $\gamma$-ray data. The table shows that the high model of  (?) and the best-fitting model of  (?) are excluded. Moreover, the models of  (?) and  (?) are found in tension at the $\sim 3$ $\sigma$ level with the Fermi-LAT observations. All these models predict a larger optical-UV intensity of the EBL than the models compatible with the LAT data.

The optical depth as a function of energy and redshift can be measured by repeating the above procedure (i.e., renormalizing the optical depth predicted by a model), but in small energy and redshift bins. In this process, the uncertainty due to the small disagreement between different EBL models, about the shape of the optical depth curve within any given bin, has been included in the final uncertainty of the optical depth. The final uncertainty includes also the 10 % systematic uncertainty discussed below. The redshift bins were chosen so that similar values of $TS_{\rm EBL}$ were obtained in all the bins. Figure S3 shows measurements of the optical depth $\tau_{\gamma\gamma}$ due to EBL absorption in different redshift and energy bins. It is apparent from the figure that most of the constraining power is around $\tau_{\gamma\gamma}\approx 1$. Formally the $\tau_{\gamma\gamma}(E,z)=1$ value marks the cosmic $\gamma$-ray horizon, i.e., the energy above which our Universe becomes opaque to $\gamma$ rays for a given redshift  (?, ?). The energy at which $\tau_{\gamma\gamma}(E,z)=1$ at any redshift can be found by renormalizing any EBL model to fit the data presented in Figure S3 and propagating the (statistical plus systematic) uncertainties. Figure 1 shows that Fermi LAT maps the horizon position with energy from low ($z\approx 0.03$) to high ($z\approx 3.1$) redshift. Figure 1 also shows the highest-energy photons detected from the blazars in our sample.

GRB 080916C

In order to constrain the EBL and SFH to the highest possible redshifts, we complement the blazar sample with a single gamma-ray burst (GRB), GRB 080916C, detected by Fermi LAT at $z$=4.35  (?, ?). This was an extremely luminous event, whose hard spectrum has already produced constraints on the SFH at high redshift  (?, ?). With respect to previous works, the release of the Pass 8 event-level analysis has allowed us to recover more high-energy photons, particularly one at 27.4 GeV, $\approx$146 GeV in the source frame, from GRB 080916C  (?).

The analysis is similar to the one reported by  (?). Transient-class photons between 0.1 GeV and 100 GeV were downloaded around a 10*∘* position from the burst and from the time of the GRB until 1775.9 s later. Photons detected at a zenith angle greater than 105*∘* were removed. The ROI model consists of the burst, the Galactic and isotropic templates  (?). We rely on the ‘P8R2_TRANSIENT020’ IRF.

The source intrinsic spectrum is represented (and fitted up to 10 GeV) by a single power law (with a photon index of 2.25$\pm 0.06$) employing a time-averaged analysis. No curvature is observed in the Fermi-LAT spectrum. A time-resolved analysis does not yield any difference for this work  (?). We produce a 95 % upper limit on the optical depth by adopting the same method as described above for blazars. This upper limit is $\tau_{\gamma\gamma}<0.46$ at an energy of $\approx$17 GeV and $z$=4.35 and it does not depend on the EBL model used to derive it. This upper limit is a factor of two lower than that used by  (?). This is due to the additional photons detected at $>$10 GeV and particularly to the 27.4 GeV photon. The probability that this photon belongs to the background, rather than to GRB 080916C, is only 5$\times 10^{-5}$. A so-called ‘maximally conservative upper limit’ based on the assumption that the intrinsic spectrum cannot be harder than a power law with an index $\Gamma=1.5$ is even more constraining; however, it is not adopted here. The spectrum of GRB 080916C used here is shown in Figure S4.

Tests and Simulations

Simulations of Blazar SEDs

The analysis chain described in the previous section has first been tested and optimized on Monte Carlo simulation of synthetic spectral-energy distributions (SEDs) of blazars with properties matching those of blazars observed by Fermi LAT.

The SEDs are generated from physical models of blazars’ emission that include synchrotron and synchrotron self-Compton as well as (for FSRQs) external Compton scattering and were generated with the numerical code presented in  (?) and  (?). These SEDs reproduce the range of peak frequencies very well (for both the synchrotron and $\gamma$-ray components), including peak curvatures and $\gamma$-ray photon indices observed in both Fermi-LAT BL Lacs and FSRQs. They include all known effects that contribute to determining the curvature of the intrinsic $\gamma$-ray spectrum of Fermi blazars. The crucial transition from the Thomson to the Klein-Nishina cross section as well as (mostly important for FSRQs) absorption within the broad line region (for different distances of the emission region from the broad line region) are all taken into account and contribute to determine the shape of the blazars’ spectra at high energy.

These SEDs, attenuated by the EBL for a range of redshifts similar to those of Figure S1, are then used to simulate LAT observations of these synthetic sources and have been used to optimize the analysis set-up presented above. In particular, the values of the minimum energy ($E_{\rm min}$=1 GeV) and those of $TS_{\rm c,1}$ and $TS_{\rm c,2}$ have been derived from the analysis of simulations. Figure S5 shows that the analysis chain employed in this work can effectively recover the simulated level of EBL.

Variability

Blazars are inherently variable objects (at all wavelengths) and variability may bias or complicate the measurement of the EBL attenuation. In order to cope with blazars’ variability as much as possible, a time-resolved analysis is performed for all sources that are found variable at $>10$ GeV in the recent third catalog of hard sources, 3FHL,  (?). We rely on the time bins derived by the Bayesian block analysis presented in 3FHL as these are times when the sources were found to alter their state at $>10$ GeV, which is the relevant energy range for detecting the EBL attenuation. In each time bin, the criteria reported in Table S1 are used to determine the best intrinsic spectral model. Because for a given source time-resolved spectra can be treated as independent observations, their contribution to the $TS_{\rm EBL}$ has been summed to the one of the remainder of the sample. Figure S2 includes the contribution from variable and non-variable sources. The level of EBL as determined from the variable sources alone is found to be in good agreement with the rest of the sample. Figure S6 shows the $TS_{\rm EBL}$ as a function of the $b$ parameter (for the model of  (?)) for 4 variable BL Lacs and FSRQs and how that compares to the result of a time-averaged analysis. In general, we find a time-averaged analysis works well for objects which vary primarily in flux, while a time-resolved analysis is required for all those objects experiencing also spectral variability (see right versus left plots in Figure S6). Finally, we used the Fermi All-sky Variability Analysis tool  (?, ?) to search for significant residual spectral variability within Bayesian blocks, but none could be found.

Systematic Uncertainties

In order to gauge the systematic uncertainties of this analysis we have performed the tests reported below:

• •

Instead of using a variable maximum energy up to which to fit the intrinsic spectrum (chosen to be the energy corresponding to $\tau_{\gamma\gamma}<0.1$ for the model of  (?)), we use a constant maximum energy of 10 GeV for all sources. Repeating the entire analysis we find $b=1.09\pm 0.08$ in agreement with the results presented in the main text.

• •

A similar result as above has been obtained using a maximum energy, to measure the intrinsic spectrum, defined as that obtained when \tau_{\gamma\gamma}$$<0.05 ($b=1.07\pm 0.08$ for the model of  (?)).

• •

We use the IRF bracketing method as described in  (?). By deriving two different sets of IRFs and repeating the entire analysis we find that the systematic uncertainty in the optical depth $\tau_{\gamma\gamma}$ is of the order $\sim$7 %.

The results presented in the above sections are fully confirmed and the systematic uncertainty on the optical depth $\tau_{\gamma\gamma}$ due to changing the energy threshold to characterize the intrinsic spectrum and IRF are, together, $\lesssim$10 %. A systematic uncertainty of $0.1\times\tau_{\gamma\gamma}$ (added in quadrature) has been included in the uncertainties reported in Figures 1 and S3 and propagated to all results that use those data.

Reconstructing the evolving EBL

The optical depth for a $\gamma$ ray of observed energy $E_{\gamma}$ originating in a source at redshift $z_{s}$ is related to the evolving number density of EBL photons, $n_{\rm EBL}(\epsilon,z)$,  (?, ?, ?):

[TABLE]

where the rest-frame energy of $\gamma$ rays and EBL photons are denoted by $\epsilon_{\gamma}=E_{\gamma}(1+z_{s})$ and $\epsilon_{\rm EBL}=E_{\rm EBL}(1+z_{s})$ respectively, $\mu=\cos{\theta}$ denotes the angle of incidence between the two photons, and $|dt/dz|^{-1}=H_{0}(1+z)\sqrt{\Omega_{M}(1+z)^{3}+\Omega_{\Lambda}}$, where the Hubble and the cosmological parameters are $H_{0}=70$ km s*-1* Mpc*-1*, $\Omega_{M}$=0.3, and $\Omega_{\lambda}=$0.7. The cross section for the photon-photon interaction appearing in the last integral in Equation S5 is

[TABLE]

with

[TABLE]

where $m_{e}c^{2}$ is the electron rest mass. In other words, for a given cosmology, the SED and evolution of the EBL uniquely specify the optical depth at all redshifts. Conversely, we can use the measured optical depths $\tau_{\gamma\gamma}(E_{\gamma},z)$ to reconstruct $n_{\rm EBL}(\epsilon,z)$.

The physical properties of galaxies, such as star-formation rate, stellar mass and metallicity, are encoded in their SED. Rather than the EBL, which is accumulated over cosmic time, it is more informative to study the instantaneous SED of the galaxy population as a whole i.e., the cosmic emissivity. The buildup of the EBL is related to the volume emissivity $j(\epsilon,z)$ (or equivalently, luminosity density) via:

[TABLE]

where $\epsilon^{\prime}=\epsilon(1+\bar{z})$ is the rest-frame energy at $\bar{z}$. The EBL spectral intensity (see Figure 2) can be found from the number density by $\epsilon I(\epsilon,z)=\frac{c}{4\pi}\epsilon^{2}n_{\rm EBL}(\epsilon,z)$.

Model for the Cosmic Emissivity

Our task is to reconstruct $j(\epsilon,z)$ based on the measured optical depths reported in Figure S3 without making assumptions on galaxy properties or their stellar population. We represent $j(\lambda)$ as the sum of several log-normal templates with a fixed peak position:

[TABLE]

where we fix $\sigma=0.2$, $\lambda_{i}=[0.17,0.92,2.2,8.0]\,{\rm\upmu m}$ and the amplitudes $a_{i}$ are left free to vary. We find that four log-normal templates allow for a sufficiently flexible spectral shape from UV to the mid-IR. A Lyman-break is imposed by cutting off the spectrum at $\epsilon>13.6$ eV where neutral hydrogen becomes opaque. We have chosen the fixed locations ($\lambda_{i}$) and width ($\sigma$) of the templates such that common features in galaxy SEDs, a flat far-UV continuum and a 4000Å/Balmer break, are easily captured. Each template is allowed to evolve independently with redshift based on a function similar to the SFH parametrization of  (?) leading to the full expression

[TABLE]

At each of the fixed wavelengths $\lambda_{i}$, one parameter controls the amplitude, $a_{i}$, and three control the evolution, $b_{i}$, $c_{i}$ and $d_{i}$, yielding a total of $4\times 4=16$ free parameters. The number of optical depth data points is 60.

To explore the sensitivity to different functional forms for the evolution, we also test the parametrization from  (?):

[TABLE]

with free parameters $a_{i}$, $b_{i}$, $c_{i}$, and $d_{i}$ which we display alongside our main results for the SFH in Figure S11.

Markov Chain Monte Carlo: Setup

We use the MCMC code emcee  (?), a Python implementation of an affine invariant MCMC ensemble sampler  (?), to constrain the parameters controlling the emissivity. The likelihood function is estimated as ${\mathscr{L}}\propto\exp{(-\chi^{2})}$ where $\chi^{2}$ is given by

[TABLE]

where there are $N$ energy ($E_{i}$) bins, $M$ redshift ($z_{j}$) bins, $\tau_{\rm data}(E_{i},z_{j})$ is the measured absorption optical depth presented in Figure S3, $\tau_{\rm model}(\vec{\theta}|E_{i},z_{j})$ is the model absorption optical depth with parameters $\vec{\theta}$, and $\sigma_{i,j}$ is the (statistical plus systematic) uncertainty on the absorption optical depth measurements.

We choose flat priors on all parameters $\log{a_{i}}/{\rm(ergs^{-1}Mpc^{-3}Hz^{-1})}=[22,32]$, $b_{i}=[-2,10]$, $c_{i}=[1,7]$, $d_{i}=[0,20]$. We restrict the range of the evolution parameter controlling the location of the peak (or curvature) to $c_{i}=[1,7]$ since our dataset has limited constraining power for shape changes at redshifts much larger than our sample coverage ($0<z<4$). Note however, that this does not force the presence of a peak and a turnover as the function reduces to a power-law $\propto(1+z)^{b_{i}}$ when $d\rightarrow 0$.

With the emissivity specified as a function of wavelength and redshift, we calculate the resulting EBL and optical depth according to Equations S5–S7 respectively for each proposed step in the MCMC. Each calculation of $\tau_{\gamma\gamma}$ involves integrating over wavelength, redshift and angle of incidence, but we only require $\tau_{\gamma\gamma}$ at six energies, for every redshift, making it computationally manageable. Our final results are based on MCMC chains from 120 walkers exploring the parameter space in 10,000 steps each. This results in 1,140,000 steps after a burn-in of 500 steps for each walker.

Results and Validation

In Figure S7 we display the 68% and 95% confidence regions for the total cosmic emissivity in several redshift bins. The Fermi-LAT dataset is tightly constraining at UV, optical and, at low-$z$, also near-IR wavelengths. The confidence regions get broader towards mid-IR wavelengths due to the energy range of Fermi LAT being limited to $<1$ TeV. Figure  S7 also shows that the Fermi-LAT dataset provides the tightest constraints around $z\simeq 0.5-1.5$ as the opacity increases for larger distances traveled. At $z\gtrsim 2$ we are limited by the number of bright blazars with substantial emission above $\sim$10 GeV in our sample. Comparing our cosmic emissivity with measurements of integrated galaxy luminosity functions shows that our results are in good overall agreement across the wavelength range. This implies that the bulk of the EBL is already accounted for by galaxy surveys.

We have validated this reconstruction method by creating ten sets of fake $\tau(E,z)$ data points in the same energy and redshift bins, and possessing the same fractional uncertainties, as the original dataset. The simulated datasets are generated by drawing random sets of values for the 16 parameters of Equation S9, within their assumed priors, and calculating the optical depths at each energy and redshift bin. Our reconstruction recovers the fake EBL in all cases within the derived 2 $\sigma$ uncertainty region. We see no evidence for a systematic over- or underestimation of the emissivity at any particular wavelength. The largest offsets arise at $>1{\rm\upmu m}$ where the dataset is less constraining.

Comparison with EBL models and data

The literature offers different approaches to estimate the evolving EBL . There are methodologies that are observationally motivated  (?, ?, ?, ?, ?), physically motivated  (?, ?, ?, ?), and based on semi-analytical models of galaxy formation such as  (?, ?). Typically these models are constructed in such a way that the lower redshifts and, in general, the optical/near-IR peak are better constrained. Figure S8 shows our reconstructed EBL spectral intensities in comparison with some of the models.

The reconstructed EBL follows galaxy counts  (?, ?) leaving little room for substantial contributions not resolved by deep galaxy surveys. This is in conflict with several direct measurements of EBL (which may be contaminated by foregrounds, see e.g.  (?, ?, ?)) and in tension with some models proposed to explain the anisotropies measured in diffuse light  (?).

Relative to EBL models, in the local Universe, we find that our estimate roughly follows the median of existing models. The models by  (?),  (?),  (?),  (?), and  (?) reproduce our results quite well. The fiducial model by  (?) tends to follow the upper region of our 1$\sigma$ band. Finally, the baseline model by  (?) provides too much UV and too little near-IR.

The strategy of using the observation of $\gamma$-ray photons to derive constraints on the background has been used extensively in recent years. Early attempts, characterized by scarcer $\gamma$-ray data, only allowed intensity upper limits as a consequence of assumptions on the intrinsic spectra of blazars  (?, ?, ?). These results were followed by more sophisticated approaches based on more and better data that allowed the EBL detection and study, both with the LAT at somewhat lower energies and thus, larger redshifts  (?, ?, ?, ?), and with Imaging Atmospheric Cherenkov Telescopes (IACTs). These results from IACTs mostly constrain the local Universe  (?, ?, ?), although the MAGIC collaboration also probed the EBL at $z\sim 1$ with the detection of two blazars  (?, ?). Notably, our derived EBL at $z=0$ is even closer to the integrated counts compared to previous $\gamma$-ray derived EBL measurements.

Table S4 reports the EBL spectral intensities at several redshifts as displayed in Figure 2. Using Equation S5, we can calculate the optical depth curves as a function of energy and redshift from our reconstructed EBL. These can be used to correct spectra of cosmological $\gamma$-ray sources for EBL absorption in order to study physical properties of the source and are provided \hrefhttps://figshare.com/s/14f943002230d69a4afdonline.

Implications for the high-$z$ Universe

In principle, there is always some constraining power beyond the maximum redshift of the sample of $\gamma$-ray sources. This is because the $\gamma$ rays coming from $z_{\rm max}$ start interacting with EBL photons which were built up at still earlier times and the rate of the interactions is related to $n_{\rm EBL}\propto(1+z)^{3}$. Focusing on the UV, which is important for cosmic re-ionization, Figure S9 suggests rather minimal UV emissivity at $z>4$ with respect to measurements from Lyman-break galaxy surveys. However, it is possible that the confidence regions at these redshifts may be artificially narrow due to the lack of flexibility in the parameterized shape of the evolution. In order to test the robustness of the constraints at high-$z$ we have re-run the MCMC and included an additional term in Equation S9:

[TABLE]

centered at $z_{0}=6$ with $\sigma=0.08$. We find that, while this reveals a relatively unconstrained lower limit for the UV emissivity, the upper limit remains robust at $\sim 3.2(5.3)\times 10^{26}$ erg s*-1* Mpc*-3* Hz*-1*, 1$\sigma$(2$\sigma$), at $z=5-6$. In Figure 4, we compare this value with the integrated UV luminosity functions from the Hubble Frontier Fields (HFF) program, which targets extremely faint galaxies behind strong gravitational lenses, reaching $M_{\rm AB}\sim-13$. Some HFF analyses have found evidence for a turnover in the steep faint-end of the luminosity function (LF)  (?, ?), whereas others do not see such a feature  (?, ?). The conflicting results at $M_{\rm AB}\gtrsim-15$ could be due to uncertainties in the magnification factor determined by lens models.

Our constraints limit how far a steep faint-end slope can be extrapolated. In fact, Figure 4 shows that the emissivities from the integrated UV luminosity functions are already close to our derived upper limits, but are all compatible within 2 $\sigma$. They favor a turnover of the LF at $M_{\rm AB}\sim-14$ in agreement with  (?) and  (?). The UV emissivity implied by  (?) for example (with no turnover), would reach the 2$\sigma$ upper limit if extrapolated to $M_{\rm AB}\sim-10$.

In Figure 4, we also show the UV emissivity necessary to sustain a reionized Universe at $z=6$. The required emissivity (at 0.15 ${\rm\upmu m}$) can be shown to be  (?):

[TABLE]

Here, $\Omega_{b}$ is the cosmic baryon density, $h_{70}$ is Hubble parameter in units of 70 ${\rm km~{}s^{-1}Mpc^{-1}}$, $C$ is the clumping factor of ionized hydrogen and $f_{\rm esc}$ is the mean escape fraction of ionizing photons. The parameter $\epsilon_{53}$ is the number of Lyman continuum photons per unit of forming stellar mass in units of $10^{53}$ photons $\cdot{\rm s^{-1}(M_{\odot}\cdot yr^{-1})^{-1}}$. For this we follow  (?) exploring values of $\epsilon_{53}$ based on stellar population synthesis models assuming a Salpeter IMF and a constant star formation rate. The width of the grey regions in Figure 4 correspond to the range $0.9<\epsilon_{53}<1.4$ when the metallicity is varied from $0.02Z_{\odot}$ to $1.0Z_{\odot}$ (where $Z_{\odot}$ is the solar metallicity). We display the emissivity for a reasonable assumption of $C/f_{\rm esc}=30$, showing that our constraints accommodate a scenario in which the Universe is reionized at $z=6$.

Our constraints at $z>4$ come almost entirely from GRB 080916C which provides a strong upper limit to the optical depth at $z=4.35$ whereas the blazar sample alone ($z<3.1$) has lower constraining power. This is a benefit of detecting more high-$z$ $\gamma$-ray sources as probes of the epoch of re-ionization  (?).

The Star-Formation History

We derive the SFH from our constraints on the far-UV emissivity in a similar manner to galaxy surveys that measure the rest-frame UV emission  (?, ?, ?, ?, ?). The conversion into star-formation rate (SFR) requires two assumptions: i) the amount of UV emission expected per unit SFR, $\mathscr{K}_{\rm UV}$, which is dictated by the initial mass function (IMF) of choice, and ii) the mean dust extinction within the host galaxies, $A_{V}$, since photons become a part of the EBL only if they escape their progenitor galaxies. For the former quantity, we assume $\mathscr{K}_{\rm UV}=7.25\times 10^{-29}{\rm M_{\odot}~{}yr^{-1}~{}erg^{-1}~{}s~{}Hz}$ which is consistent with a Chabrier IMF  (?). Our results on the SFH can be re-scaled by constant factor of 1.6 to represent a Salpeter IMF  (?).

For the dust extinction correction, we rely on measured values of the mean $A_{V}$ from the literature and fit its evolution with redshift using the following parametrization: $A_{V}\propto\frac{(1+z)^{f}}{1+\left(\frac{1+z}{c}\right)^{d}}$. The result is shown in Figure S10. The measured values of $A_{V}$ are based on different methods. For instance, these come from: measured UV continuum slopes  (?, ?), stellar population synthesis SED fitting  (?, ?) and comparison of the integrated UV and IR luminosity functions  (?, ?). We choose to use only those data that are measured from a large sample where robust uncertainty estimation is provided. Studies that assume or estimate values of $A_{V}$ do not contribute to the fit but are shown in Figure S10 for reference. We obtain the evolving extinction:

[TABLE]

The SFH is then calculated as:

[TABLE]

where $j_{\rm UV}(z)$ is our reconstructed emissivity at 1600Å.

The confidence regions for the cosmic SFH are shown in Figure 3 in the main paper along with data points from UV-derived measurements  (?). We also display the same result in Figure S11 showing data from various studies using different tracers of SFR, including limits from $\gamma$-ray constraints of the EBL  (?). At low and intermediate redshifts, our results are in good agreement with (albeit a little bit above) independent measurements from galaxy surveys. At $z>3$, our results are in agreement, within the uncertainties, but favor a rather low SFH. As discussed in the previous subsection, this is primarily driven by GRB 080916C. More importantly, because the SFH derived from $\gamma$-ray absorption complements traditional methods that probe the SFH from sources resolved in surveys, our results imply that the bulk of star formation across cosmic time is already accounted for by surveys.

Stellar Population Model Method

The cosmic SFH, $\dot{\rho}(z)$, is the starting point in the EBL model of  (?, ?) making it a useful model for further exploration of the parameter space that is made possible by the $\gamma$-ray optical depth data. The model assumes that stars emit as blackbodies, with their temperatures, luminosities, and time evolution determined from formulae given by  (?). The radiation emitted by stars is convolved with an IMF and star-formation rate density parameterization to get the luminosity density $j(\epsilon;z)$. The fraction of light that escapes dust extinction ($f_{\rm esc,dust}$) is based on the extinction curve from  (?), which was derived from a fit to the luminosity density data in the local Universe. We let the dust extinction evolve according to Equation S14. The infrared portion of the EBL is computed assuming that all the energy absorbed by dust is re-radiated in the infrared. The SFH and IMF model parameters were chosen to reproduce the luminosity density data available at the time. Once $j(\epsilon;z)$ is calculated, the EBL number density and $\gamma$-ray absorption optical depth are computed from Equations S5–S7 above.

Using the methodology of  (?), we have performed an MCMC fit to the $\gamma$-ray optical depth data. We parameterize the SFH and let the parameters vary, calculating the resulting EBL and optical depths in each step. A similar MCMC model fit, but limited to $z\geq 2$, was done by  (?) to the earlier EBL absorption data from  (?). We use the standard parameterization for the SFH (Equation S9), but also consider evolution according to Equation S10. The SFH result from our MCMC fits, reported in Figure 3 (as the green confidence region), are consistent with the SFH used for the “model C” of  (?), which relied on the  (?) parametrization with free parameters given by  (?), at all values of $z$, but the confidence interval is particularly narrow up to $z\leq 2.5$. Table S5 reports the values of the SFH obtained from both methods as displayed in Figure 3.

The contribution of Active Galactic Nuclei

The methods employed here to derive the SFH of the Universe rest on the assumption that most of the EBL is the product of stellar emission. While this is almost certainly true in the IR and optical bands, there could be a non-negligible contribution of active galactic nuclei (AGN) to the global UV background  (?). We estimate this contribution by considering measurements of the integrated quasar luminosity function. A fitting formula for the resulting emissivity as a function of redshift was provided by  (?) at a rest-frame wavelength of 912Å. We convert this to 0.16 $\upmu$m (the wavelength used to estimate our SFH) adopting the same power law spectrum $\propto\lambda^{-0.61}$  (?) and show the ratio of the AGN to total (as estimated in this work) emissivities at 0.16 $\upmu$m in Figure S12. This shows that the contribution from known AGN is no more than a $few$ percent. This is in agreement with the recent estimate of the AGN contribution to the EBL  (?).