THE NuSTAR EXTRAGALACTIC SURVEYS: THE NUMBER COUNTS OF ACTIVE GALACTIC NUCLEI AND THE RESOLVED FRACTION OF THE COSMIC X-RAY BACKGROUND

A complete census of accreting supermassive black holes (SMBH) throughout cosmic time is necessary to quantify the efficiency of accretion, which is believed to drive the majority of SMBH growth (e.g., Soltan 1982; Yu & Tremaine 2002; Di Matteo et al. 2008, Merloni & Heinz 2008). X-ray emission is nearly universal from the luminous active galactic nuclei (AGNs) that signal the most rapid SMBH growth phases, making surveys in the X-ray band particularly efficient at identifying accreting SMBH. Unlike optical and infrared light, X-rays are not diluted by host-galaxy emission, which is generally weak above ∼1 keV. X-rays are also penetrating, and hard (≳10 keV) X-rays are visible through columns up to ${N}_{{\rm{H}}}\sim {10}^{25}\,{\mathrm{cm}}^{-2}$. For even higher columns, AGNs can be identified through scattered X-rays; though, the emission becomes progressively weaker with increasing columns.

Cosmic X-ray surveys with Chandra and XMM-Newton have provided measurements of the demographics of the AGN population and its evolution in the 0.1–10 keV band out to large cosmic distances (see Brandt & Alexander 2015, for a recent review). These surveys are sufficiently complete over a broad enough range of luminosity and redshift that many fundamental questions regarding AGN evolution can be addressed. In the deepest fields, $\gt 80$% of the 2–10 keV cosmic X-ray background (CXB) has been resolved into individual objects (Hickox & Markevitch 2006; Brandt & Alexander 2015).

At X-ray energies above 10 keV, the observational picture is far less complete. Coded-mask instruments, such as INTEGRAL and Swift/BAT, have probed the demographics of hard X-ray emitting AGNs in the very local universe, to redshifts of $z\lesssim 0.1$ (Tueller et al. 2008; Beckmann et al. 2009). The fraction of the CXB resolved by these instruments at its peak intensity (20–30 keV) is ∼1% (Krivonos et al. 2007; Ajello et al. 2012; Vasudevan et al. 2013). Thus, until now, samples of AGNs selected at $\gt 10\,\mathrm{keV}$, which are inherently less biased by obscuration than those at lower energy, could not be used to probe AGN demographics and, in particular, the evolution of highly obscured to Compton-thick (${N}_{{\rm{H}}}\gtrsim {\sigma }_{{\rm{T}}}^{-1}\sim {10}^{24}\,{\mathrm{cm}}^{-2}$ ) sources.

There is, therefore, strong motivation for extending sensitive X-ray surveys to high energy. Simple extrapolations of AGN populations detected by Chandra and XMM-Newton to higher energies based on average spectral properties fail to reproduce the shape and intensity of the CXB at 30 keV (e.g., Gilli et al. 2007). This indicates either that spectral models based on the small samples of high-quality 0.1–100 keV measurements used in CXB synthesis models fail to capture the true spectral complexity of AGNs, or that an additional highly obscured AGN population is present in the redshift range of $0\lesssim z\lesssim 2$. It is likely that both factors are important at some level. A higher fraction of reflected emission than typically assumed, which hardens the emission above 10 keV, appear to be present in moderately obscured, $23\lt \mathrm{log}({N}_{{\rm{H}}}$/cm⁻²) <24, AGNs (Ricci et al. 2011). Even at moderate redshifts, reflection fractions are difficult to constrain with data restricted to $\lt 10\,\mathrm{keV}$ (e.g., Del Moro et al. 2014). In addition, it is difficult at low redshifts to properly measure high obscuring columns, which can lead to large errors in estimating intrinsic AGN luminosities (e.g., Lansbury et al. 2015). It is also challenging to identify Compton-thick objects in the range $0\lesssim z\lesssim 2$ with the limited bandpass of Chandra and XMM-Newton (at $z\gtrsim 2$ these missions sample rest-frame energies above 20 keV). Thus both to characterize highly obscured AGNs and to constrain their evolution at $z\lesssim 2$ requires sensitive surveys at energies above 10 keV.

