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

See allHide authors and affiliations

Science  30 Nov 2018:
Vol. 362, Issue 6418, pp. 1031-1034
DOI: 10.1126/science.aat8123

Gamma rays reveal the Universe's history

How many stars have formed in the Universe, and when did they do so? These fundamental questions are difficult to answer because there are systematic uncertainties in converting the light we observe into the total mass of stars in galaxies. The Fermi-LAT Collaboration addressed these questions by exploiting the way that gamma rays from distant blazars propagate through intergalactic space, which depends on the total amount of light emitted by all galaxies. The collaboration found that star formation peaked about 3 billion years after the Big Bang (see the Perspective by Prandini). Although this is similar to previous estimates from optical and infrared observations, the results provide valuable confirmation because they should be affected by different systematic effects.

Science, this issue p. 1031; see also p. 995


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 measured this attenuation using 739 active galaxies and one gamma-ray burst detected by the Fermi Large Area Telescope. This allowed 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 ~ 2. Upper limits of the EBL at the epoch of reionization suggest a turnover in the abundance of faint galaxies at z ~ 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)—that is, the stellar birth rate as a function of the Universe’s age—summarizes the history of stellar formation since the Big Bang (1). 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), 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 (24). 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 thereby estimating the star formation associated with those galaxies. This difficulty becomes particularly relevant within the first billion years after the Big Bang when a large population of faint, still undetected, galaxies existed (5). These galaxies are expected to have driven the reionization 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 1 billion years from the present age) star formation measured using spaceborne UV observatories is based on surveys extending over small solid angles (6) and is 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 (79). Accurate measurements of this diffuse all-sky background at UV to IR wavelengths, and particularly its buildup over time, have only just become possible (10).

We present an alternative approach to measuring the SFH, based on the attenuation that the EBL produces in the γ-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 γγ → e+e), effectively being absorbed as a result of the interaction (11). Above a given threshold energy, the attenuation experienced by every γ-ray source at a similar distance depends on the number density of the EBL target photons integrated along the line of sight; observations of γ-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 analyzed γ-ray photons detected by the Large Area Telescope (LAT) instrument onboard 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 supermassive 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 to 11.6 billion years according to the standard cosmological model (12). We performed a likelihood analysis to find the EBL attenuation experienced by all blazars while simultaneously optimizing the spectral parameters independently for each blazar (13). This can be accomplished individually for each source by defining a region of interest that comprises all γ rays detected within 15° of the source position and creating a sky model that includes all sources of γ 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 (14) and extrapolated to higher energy, where the γ rays are expected to be attenuated by the EBL.

Potential EBL absorption is added to the fitted spectra as follows:Embedded Image(1)where (dN/dE)obs and (dN/dE)int are the observed and intrinsic blazar spectra, respectively; τγγ(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 that no EBL attenuation is present in the spectra of blazars, whereas b ≈ 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 ≳16 standard deviations (σ) 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 γ rays, known as the cosmic γ-ray horizon (Fig. 1). With the optical depths measured in six energy bins (10 to 1000 GeV) across 12 redshift bins (14), we are able to reconstruct the intensity of the EBL at different epochs (Fig. 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 ~90% of cosmic time.

Fig. 1 The cosmic γ-ray horizon.

Measurement of the cosmic γ-ray horizon (τγγ = 1, i.e., the point after which the Universe becomes opaque to γ rays) as a function of redshift (blue stars). Horizontal and vertical blue lines represent the redshift bin size and the uncertainty on the energy, respectively. Also shown are predictions from three different EBL models (29, 37, 38) for comparison. The gray points show the highest-energy photon (HEP) detected from each blazar considered in this work. Gyr, billions of years.

Fig. 2 The spectral intensity of the EBL at redshifts 0, 1, 2, and 3.

(A) Today’s Universe (z = 0). The blue area shows the 1σ confidence regions based on the reconstructed cosmic emissivity (14). Data from other γ ray–based measurements are shown with orange symbols (3942); integrated galaxy counts are displayed with green symbols (1520). Horizontal yellow lines represent the wavelength range; vertical yellow and green lines represent the uncertainty on the intensity of the EBL. The downward-pointing arrows show upper limits on the intensity of the EBL from (42). (B to D) At higher redshifts (z = 1, 2, and 3, respectively), the EBL is shown in physical coordinates. Figure S8 includes a more complete set of measurements from the literature.

At z = 0, the energy spectrum of the EBL is close to the one inferred by resolving individual galaxies in deep fields (15). 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 Fig. 2). The intensity of the UV background in the local Universe remains uncertain, with independent studies reporting differing values (1618). Our determination of Embedded Image nW m–2 sr–1 [1σ (2σ)] at 0.2 μm favors an intermediate UV intensity in agreement with (18). In the NIR, our measurement of Embedded Image nW m–2 sr–1 [1σ (2σ)] at 1.4 μm is consistent with integrated galaxy counts (19, 20), leaving little room for additional components, contrary to some suggestions (21, 22). This notably includes contributions from stars that have been stripped from galaxies, as the technique presented here is sensitive to all photons (23, 24).

At any epoch, the EBL is composed of the emission of all stars (25) 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 the star formation rate (26), taking into account the mean dust extinction within galaxies (10, 27, 28). The second approach uses a physical EBL model (29) 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 (14). The two approaches yield consistent results for the SFH, which is well constrained out to a redshift of z ≈ 5—that is, to the epoch 1.5 billion years after the Big Bang (Fig. 3).

Fig. 3 Cosmic star formation history as constrained from optical depth data.

The shaded regions correspond to the 1σ confidence regions on the star formation rate density as a function of redshift, Embedded Image, obtained from two independent methods, based on a physical EBL model (green) and an empirical EBL reconstruction [blue; see (14)]. The data points show the star formation history derived from UV surveys at low z and deep Lyman-break galaxy surveys at high z [see (1) and references therein]. Horizontal and vertical yellow lines represent the redshift bin size and the uncertainty on the cosmic star formation history. Figure S11 includes a more complete set of data from different tracers of the star formation rate.

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

At z ≳ 6, the far-UV background (photon energy > 13.6 eV) is responsible for the reionization of the neutral hydrogen in the Universe, but the nature of ionizing sources has not been conclusively identified. One possibility is that ultrafaint galaxies existing in large numbers can provide the required ionizing photons (31, 32). 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 MAB ≳ –15), with some suggesting a continued steep faint-end slope (33, 34) and others claiming a turnover (35, 36). Our upper limits at z = 5 to 6 on the UV emissivity ρUV < 3.2 (5.3) × 1026 erg s–1 Mpc–3 Hz–1 [1σ (2σ)] (Fig. 4) suggest a turnover of the luminosity function at MAB ~ –14, in agreement with previous claims (35, 36). This still provides abundant photons to drive the reionization.