The Nuclear Spectroscopic Telescope Array (NuSTAR), the first focusing high-energy X-ray (3–79 keV) telescope in orbit (Harrison et al. 2013), has been executing a series of extragalactic surveys as part of its core program, with the aim of measuring the demographics and properties of obscured AGNs. Through contiguous surveys in fields with existing multi-wavelength data combined with dedicated spectroscopic followup of sources serendipitously identified in individual fields, the NuSTAR extragalactic surveys have improved the sensitivity limits in the hard, 8–24 keV band by two orders of magnitude relative to INTEGRAL or Swift/BAT, and have probed a wide range in redshift, out to $z\gtrsim 2$.

In this paper, we present the X-ray number counts (log N–log S) of AGNs, including the first sensitive measurements in the 8–24 keV band. We reach fluxes of S(8–24 keV)$\,\sim \,3\times {10}^{-14}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$, a factor of 100 deeper than previous measurements in this band. We compare the number counts to predictions from X-ray background synthesis models, and to extrapolations from Chandra and XMM-Newton surveys, and compare the intensity of resolved sources to that of the CXB. A companion paper (Aird et al. 2015a) presents direct constraints on the $\gt 10\,\mathrm{keV}$ AGN X-ray luminosity function. We adopt a flat cosmology with ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.7$ and h = 0.7, and quote 68.3% (i.e., 1σ equivalent) errors unless otherwise noted.

For the results in this paper, we include all components of the NuSTAR extragalactic survey program that were completed and analyzed prior to 2015 March. This includes coverage of well-studied contiguous fields as well as identification and followup of sources serendipitously detected in all NuSTAR fields. Figure 1 shows the area as a function of depth for the survey components included in this work. At the shallow end, NuSTAR surveyed 1.7 deg² of the Cosmic Evolution Survey field (COSMOS; Scoville et al. 2007) to a depth of S(8–24 keV) = $1.3\times {10}^{-13}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$. The catalog and results from this survey are presented in Civano et al. (2015; hereafter C15). At the deep end, NuSTAR surveys have reached depths of S(8–24 keV) = $2.5\times {10}^{-14}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$ in the Extended Chandra Deep Field South (ECDFS; Lehmer et al. 2005) over an area of 0.3 deg² and in the deep Chandra region of the Extended Groth Strip (Goulding et al. 2012; Nandra et al. 2015) over an area of 0.23 deg². The ECDFS source catalog is presented in Mullaney et al. (2015; hereafter M15) and the EGS catalog will be presented in J. Aird et al. (2016, in preparation). The serendipitous survey adopts a similar approach used with the Chandra SEXSI survey (Harrison et al. 2003), where fields are searched for point sources not associated with the primary target, and subsequent spectroscopic followup with the Keck, Palomar 200-in, Magellan, and NTT telescopes provides redshifts and AGN classifications. Preliminary results from this survey were published in Alexander et al. (2013), and the 30-month catalog and results will be published in G. Lansbury et al. (2016, in preparation).

Zoom In Zoom Out Reset image size

The overall sample consists of 382 unique sources detected in the full (3–24 keV), soft (3–8 keV), or hard (8–24 keV) bands, of which 124 are detected in the hard band and 226 are detected in the soft band. Figure 2 shows the distribution of rest-frame 10–40 keV luminosity versus redshift for the sources in our sample. Luminosities in this plot are derived from the 8–24 keV count rates (if detected in that band) assuming an unabsorbed X-ray spectrum with a photon index ${\rm{\Gamma }}=1.8$ folded through the NuSTAR response. If the source is not detected in the 8–24 keV band, we calculate luminosities using the 3–24 keV or (if not detected at 3–24 keV) the 3–8 keV count rates. The median redshift for the entire NuSTAR sample is $\langle z\rangle =0.76$, with a median luminosity of $\langle \mathrm{log}({L}_{10-40\mathrm{keV}}$/erg ${{\rm{s}}}^{-1})\rangle =\,44.37$. For comparison, we plot the distribution of AGNs in the Swift/BAT 70-month catalog (Baumgartner et al. 2013). For reference, the dashed line in Figure 2 shows the location of the knee in the luminosity function as a function of redshift from the Chandra-based study of Aird et al. (2015b), extrapolated to the 10–40 keV band. Together, these surveys cover a broad range of luminosity and redshift.

Zoom In Zoom Out Reset image size

2.1. Contiguous Survey Fields

2.2. Serendipitous Survey

To determine the number counts (log N–log S), we adopt a Bayesian approach (see Georgakakis et al. 2008, Lehmer et al. 2012) that assigns a range of possible fluxes to a given source based on the Poisson distribution, which we then fold through the differential number counts distribution. We assume that the differential number counts are described by a single power-law function:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dS}}=K{\left(\displaystyle \frac{S}{{10}^{-13}\mathrm{erg}{{\rm{s}}}^{-1}{\mathrm{cm}}^{-2}}\right)}^{\beta },\end{eqnarray} \tag{ 1 }$$

where K is the normalization at $S={10}^{-13}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$ and β is the slope. Folding the Poisson likelihood for each individual source through the differential number counts given by Equation (1) accounts for the Eddington bias, allowing for the fact that a detection is more likely to be due to a positive fluctuation from a source of lower flux than vice versa. We limit the range of possible source fluxes to a factor of three below the nominal flux limit²⁹ to prevent the probability distribution from diverging at the faintest fluxes as a result of our assumed power-law function.

We optimize the values of the parameters describing the power-law model for the differential number counts by performing an unbinned maximum likelihood fit for all sources detected in a given band (see Georgakakis et al. 2008). We also estimate differential source number counts in a number of fixed-width bins in flux using the ${N}_{\mathrm{obs}}/{N}_{\mathrm{mdl}}$ method of Miyaji et al. (2001), as expanded on in Aird et al. (2010), to account for flux probability distributions. The binned estimate of the differential source number counts is then given by

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{dN}}{{dS}}\right]}_{\mathrm{bin}}={\left[\displaystyle \frac{{dN}}{{dS}}\right]}_{\mathrm{mdl}}\displaystyle \frac{{N}_{\mathrm{obs}}}{{N}_{\mathrm{mdl}}},\end{eqnarray} \tag{ 2 }$$

where ${\left[\tfrac{{dN}}{{dS}}\right]}_{\mathrm{mdl}}$ corresponds to the power-law model for the differential number counts evaluated at the center of the bin, N_mdl is the predicted number of sources in a bin (found by folding the model through the sensitivity curve), and N_obs is the effective observed number of sources, allowing for the distribution of possible fluxes (thus, a single source can make a partial contribution to multiple bins). We estimate errors based on Poisson uncertainties in N_obs as given by Gehrels (1986).

Figure 3 (left) shows the resulting differential source number densities based on the NuSTAR samples from the four survey components combined, in the 3–8 keV band. We also compare the NuSTAR measurements to the best-fit broken power-law functions determined from Chandra surveys by Georgakakis et al. (2008) and a joint analysis of XMM-Newton and Chandra surveys by Mateos et al. (2008). We convert the Chandra (4–7 keV) and XMM-Newton (2–10 keV) number counts to match the 3–8 keV NuSTAR band assuming a ${\rm{\Gamma }}=1.8$ power-law X-ray spectrum. Because the bands are largely overlapping, the choice of photon spectral index does not significantly affect the results. To ease comparison between the different results, we have scaled dN/dS by the Euclidean slope, ${S}^{-2.5}$. The NuSTAR measurements constrain the slope well over the range of $-14\lt \mathrm{log}$ (S/$\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$) < −12.5. The best-fit values for the power-law parameters in this band are log $K=13.84\pm 0.04$, $\beta =-2.81\pm 0.08$.