Fig. 4 Upper limits on the UV luminosity density of galaxies at z ~ 6.

The 1σ and 2σ limits are shown as dashed horizontal lines, light blue and dark blue, respectively. The solid curves show the z ~ 6 UV emissivity from (3336) of the Hubble Frontier Fields (HFF) program as a function of the lower integration limit of the UV luminosity function. The dotted lines correspond to extrapolations beyond the limiting magnitude of the HFF analyses. The data from (35) are those derived using the GLAFIC magnification model. The lines of (34) and (36) have been shifted up by 0.15 dex to account for evolution of their combined sample (z ~ 6 to 7) to z ~ 6. The gray area corresponds to the luminosity required to keep the Universe ionized at z = 6 assuming C/fesc = 30, where C is the clumping factor of ionized hydrogen and fesc is the mean escape fraction of ionizing photons (14).

Supplementary Materials

Authors and Affiliations

Materials and Methods

Figs. S1 to S12

Tables S1 to S5

Data S1

References (43115)

References and Notes

  1. We adopted the following values for the Hubble constant and cosmological parameters: H0 = 70 km s–1 Mpc–1, ΩM = 0.3, and ΩΛ = 0.7.
  2. See supplementary materials.
  3. The contribution of active galactic nuclei is small in comparison (14).
  4. The templates used are gll_iem_v06.fits and iso_P8R2_SOURCE_V6_v06.txt, available at
  5. We used the templates gll_iem_v06.fits and iso_P8R2_TRANSIENT020_V6_v06.txt, available at
Acknowledgments: Funding: Supported by NSF and NASA through grants AST-1715256, NNX16AR72G, and 80NSSC17K0506 (M.Aj., V.P., and A.De.); Icelandic Research Fund grant 173728-051 (K.H.); the Chief of Naval Research and a grant of computer time from the Department of Defense High Performance Computing Modernization Program at the Naval Research Laboratory (J.F.); the Juan de la Cierva program from the Spanish MEC (A.Do.); Wallenberg Academy (J.C.); and Italian Ministry of Education, University and Research (MIUR) contract FIRB-2012-RBFR12PM1F (M.R.). 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 NASA and the U.S. Department of Energy (DOE); the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules (France); the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare (Italy); the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK), and Japan Aerospace Exploration Agency (JAXA) (Japan); and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board (Sweden). Additional support for science analysis during the operations phase came from the Istituto Nazionale di Astrofisica, Italy, and the Centre National d’Études Spatiales, France. This work was performed in part under DOE contract DE-AC02-76SF00515. Author contributions: M.Aj. designed the project and wrote most of the paper; V.P. performed the analysis of the γ-ray data; K.H. derived the constrains on the EBL and the star formation history and wrote all the corresponding text; J.F. derived the results of the stellar population method and wrote the corresponding text; A.De. tested all the EBL models; A.Do. provided all the EBL-related data reported in the figures and wrote the corresponding text; and all co-authors have read, provided comments on, and approved the manuscript. Competing interests: All co-authors declare that there are no competing interests. 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 for analysis are available from the Fermi Science Support Center, The reconstructed optical depth templates, EBL, and SFH are also available at The tool to produce physical models of blazars’ SEDs is available at

Stay Connected to Science

Navigate This Article