Zoom In Zoom Out Reset image size

The NuSTAR results generally agree with Chandra and XMM-Newton; though, the slope of the number counts with flux is somewhat steeper. Below $S\sim {10}^{-14}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$, the NuSTAR measurements are poorly constrained, so that the location of the break in the number counts seen by Chandra and XMM-Newton cannot be independently confirmed. The difference in slope of the number counts with flux between NuSTAR and the soft X-ray telescopes is due the different instrument responses and the corresponding uncertainties in converting count rates to fluxes based on simple spectral models. To verify this, we folded the population synthesis model of Aird et al. (2015b) through the NuSTAR response function to predict the observed count-rate distribution. We convert the count rates to fluxes, applying our standard count rate to flux conversion factor. We repeated this exercise adopting the Chandra response function, applying the count rate to flux conversion factor from Georgakakis et al. (2008). For moderately to heavily absorbed sources (${N}_{{\rm{H}}}\sim {10}^{23}$ cm⁻²), the expected count rates are more strongly suppressed using the Chandra response (because the NuSTAR instrument response is more strongly weighted to higher energies). The intrinsic fluxes are therefore underestimated for such sources with Chandra when a single conversion factor is assumed. We have verified that this effect leads to a slope difference in the expected log N–log S at ${f}_{3\mbox{--}8\mathrm{keV}}\sim {10}^{-14}\mbox{--}{10}^{-13}$ erg s⁻¹ cm⁻² that matches the discrepancy seen between the NuSTAR measurements and the Chandra and XMM-Newton measurements shown in Figure 3 (left). Figure 4 (left) shows the integral number counts in the 3–8 keV band.

Zoom In Zoom Out Reset image size

Figure 3 (right) shows the differential source number counts in the 8–24 keV band, along with the best-fit power law, parametrized by log $K=14.3\pm 0.04$, and $\beta \,=-2.76\pm 0.10$. We also plot extrapolations of the Mateos et al. (2008) Chandra+ XMM-Newton counts and the Georgakakis et al. (2008) Chandra counts to the harder, largely non-overlapping 8–24 keV band (dotted and solid lines in Figure 3). The dotted line shows the extrapolation assuming ${\rm{\Gamma }}=1.8$, which is typical of unabsorbed AGNs, and systematically under-predicts the NuSTAR measurements. Using a maximum likelihood analysis, we find that a spectral photon index of ${\rm{\Gamma }}=1.48$ (solid line) provides the best match between the extrapolation of the Georgakakis et al. (2008) model and the NuSTAR data, indicating that a substantial population of absorbed and/or hard-spectrum sources is required to reproduce the NuSTAR measurements.

In combination, NuSTAR and Swift/BAT sample the hard X-ray AGN population over a wide range in flux and redshift. Swift/BAT has measured the number counts at fluxes of $S\gtrsim 3\times {10}^{-12}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$ in the 15–55 keV band for the local ($z\lesssim 0.1$) AGN sample (Ajello et al. 2012). Although there is a gap in the flux range probed by BAT and the $S\lesssim 3\times {10}^{-13}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$ population probed by NuSTAR, extrapolating the BAT differential number counts to fainter fluxes indicates whether there is any evolution between the low-redshift BAT AGN population and the higher-redshift NuSTAR sources.

Figure 5 shows the NuSTAR 8–24 keV number counts together with the BAT measurements from Ajello et al. (2012), which are well described by a single power-law dN/dS with slope $\beta \approx 2.5$ (bold dashed line). We have converted the BAT 15–55 keV fluxes to the 8–24 keV band assuming a photon index of ${\rm{\Gamma }}=1.7$, the value that provides the best fit to the average BAT AGN spectral model from Burlon et al. (2011). The hatched region indicates the measurement uncertainties from Ajello et al. (2012), showing that the formal error is small, about the width of the dashed line. It is clear that the extrapolation of the BAT counts to lower fluxes where the counts are well-constrained by NuSTAR ($f\sim {10}^{-13.5}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$) significantly under-predicts the NuSTAR measurements. This disagreement is not surprising, since there is strong evolution in the AGN population between the local ($z\sim 0.1$) BAT sample and the higher-redshifts probed by NuSTAR.

Zoom In Zoom Out Reset image size

The existence of a population of obscured and Compton-thick objects beyond those directly resolved below 10 keV has long been postulated by CXB synthesis models. These models reproduce the hard spectrum of the CXB using measured AGN X-ray luminosity functions (XLFs), which are well constrained (except in the very local universe) only below 10 keV, together with models for the broadband (∼1–1000 keV) AGN spectra and estimates of the Compton-thick AGN fraction and its evolution (e.g., Gilli et al. 2007, Treister et al. 2009, Ueda et al. 2014). The XLF measurements, assumptions about AGN spectral evolution, and Compton-thick fractions differ significantly among models.

Figure 5 compares the NuSTAR counts to three different CXB synthesis models: the model from Gilli et al. (2007; dashed line with circles), the model of Ueda et al. (2014), and an an updated version of the Ballantyne et al. (2011) model. The details of the Gilli et al. (2007) and Ueda et al. (2014) models are documented in the relevant publications. The updated Ballantyne model differs from Ballantyne et al. (2011) in that it uses the new Ueda et al. (2014) luminosity function, a better spectral model from Ballantyne (2014), the Burlon et al. (2011) ${{\rm{N}}}_{{\rm{H}}}$ distribution, and the redshift evolution and obscured fraction from Ueda et al. (2014). Other aspects, including the normalization of the Compton-thick fraction are the same as in Ballantyne et al. (2011). We include this model because it fits the CXB spectrum even after including the effects of blazers.

All three models are in good agreement with the measured NuSTAR counts. However, all of the models lie significantly above the Swift/BAT measurements at $S\gtrsim {10}^{-11.5}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$. The discrepancy between the BAT measurements and the models cannot be accounted for by uncertainties in the spectral shape; a spectral index of ${\rm{\Gamma }}=2.1$ when converting to the 8–24 keV band is required to make the Swift/BAT 15–55 keV measurements agree with the model predictions. This is significantly softer than the average measured spectral shape, even for unobscured AGNs, and can thus be ruled out. Therefore, this discrepancy appears to be due to evolution in the hard XLF, absorption distribution, or spectral properties of AGNs between the very local Swift/BAT sample and the more distant ($z\sim 0.5\mbox{--}1$) NuSTAR sample that is not fully accounted for in these population synthesis models (see also Aird et al. 2015b).

We have presented measurements of the number counts of AGN with NuSTAR in two bands, from 3–8 keV and 8–24 keV. The data span a broad range in flux, with good constraints covering ${10}^{-14}\,\lesssim$ S(8–24 keV; $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$) $\lesssim \,{10}^{-13}$. The 3–8 keV differential source number densities are in agreement with measurements from Chandra and XMM-Newton; though, the slope measured by NuSTAR is somewhat steeper. The slope difference results from the fact that the NuSTAR effective area curve is weighted to significantly higher energies compared to XMM-Newton and Chandra, which results in the observed discrepancy when using a simple counts to flux conversion factor.

In the 8–24 keV band, we present the first direct measurement of the AGN number counts that includes data above ∼10 keV and reaches down to flux levels $\sim 3\times {10}^{-14}$ $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$. In order to match the NuSTAR number counts, the flux measurements from Chandra and XMM-Newton must be extrapolated to higher energy using a spectrum with photon index ${\rm{\Gamma }}=1.48$. This photon index is significantly harder than the Γ = 1.7–1.9 that characterizes the unobscured AGN population and is also significantly harder than the average spectral index that characterizes the Swift/BAT AGN sample.

The NuSTAR number counts are in good agreement with predictions from population synthesis models that explain the hard spectrum of the CXB using different assumptions about obscuration, the Compton-thick sample, and the spectra shape of AGN in the hard X-ray band (see Figure 5). This directly confirms the existence of a population of AGNs with harder spectra than those typically measured below 10 keV. The spectral hardness could be due either to increased reflection or near Compton-thick absorption. The updated Ballantyne model, for example, uses an AGN spectral model with a strong Compton reflection component with a relative normalization component of R = 1.7. Our measurements of the rest-frame 10–40 keV XLF Aird et al. (2015a) also indicate that a significant population of AGNs with hard X-ray spectra is required to reconcile our NuSTAR data with prior, lower-energy XLF measurements. To what extent this results from obscuration versus higher levels of reflection will be determined by the spectral analysis of NuSTAR sources in the survey fields (A. Del Moro et al. 2016, in preparation, L. Zappacosta et al. 2016, in preparation), and from spectral modeling of high quality data from local AGN samples (M. Baloković et al. 2016, in preparation).

The NuSTAR  8–24 keV number counts lie significantly above a direct extrapolation with flux of the number counts at brighter fluxes, sampled by the Swift/BAT survey in the 15–55 keV band (Figure 5). This discrepancy is not surprising given the known evolution of the AGN population between the low redshift Swift/BAT AGNs and the higher-redshift NuSTAR sample It is interesting that the BAT data are in tension with CXB synthesis models (see Figure 5), which are in good agreement with the NuSTAR measurements. The most natural explanation for the difference is an evolution in the hard XLF, absorption distribution, or spectral properties of AGNs between the very local objects seen by BAT and the more distant ($z\sim 0.5\mbox{--}1$) NuSTAR sample that is not accounted for in the current population synthesis models.

Table 1 provides CXB fluxes measured by hard X-ray instruments in the 20–50 keV and 8–24 keV bands. To convert the intensities to the 8–24 keV band for those instruments that do not cover this energy range (BeppoSAX, INTEGRAL, Swift/BAT), we used Equation (5) of Ajello et al. (2008), which parametrizes the CXB spectrum between 2 keV and 2 MeV based on a fit to available data. The NuSTAR extragalactic surveys have reached depths of S(8–24 keV)$\,\approx \,3\times {10}^{-14}$ $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$. Comparing to the total integrated flux of the CXB as measured by collimated and coded aperture instruments shown in Table 1, this corresponds to a resolved fraction of 33%–40% in the 8–24 keV band (Figure 6), with an additional statistical uncertainty of 5%. Even for the highest measured CXB flux from Churazov et al. (2007) NuSTAR still resolves 33%, which is a significant advance compared to the 1%–2% resolved to-date by coded-mask instruments above 10 keV (Krivonos et al. 2007; Vasudevan et al. 2013).

Zoom In Zoom Out Reset image size

The resolved fraction of the CXB we measure is in good agreement with pre-launch predictions based on CXB synthesis models (Ballantyne et al. 2011). At the current depth, the NuSTAR surveys do not probe the break in the number counts distribution expected, based on extrapolations from Chandra and XMM-Newton, to occur at S(8–24 keV) ∼ 10⁻¹⁴ $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$. Reaching these depths will be challenging because the deep fields are currently background dominated, such that sensitivity improves only as the square root of observing time improves. However, additional exposure in the ECDFS is planned, along with expansion of the deep surveys to cover the CANDELS/UDS field (Grogin et al. 2011; Koekmoer et al. 2011), and continuation of the serendipitous survey, which will better constrain the slope of the number counts distribution above the break and improve spectral constraints on the resolved AGN population.

This work was supported under NASA Contract No. NNG08FD60C, and made use of data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration. We thank the NuSTAR Operations, Software and Calibration teams for support with the execution and analysis of these observations. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA). J.A. acknowledges support from ERC Advanced Grant FEEDBACK at the University of Cambridge and a COFUND Junior Research Fellowship from the Institute of Advanced Study, Durham University.