{'fraction_non_alphanumeric': 0.04323711363355676, 'fraction_numerical': 0.02973042765558512, 'mean_word_length': 4.143385112894843, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Detailed Shape and Evolutionary Behavior of the X-ray Luminosity Function of Active Galactic Nuclei

We construct the rest-frame 2--10 keV intrinsic X-ray luminosity function of Active Galactic Nuclei (AGNs) from a combination of X-ray surveys from the all-sky Swift BAT survey to the Chandra Deep Field-South. We use ~3200 AGNs in our analysis, which covers six orders of magnitude in flux. The inclusion of the XMM and Chandra COSMOS data has allowed us to investigate the detailed behavior of the XLF and evolution. In deriving our XLF, we take into account realistic AGN spectrum templates, absorption corrections, and probability density distributions in photometric redshift. We present an analytical expression for the overall behavior of the XLF in terms of the luminosity-dependent density evolution, smoothed two power-law expressions in 11 redshift shells, three-segment power-law expression of the number density evolution in four luminosity classes, and binned XLF. We observe a sudden flattening of the low luminosity end slope of the XLF slope at z>~0.6. Detailed structures of the AGN downsizing have been also revealed, where the number density curves have two clear breaks at all luminosity classes above log LX>43. The two break structure is suggestive of two-phase AGN evolution, consisting of major merger triggering and secular processes.


Introduction
The AGN/QSO luminosity function and its evolution with cosmic time are key observational quantities for understanding the origin of supermassive black holes (SMBHs) and accretion onto them. This is one of the most important observational product along with the demography of SMBHs, which are made by scaling relations with host galaxy properties (e.g. Magorrian et al. 1998;Häring & Rix 2004;Gültekin et al. 2009; see Kormendy & Ho 2013 for recent review) as well as emerging AGN clustering measurements, which locate typical masses of dark matter halos (DMHs) that the AGNs reside (e.g., Mullis et al. 2004;Yang et al. 2006;Coil et al. 2009;Gilli et al. 2009;Cappelluti et al. 2010;Krumpe et al. 2010Krumpe et al. , 2012 or the halo occupation distribution (e.g., Miyaji et al. 2011;Allevato et al. 2012;Richardson et al. 2013;Chatterjee et al. 2013).
X-ray surveys are practically the most efficient means of finding active galactic nuclei (AGNs) over a wide range of luminosity and redshift. Xray emission from almost all extragalactic point Xray sources with luminosities of L X 10 42 erg s −1 are considered to originate from SMBH accretion, because X-ray emission from the star-formation origin (through supernova remnants and X-ray binaries) of galaxies with the highest star formation rates typically has a luminosity lower than this (Persic et al. 2004;Ranalli et al. 2005;Ptak et al. 2007). Enormous efforts have been made by several groups to follow up the survey X-ray sources with major optical telescopes around the globe, so that we have fairly complete samples of X-ray selected AGNs (Brandt & Hasinger 2005, for review). One of the most important products of such X-ray AGN surveys is the X-ray luminosity function (XLF) and its evolution with cosmic time. Since X-ray emission marks the accretion activity onto the SMBHs, the XLF and its evolution gives a bird's eye view of the moments of SMBH growth process.
The launch of ROSAT and the enormous optical followup efforts of the X-ray sources detected by this satellite enabled us, for the first time, to probe the X-ray emitting AGN populations to highly cosmological distances. From the mid-1990's to mid-2000's, the evolution of the XLF in the soft band (0.5-2 keV) was investigated by numerous authors (e.g., Boyle et al. 1993;Page et al. 1997;Jones et al. 1997;Miyaji et al. 2000aMiyaji et al. , 2001Hasinger et al. 2005). While the soft-band surveys found predominantly unob-scured (type 1) AGNs except at very high redshifts, imaging surveys at 2 < E [keV] 10 with ASCA, XMM-Newton and Chandra enabled us to expanded the investigation to include obscured (type 2) AGNs, predominantly within the Compton-thin limit (log N H [cm −2 ] 24), providing a view of the SMBH accretion process with much better sampling. Analysis of the 2 < E [keV] 10 XLF has been made by various groups using different samples with different approaches and different levels of sophistication (e.g., Ueda et al. 2003;Cowie et al. 2003;La Franca et al. 2005;Silverman et al. 2008;Ebrero et al. 2009;Yencho et al. 2009;Aird et al. 2010;Ueda et al. 2014). Many groups took a simple approach and constructed the 2-10 keV XLF assuming a simple power-law spectrum without absorption corrections. From the combination of AGNs from HEAO-1, ASCA and early Chandra survey samples, Ueda et al. (2003) constructed an absorption-corrected 2-10 keV XLF with a full N H function, as well as a new AGN population synthesis model of the X-ray background. This work was recently revised by Ueda et al. (2014) (herefter, U14), which included both soft and hard samples from diverse sources derived from ROSAT, MAXI, ASCA, XMM-Newton and Chandra surveys. U14 included sophisticated spectral model templates, a revised population synthesis model, discussions on possible contributions from Compton-thick AGNs as well as the growth of SMBHs. A few works have included absorption into analysis (La Franca et al. 2005;Ebrero et al. 2009). Aird et al. (2010) took a simpler approach in spectral assumptions, but used sophisticated statistical techniques and took into account the probability density distribution function (PDF) of photometric redshifts and incompleteness corrections, with Bayesian and Markov Chain Monte-Carlo methods in their model parameter search. In an attempt to make the XLF construction as free from absorption corrections as possible, S. Fotopoulou et al. (2015, submitted to A&A) restricted their analysis to the 5-10 keV band and use the full Bayesian approach in obtaining parameter constraints.
Most authors fitted their XLFs to simple analytical formulae. As a common feature of most such models, the XLF at a given redshift has a (smoothed) two power-law form with a shallower slope at the low luminosity end and a steeper slope at the high luminosity end. This is quite unlike the Schechter function, which is customarily used to describe the optical luminosity function of galaxies and has an exponential decline at higher luminosities. As for the XLF evolution, early evolution models with the pure luminosity evolution scheme (PLE; e.g., Boyle et al. 1993;Page et al. 1997) have been superseded by more complex expressions. Many authors have used the Luminosity-Dependent Density Evolution (LDDE) (e.g., Schmidt & Green 1983;Miyaji et al. 2000a;Ueda et al. 2003;Hasinger et al. 2005) model of one form or another. Yencho et al. (2009) and Aird et al. (2010) considered the models where a smoothed two power-law XLF evolves both in the luminosity and density directions, keeping the XLF shape the same. These are called the (Independent) Luminosity And Density Evolution (called the ILDE or LADE) models. Because the change of at least the low luminosity slope with redshift is observed to be present, generally the LDDE model is preferred, as critically compared by S. Fotopoulou et al. (2015, submitted to A&A).
When AGN number densities are viewed as a function of redshift in different luminosity classes, the AGNs evolve in an "anti-hierarchical" manner, or show "AGN downsizing" (e.g., Ueda et al. 2003;Cowie et al. 2003;Hasinger et al. 2005), i.e., the peak of the AGN number density appears at lower redshifts for low luminosity AGNs than those with high luminosities. This was in contrast to an early theoretical prediction based on AGN triggering by hierarchical merging (Wyithe & Loeb 2003), in which it was predicted that lower luminosity AGNs peak at higher redshifts (or emerge earlier in the cosmic history). In recent years, there have been growing support for the hypothesis that the AGN population might be composed of two populations that have been triggered by different mechanisms. One is a major merger driven triggering, which mainly accounts for high luminosity QSOs and dominant in high redshifts. The other is a secular process, which might include fly-by encounters, minor mergers, and/or disk instabilities, and is dominant in intermediatelow luminosity Draper & Ballantyne (e.g. 2012); Georgakakis et al. (e.g. 2014). A number of theoretical works, involving cosmological N-body and/or N-body+hydrodynamical simulations and semi-analytical treatments in identifying AGNs, have reproduced and/or explained this AGN down-sizing effect (e.g., Di Matteo et al. 2003;Hopkins et al. 2005;Croton 2009;Marulli et al. 2008;Degraf et al. 2010;Fanidakis et al. 2012).
Some models also include "radio mode" (a.k.a. "hot halo" mode). Fanidakis et al. (2012) put emphasis on dividing the AGN activity into "starburst mode", which includes SMBH accretion following both merger-driven and disk instability driven starbursts, and the "hot halo" mode. By also including empirical obscuration effect, they reproduced the AGN down sizing effect. With this general scheme, Fanidakis et al. (2013) explained the tendency that intermediate X-ray luminosity AGNs are associated with more massive dark matter halos (DMHs) than luminous QSOs (Miyaji et al. 2007;Krumpe et al. 2010;Allevato et al. 2011;Krumpe et al. 2012). In order to critically compare these lines of theories and observation, refinements in calculations of the XLF and investigating its detailed evolutionary behavior continues to be an important observational task.
In this work, we use AGNs selected at E > 2 keV from a collection of surveys. Our aim is to construct and investigate the detailed behavior of the XLF. In this particular work, we make our best effort to construct the luminosity function of unabsorbed and absorbed AGNs within the Compton-thin limit, i.e. absorbing column densities of log N H [cm −2 ] 24. Above this column density, Compton scattering causes photons to travel along longer paths within the absorbing medium and are subject to much higher chance of photo-electric absorption. In such Comptonthick (CTK) AGNs, X-ray emission at E 10 keV, which is observable with ASCA, XMM-Newton, and Chandra, is highly suppressed unless they are at very high redshifts. Thus AGNs detected in the E 10 keV surveys are considered to be dominated by Compton-thin (CTN) AGNs, with the exception of a small number of X-ray sources, for which X-ray emission is dominated by in scattered and reflected components (Brunner et al. 2008;Brightman & Ueda 2012;Brightman et al. 2014). In the local universe, higher energy X-ray surveys (20 E 200 keV) such as those available from Swift BAT (e.g., Tueller et al. 2010;Ajello et al. 2012), INTE-GRAL (e.g., Beckmann et al. 2009;Treister et al. 2009) can detect modestly Compton-thick (CTK) AGNs 24 log N H [cm −2 ] 25 and quantify their space density and emissivities. NuSTAR (e.g., Alexander et al. 2013) is extending this to higher redshifts (z ∼ 0.5 − 1).
We present here a new estimate of the X-ray luminosity function evaluated at the traditional 2-10 keV rest frame. Complementary to a recent work by Ueda et al. (2014), whose emphasis is on the X-ray population synthesis of the X-ray background as well as discussions on the role of the CTK AGNs, our emphasis is on the detailed behavior of the shape and evolution of the XLF as well as discussions of possible systematic errors. In particular, in order to assess the uncertainties due to photometric redshifts, we include the PDFs of photometric redshifts (photo-z's) into analysis and compare with the results obtained by using single best-fit photo-z's. We focus on the X-ray luminosity function of unabsorbed and Comptonthin (CTN) AGNs because of the wide availability of E 10 [keV] AGNs. Also X-ray luminosity L X , which is the independent variable of the XLF presented in this paper, is primarily defined as the intrinsic 2-10 keV luminosity, i.e., before absorption and not including the reflection component.
In this work, the sample is greatly expanded by the addition of the datasets from XMM-COSMOS ) and C-COSMOS (Elvis et al. 2009), which are parts of the COSMic evolution Survey (COSMOS) ). In Sect. 2, we explain our COSMOS sample, including the construction of our combination of the XMM-COSMOS and C-COSMOS datasets. The details of the selection criteria for other samples from the literature are explained in Sect. 3. At the end of this section, we explain our incompleteness correction and derive the N (> S) function for extragalactic X-ray sources in each sample.
In sect. 4, we explain the computation of our estimated XLFs. Model parameter estimations using a Maximum-likelihood fitting method, which also include absorption and the probability density distributions (pz-PDFs) of photometric redshifts, are explained. This section also covers the N obs /N mdl estimation of the binned XLF in the presence of the absorption and pz-PDFs.
Sect. 5 presents our global expression of the Luminosity-Dependent Density Evolution (LDDE). In Sect. 6, we take a closer look at detailed behavior of the XLF in redshift divided shells and the number density/emissivity evolution in divided in luminosity classes. We also present binned XLFs. The results are discussed in Sect. 7 and a concluding summary is given in Sect. 8. Throughout this work we use a Hubble constant H 0 = 70h 70 km s −1 Mpc −1 Ω m = 0.3, and Ω Λ = 0.7. The h 70 dependence is explicitly stated. If units are omitted, X-ray luminosities (L with any sub/super scripts) are measured in units of h −2 70 erg s −1 , and column densities N H in cm −2 . The symbol log signifies a base-10 logarithm and ln a natural logarithm.

The XMM-COSMOS Sample
As a part of the COSMOS collaboration, the XMM-COSMOS program (PI=G. Hasinger Hasinger et al. 2007) included observations the entire 2.14 deg 2 COSMOS field with the EPIC camera on board the XMM-Newton observatory over three observing cycles, AO-3, 4 and 6. X-ray point source catalogs have been produced using a sophisticated source detection and characterization procedure (Cappelluti et al. , 2009. The procedure was developed as a part of the XMM-SAS package 1 . In the combined dataset, which is composed of 55 XMM-Newton observations with a total exposure of ∼1.5 Ms, 1887 unique X-ray sources have been detected. Optical counterparts and their multiwavelength properties were published for the X-ray source detected in the first 12 XMM-Newton observations by Brusa et al. (2007) and all 53 successful observations made before 2007 by Brusa et al. (2010) respectively.
In this work, we use X-ray sources detected in the hard (EPIC energy channel between 2 and 8 keV) band in the 53 XMM-COSMOS observation data with a maximum-likelihood of M L > 10, which corresponds to the false detection probability of e −10 ≈ 4.5 × 10 −5 . Brusa et al. (2010) published a multiwavelength catalog for the 53-field XMM-COSMOS sources with a primary 2-10 keV flux cut of S X ≥ 3 × 10 −15 erg s −1 cm −2 , which is our XMM-COSMOS base catalog. After this flux cut, we are left with 923 hard-band selected XMM-COSMOS sources. We use the EPIC 2-8 keV countrate to 2-10 keV flux energy conversion factor (hereafter referred to as ECF) obtained by assuming a power-law spectrum with photon index Γ = 1.7 absorbed by a Galactic absorbing column density of N H = 2.5 × 10 20 cm −2 (Cappelluti et al. 2009).

The Chandra-COSMOS Sample
With its unprecedented spatial resolution, the Chandra X-ray observatory (CXO) provides superb point source sensitivity in deep surveys. The Chandra COSMOS program (C-COSMOS, PI=M. Elvis;Elvis et al. 2009) is a large Chandra program with a total exposure of ∼ 1.8 Ms with 49 observations with the ACIS-I instrument. The observed fields are arranged in 7×7 overlapping tiles, covering a central 0.9 deg 2 of the COSMOS field. Inclusion of the information obtained by the recent Chandra Legacy COSMOS survey (2.5 Ms, PI=F. Civano) will be a topic of a future work. A total of 1781 unique X-ray sources have been detected with a likelihood threshold of M L > 10.8 in any of the full (0.5-7keV), soft (0.5-2 keV) or hard (2-7 keV) bands, out of which 1017 are hard band detected (Elvis et al. 2009;Puccetti et al. 2009).
In Elvis et al. (2009), the 2-10 keV flux is given for an ECF calculated assuming a Γ = 1.4 power-law. To treat the sample in the same way as XMM-COSMOS, we have converted the 2-10 keV flux using the ECF assuming Γ = 1.7 for our further analysis. Unlike the case of XMM-COSMOS, Civano et al. (2012) did not impose any further flux limit in their multiwavelength identification catalog. However, we have imposed a 2-10 keV flux (Γ = 1.7) limit of S X > 1.6 × 10 −15 erg s −1 cm −2 , below which the log N − log S curve spuriously rises due to the Eddington bias. This removes the 24 faintest hard-band detected sources. Note that the log N −log S curve in Fig. 9 of Elvis et al. (2009) effectively plots down to this flux level. This leaves 993 C-COSMOS hard-band sources.

Sensitivity map combination
In order to take advantage of both the XMM and C-COSMOS surveys, we have made a combined sample, which we call the XC-COSMOS sample, by selecting X-ray sources from the more sensitive survey out of the two at any given source position. For this purpose, we have compared the hard X-ray sensitivity map for the 53-field XMM-COSMOS survey (Cappelluti et al. 2009) with that for C-COSMOS (Puccetti et al. 2009). The XMM-COSMOS sensitivity map has been truncated at the lowest limiting flux of the optical catalog described above. In order to compare both maps on the same grounds, we scaled the C-COSMOS sensitivity map using the ECF for the photon index of Γ = 1.7. We further smoothed the C-COSMOS sensitivity map with a Gaussian filter with σ = 16 ′′ and combined with the XMM-COSMOS sensitivity map by choosing the smaller limiting flux at each point. We have also generated a mask file that indicates whether the XMM-COSMOS or C-COSMOS source should used at each position. In effect, the C-COSMOS sources should be used in most of the region covered by C-COSMOS, except for small areas near corners and edges. The combined sensitivity map is shown in Fig. 1. In the combined sample, there are 1400 extragalactic point sources (mostly AGNs), 18 Stars, and 29 unidentified sources. Overall the fraction of sources identified with AGNs/galaxies with spectroscopic or photometric redshifts and stars is 98%. Note that the XC-COSMOS AGNs comprises as many as ∼ 43% of all AGNs in our global AGN sample defined below.
One of the greatest advantage of the COSMOS survey is the availability of multi-epoch, homogenized, deep 31 band photometry, covering from the UV to the MIR wavelengths, including 12 intermediate bands from Subaru (Taniguchi et al. 2007). These data, along with the correction for variability, and the use of hybrid AGN templates, have allowed the COSMOS survey to reach a mean photometric redshift accuracy of σ ∆z/(1+zspec) ≈ 0.015 with a low fraction of outliers (∼5%) for both the XMM and Chandra COSMOS samples (Salvato et al. 2009(Salvato et al. , 2011. In this work, the newest photometric redshifts by Salvato et al. (2011) have been used when there is no available spectroscopic redshift. The associated probability density distributions of photometric redshifts are also used. In addition to the spectroscopic redshifts included in Brusa et al. (2010), newly available proprietary spectroscopic redshifts from spectroscopic survey programs within the COSMOS consortium have been included.

The Swift BAT Survey Sample
Based on the Swift BAT slew scan data accumulated over 60 months, Ajello et al. (2012) produced a catalog of AGNs detected in the 15-55 keV band, providing the current most complete census of nearby bright AGNs over the entire sky. These include unabsorbed and absorbed AGNs up to moderate CTK AGNs. For CTN AGNs, photoelectric absorption affects very little to the 15-55 keV flux and thus there is no need for complicated absorption corrections. However, since we are evaluating intrinsic XLF at 2-10 keV, we need spectral template for unabsorbed AGNs to convert from 15-55 keV to 2-10 keV. In our simplest treatment, we assume a power-law spectrum with a photon index of Γ = 1.7, which is approximately (within ∼ 4%) equal to the corresponding ratio for the unabsorbed AGN template spectrum (see Sect. 4.3). Under this assumption, we can convert the 15-55 keV flux to 2-10 keV flux by multiplying by a factor of 0.705. In the limiting flux-survey area curve from Fig. 1 of Ajello et al. (2012), the flux is converted to the 2-10 keV value in this way. In our full treatment in Sect. 4.3, we use a luminosity-dependent effective photon index for the conversion.

ASCA LSS/MSS
The Advanced Satellite for Cosmology and Astrophysics (ASCA) was the first X-ray observatory with focusing optics capable of making spectrocopic imaging in the 2-10 keV band. The combination of ASCA Large Sky Survey (ALSS) and the ASCA Medium Sensitivity survey (AMSS) fills the gap between all-sky surveys and deeper surveys with XMM-Newton and Chandra, occupying a unique regime in the z − L X space at E > 2 keV.
The ALSS covers a continuous area of 5.5 deg 2 near the North Galactic Pole with a flux limit of S X ∼ 1 × 10 −13 [erg s −1 cm −2 ] (2 -10 keV) (Ueda et al. 1998(Ueda et al. , 1999. Thirty AGNs detected with the SIS instrument are completely identified by (Akiyama et al. 2000). The AMSS is based on a serendipitous X-ray survey with the GIS instrument (Ueda et al. 2001(Ueda et al. , 2005 and has a survey area of ∼ 85 deg 2 . The flux limit ranges from 5 × 10 −12 to 3 × 10 −13 erg s −1 cm −2 . The identification catalogs in the northern part of the AMSS (Akiyama et al. 2003, AMSSn;) and the southern part (AMSSs; Akiyama et al. in prep.) contain 74 and 20 spectroscopically identified non-blazer AGNs respectively. Three X-ray sources are left unidentified and the identification completeness is 97% for the combined ALSS and AMSS sample.
The 2-10 keV fluxes of these sources have been obtained from the ASCA GIS count rate assuming a Γ = 1.7 power-law using the GIS response matrix.

The XMM-Newton Hard Bright Serendipitous Sample (HBSS)
The Hard Bright Serendipitous Sample (HBSS Della Ceca et al. 2004) is a subsample detected in the 4.5 -7.5 keV band from the larger survey of the XMM-Newton Bright Survey (XBS). Della Ceca et al. (2008) defined a complete fluxlimited sample of 67 sources with MOS2 count rates larger than 0.002 [cts s −1 ] (4.5 -7.5 keV) over 25 deg 2 of the sky. The MOS2 count rates in the HBSS were converted to observed fluxes in the 2-10 keV band, assuming a photon power-law index of 1.7 (using an ECF of 1.044 × 10 −11 erg cm −2 cts −1 . The corresponding flux limit is 2.1 × 10 −13 erg s −1 cm −2 . These sources are completely identified spectroscopically except for two, making the identification completeness of ∼ 97%. The 65 identified sources consist of 62 extragalactic point sources (mostly AGNs), 2 stars and 1 cluster of galaxies.

The XMM-Newton Medium-sensitivity Survey (XMS)
The XMM-Newton Medium-sensitivity Survey (XMS; Barcons et al. 2007) is a serendipitous Xray survey and optical identification program of sources with intermediate X-ray fluxes discovered in 25 XMM-Newton high Galactic latitude fields covering a sky area of ∼3 deg 2 . For this analysis the XMS-H sample, selected in the 2-10 keV band, is used. Barcons et al. (2007) assumed a spectral index of Γ = 1.7 for the count-rate to flux conversion for the hard band, and therefore we use the 2-10 keV flux as shown in their dataset. In our analysis, we use the the sub-sample of XMS-H defined by Hasinger (2008).
The original Barcons et al. (2007) sample contains 159 sources with 2-10 keV fluxes > 3.3 × 10 −14 [erg cm −2 s −1 ] with a spectroscopic identification fraction of 83% (27 unidentified sources). However, the actual spectroscopic completeness limit varies from field to field. Therefore Hasinger (2008) selected a subsample that is comprised of further flux cuts on a field-by-field basis. This way a cleaner XMS sample could be defined, including 128 extragalactic point X-ray sources and almost half the number of unidentified sources, i.e. achieving an identification fraction of 91%. The survey solid angle has been corrected accordingly. See Sect. 2.4 of Hasinger (2008) for further details.

Hellas2XMM
HELLAS2XMM (Baldi et al. 2002) is a serendipitous survey based on suitable XMM-Newton pointings (complementary to the XMS). The assumed spectral index for the ECF for the 2-10 keV flux in Baldi et al. (2002) was Γ = 1.7, and therefore we use their flux as provided.
As in the case of the XMS, we use the subsample compiled by Hasinger (2008). Fiore et al. (2003) presented optical identifications and spectroscopic redshifts for 122 sources selected in the 2-10 keV band in five XMM-Newton fields, covering a survey area of 0.9 deg 2 . Later Maiolino et al. (2006) presented additional redshifts from VLT NIR spectroscopy of optically extremely faint objects for two additional sources. Both Fiore et al. (2003) and later Mignoli et al. (2004) tried to estimate redshifts for the remaining, optically faint, unidentified sources. Cocchia et al. (2007) published photometry and spectroscopic redshifts in five additional HELLAS2XMM fields providing 59 new redshift identifications for the sample of 110 new sources with an additional solid angle of 0.48deg 2 . As explained in Hasinger (2008), further flux cuts were made on the field-by-filed basis. In addition, different flux cuts were applied to inner and outer regions of each field. The redshift completeness in this subsample of HELLAS2XMM is thus 84%. In this paper, we further raise the flux limit from 1.0×10 −14 (adopted by Hasinger 2008) to 1.63×10 −14 erg s −1 cm −2 to make the identification 90% complete, leaving 115 AGNs.

Lockman Hole (XMM)
Because of the extremely low Galactic column density (N H = 5.7 × 10 19 cm −2 ), Lockman hole has been selected as a target of X-ray deep survey observations, including a XMM-Newton survey Brunner et al. 2008). The base X-ray source catalog for this work is from Brunner et al. (2008) from 18 individual pointings, twelve of which are centered near the same pointing direction and the remaining pointings are spread out over about ∼ 30 ′ in right ascension. Brunner et al. (2008) concentrated on the deep central circular area with a radius of 15 ′ (= 0.20deg 2 ). Their catalog contains 409 X-ray sources, out of which 266 are detected in the hard (2-4.5 keV) band at the maximum-likelihood (ML) value of 6.0 or above. Based on simulations, they showed that the sources above this ML threshold, in combination with the flux-sensitivity curve in their Fig. 5, reproduces the log N − log S accurately.
Their listed 2-10 flux is converted from the EPIC 2-4.5 keV count rate assuming a Γ = 2.0 power-law spectrum. For our purposes, we have converted their Γ = 2.0 based 2-10 keV flux to the Γ = 1.7 based flux using a combined response matrix of the XMM-Newton EPIC PN and two MOS detectors. This corresponds to a multiplication of the Γ = 2.0 based flux by 1.16. The limiting-flux -area curve in Fig. 5 of Brunner et al. (2008) has also been converted to the Γ = 1.7-based flux.
Of the 266 hard X-ray selected sources, 60% have published or unpublished spectroscopic redshifts (including 5 stars). Including photometric redshifts by Fotopoulou et al. (2012), 98% are identified with either spectroscopic or photometric redshifts. Only 4 are left unidentified. This results in 257 hard-band selected extragalactic X-ray point sources. Probability density distributions of the photometric redshifts are also available, which we use in Sect. 4.2.

CLASXS and CLANS
The catalogs from two intermediate-depth Chandra surveys near the Lockman Hole, the Chandra Large Area Synoptic X-ray Survey (CLASXS; Yang et al 2004) and the Chandra Lockman Area North Survey (CLANS; Trouille et al. 2008) are also included in our analysis. The Xray source identifications, and limiting-flux -survey area curves of both surveys (as well as those of the CDF-N; see below) are conveniently included in Trouille et al. (2008) as a part of their "Opt-X" project. Additional spectroscopic redshifts of CLANS sources can also be found in Trouille et al. (2009). Their tables contain 2-8 keV fluxes based on the 2-8 Chandra ACIS-I count rates with the ECFs calculated with variable Γ and ACIS-I response matrices based on old calibrations. Therefore we have converted from their 2-8 keV count rate to 2-10 keV flux assuming Γ = 1.7 and with ancillary response files (ARFs) we have re-created for the ACIS-I aim point using the Chandra CALDB 4.5.7 (the latest calibration as of our analysis). For each of the CLASXS and CLANS surveys, the applied response is the exposure weighted average of the responses at the aim points of the observations belonging to the survey. The flux in their limiting flux-survey area curves ( Fig. 5 in Trouille et al. 2008) have been adjusted accordingly. The curves for the probability of detection of 30% are used, because they reproduce the log N − log S relation in the simulations by Yang et al (2004). Because their curves only apply to the sources with offaxis angles of < 8 ′ , we use sources that meet this criterion. Also, we only include X-ray sources that are detected at a signalto-noise ratio of > 3 in the hard band in order to be compatible with their flux-survey area curve. Moreover, in order to be compatible with there flux-limit -area curve, they also imposed offaxisdependent flux limits to exclude some sources for their log N − log S calculations. These imposed limits cannot be reproduced with the published information. However, this effect does not significantly affect our analysis as described below.
In order to minimize the effects of the "completeness correction" (Sect. 3.12), we impose at least 90% identification completeness for each survey. In order to achieve this, we have excluded the sources with 2-10 keV flux (Γ = 1.7 based) of S X < 6.5×10 −15 and S < 6.5×10 −15 erg s −1 cm −2 , leav-ing 250 and 119 AGNs for CLASXS and CLANS respectively. At these flux limits, the survey area is already ∼ 95% of the total geometric area with offaxis angle < 8 ′ . Thus the effects of offaxisdependent flux limits are minimal.

Chandra Deep Field North (CDF-N)
Our base X-ray source catalog for the 2Ms CDF-N is from Alexander et al. (2003). As in the cases of CLASXS and CLANS, we have used the hard-band detected sources that have signalto-noise ratio of 3 or better to match with the flux-limit survey area curve shown in their Fig.  19. In order to convert their 2-8 keV count rate to 2-10 keV flux, we have re-generated the ancillary response file for an exposure-weighted average of all the CDF-N Chandra observations used by Alexander et al. (2003), assuming a Γ = 1.7 power-law and using CALDB 4.5.7. Since their flux-limit survey area curve is based on Γ = 1.4 and old calibration, we have adjusted their curve for our flux definition.
We use the optical identifications listed in Trouille et al. (2008), which include redshifts from the literature as well as from their own program. We have further imposed a lower flux cut of S X > 8.9 × 10 −16 erg s −1 cm −2 to ensure an identification completeness of at least 90%, leaving 182 extragalactic X-ray point sources.

Chandra Deep Field South (CDF-S)
An extensive X-ray source and identification catalog of the 4 Ms Chandra Deep Field South was published by Xue et al. (2011), which also includes a list of previously published spectroscopic and photometric redshifts. Subsequently Lehmer et al. (2012) generated limiting fluxsolid angle curves and the log N − log S relation. This work took into account flux probability distributions and recovery functions for those sources with a threshold probability of false detection of P thresh = 0.004, as adopted by Xue et al. (2011).
We have chosen hard X-ray sources that have been detected in the 2-8 keV band from their catalog. Their count-rate to flux conversion is based on the effective photon index Γ from the ratio of the hard and soft band count rates, and therefore the ECF is variable. For each source, we have converted their 2-8 keV flux (based on their estimated Γ) to 2-10 keV flux based on Γ = 1.7 using the aim point response matrix of one of the CDF-S observations (OBSID=12049). As shown below, the value of the count-rate to flux conversion factor is subject to the variation of instrumental response over time as well as calibration updates, but the relative change of conversion factors among different spectra is insensitive to ACIS-I response matrix used. The 2-10 keV flux (after the conversion) of the faintest object amongst the hard-band detected sources was 8.3×10 −17 erg s −1 cm −2 .
We have used the 2-8 keV band limiting flux -solid angle curve in Fig. 2 of Lehmer et al. (2012). Their curves are drawn for Γ = 1.0 and we have converted their relation to Γ = 1.7 using the method described above.
We have further included new spectroscopic redshifts that are not included in the Xue et al. (2011) catalog. In the catalog of Xue et al. (2011), there are two objects with photometric redshifts above z ph > 7. The second solution of these objects by Luo et al. (2010) are both at z ph ≈ 3. However, a new spectrocopic redshift obtained for one of them (XID=28) is z sp = 0.69. Thus it is possible that the other also suffer a similar catastrophic photometric redshift error also. Thus we have excluded the remaining z ph > 7 object from our sample. The overall identification completeness is 97% (including photometric redshifts) and the spectrocopy completeness is 71%.

Removing Compton-thick AGNs
The 2-10 keV selected samples, especially at low redshifts, are highly selected against CTK AGNs. Even the Swift BAT sample selected at 15-55 keV suffers from suppression of the CTK population, because of the Compton scattering itself which results in subsequent lowering of the photon energy as well as longer paths in absorbing medium. However, there always is some spill over of the CTK population to the sample defined here. There are two approaches in treating the situation. One is to include the CTK population into the model and make model fits to the sample that may include CTK AGNs. In other words, the amount of CTK AGN spill-over is estimated by the model. The model of U14 introduced the parameter f CTK , which is the ratio of the number density of the CTK AGNs (24 ≤ log N H < 26) to that of absorbed CTN AGNs (22 < log N H < 24), and they then assumed the same f CTK at all luminosities and redshifts. They showed that this is at least consistent with various observational estimates of the CTK AGN population in the literature. However, U14 also estimated that the observational constraint on this parameter was 0.5 < f CTK < 1.6, carrying a factor of three uncertainty. This is consistent with recent results by Buchner et al. (2015). An alternative approach is to remove known CTK AGNs from the sample on a best effort basis and construct an XLF model that represents only the CTN AGNs. We take the latter approach.
There are a number of studies aimed at identifying CTK AGNs among the AGN samples used in this paper. Because the Swift BAT sample is selected in the 15-55 keV band, it has some sensitivity to CTK AGNs with log N H 25. Ajello et al. (2012) list 18 known CTK AGNs from their preceding work of Burlon et al. (2011) and other literature. These have been excluded from our analysis. Ajello et al. (2012) estimated that a few unknown CTK AGNs might still be present in this sample, the effect of which on our analysis is negligible. In their XMM-Newton Lockman hole sample, 11 AGNs meet the X-ray color-color criteria of CTK AGNs defined by Brunner et al. (2008). None of them is in our hard-band selected sample above our luminosity cut. Lanzuisi et al. (2015) selected 10 CTK AGNs based on spectral fits for the XMM-COSMOS sources that have 30 counts or more in the full 0.3-10 keV band. None of them are in the XMM part of the XC-COSMOS AGN sample. Note that all but one XMM-Newton AGNs in the XMM part of the XC-COSMOS sample have 0.5-10 keV counts greater than 30, thus, Lanzuisi et al. (2015) would have identified almost all CTK AGNs in it if any were present. One AGN in our CLANS sample is among the five X-ray selected CTK AGNs by Polletta et al. (2007) and this has been excluded for our analysis. Nine of the CDF-N AGNs are among the CTK AGN candidates derived from the X-ray spectral analysis by Georgantopoulos et al. (2009). These are also excluded from our analysis. Georgantopoulos et al. (2009) performed spectral analysis of the CDF-N AGNs with a 2-10 keV flux of > 1 × 10 −15 [erg s −1 cm −2 ], which is approximately the same as our flux cut to the CDF-N sources to meet our completeness criteria. Recently Brightman et al.
(2014) published a table of extensive X-ray spectral analysis results of the X-ray sources in the CDF-S, AEGIS-XD and C-COSMOS surveys with an intention to identify CTK AGNs. Because of the low background level of Chandra surveys, they attempted spectral analysis to the sources down to 10 Chandra ACIS counts. Following the criteria by Brightman et al. (2014), we consider as probable CTK AGNs those that satisfy the both of the conditions; the best-fit value of log N H is greater than > 24 and the lower bound of its 90% confidence error range is > 23.5. We exclude these from from our analysis also. The number of excluded sources by the CTK AGN criteria are 24 for the C-COSMOS and 23 for the CDF-S samples respectively. The full band counts of almost all faintest sources above our flux cuts (see above) well exceed 10 counts for both C-COSMOS and CDF-S samples. Thus CTK AGNs that may be still remain in our sample can be neglected. The CTK AGNs are still a minor population in our sample, even at the faint end. The CTK AGNs identified above comprise 19 ± 6%, 9 ± 3%, 3.6 ± 0.5%, and 1.9 ± 0.3% for the log S x ranges of (−16.

Sample Summary
In each sample except for Swift BAT, we have converted the 2-10 keV fluxes of all the sources using the energy conversion factor (ECF) assuming a power-law spectrum with a photon index of Γ = 1.7. To represent fluxes calculated in this way, we use the symbol S X,obs , or S Γ=1.7 X,obs when we want to emphasize that we use the ECF based on Γ = 1.7. We also define the observed luminosity where d L (z) is the luminosity distance at the redshift z. Implementation of the K-corrections and absorption corrections using realistic AGN spectra is discussed in Sect. 4.3. In our further analysis, we are interested in X-ray sources with intrinsic luminosities of log L X > 42.0 and we use, in effect, X-ray point sources with log L X,obs 41.5. Thus our analysis is strongly de-selecting non-AGN Xray sources, the number density of which is lower than AGNs by at least an order of magnitude at log L X 41.5 (Persic et al. 2004;Ranalli et al. 2005;Ptak et al. 2007). Our sample is plotted in the z − log L X,obs space in Fig. 2.
The properties of each sample used in our analysis is summarized in Table 1, where the name of the survey, the survey solid angle, the limiting flux corresponding to 20% of the maximum solid angle, the number of extragalactic point sources with spectroscopic redshift, those with photometric redshift, and identification completeness are listed for each survey. The numbers in this table do not include CTK AGNs that have been discussed in the previous section nor those below our luminosity cut.

Incompleteness Correction and Number counts
While the identification completeness of our selected samples are at least 90%, we further account for unidentified sources as follows. For samples with five or less unidentified sources, we have multiplied the survey area in the flux-limit -survey area curve by the identification completeness of the sample over the whole flux range, assuming that the redshift (or the source type) probability distribution of an unidentified source is the same as the redshift/source type distribution of the identified sources in the same sample.
If there are more than five unidentified sources, we have made a similar adjustment of the survey area in flux blocks. We divide the X-ray sources into flux blocks from the brightest to faintest, in such a way that a new block starts after the fifth unidentified source. The last (faintest) block may contain one to five unidentified sources. For each block, we have adjusted the survey area by multiplying by the completeness level of the block.
The assumption is that the redshift/source type probability distribution is the same as that of the identified source at similar X-ray fluxes.
The resulting incompleteness-corrected limitingflux curve of each survey as well as that for the total of all surveys are plotted in Fig. 3(a). In order to verify the integrity and consistency of our samples and the level of systematic errors and/or the effects of cosmic variance, we plot the "tilted" cu-mulative number counts S 1.5 X14,obs N (> S X ) curves in Fig. 3(b), where S X14,obs is the 2-10 keV observed flux (see above) measured in units of 10 −14 erg s −1 cm −2 , for each of our identified extragalactic X-ray point source sample. The cumulative count is calculated by: (2) where Ω c (S X,obs,i ) is the corrected survey area at the limiting flux of the i-th source S X,obs,i . One sigma errors, calculated by We follow the same procedure as our previous work (Miyaji et al. 2000a(Miyaji et al. , 2001Ueda et al. 2003Ueda et al. , 2014 and determine the best-fit parameters of parameterized model of the XLF.
As our maximum-likelihood estimator, we use either one of the following two forms: where i goes through each AGN in the sample and N (log L X , z) (N-function) is the expected number density of AGNs in the sample per logarithmic luminosity per redshift, calculated from a parameterized analytic model of the XLF: where d A (z) is the angular diameter distance, χ is the radial comoving distance (e.g., Sect 4.3 of  The completeness-corrected survey area curve as a function of the 2-10 keV flux limit. A photon index of Γ = 1.7 is assumed for all count-rate to flux conversions. The total area as well as those for individual surveys are shown as labeled. (b) The "tilted" cumulative number counts S 1.5 X14,obs N (> S X,obs ) for individual samples with 1σ errors. The multiplication factor S 1.5 X14,obs is applied to the cumulative number to make the Euclidean slope flat. The meaning of line colors and styles corresponds to those in the the flux limit -survey area curves in panel (a). For visibility, the data points corresponding to the 4 brightest objects in each sample are not displayed. (Color versions of these figures are available in the electronic version.) Schneider 2006), and A(S X,obs ) is the survey area as a function of limiting observed X-ray flux S X,obs as shown in Fig. 3(a). For a power-law spectrum with a photon index of Γ, d L (z) is the luminosity distance.
The minimization of L with respect to model parameters gives the best-fit model. With Eq. 5, employed in our previous works, the parameter that represents the global normalization cannot be a fitting parameter, since within the estimator itself, the model number density is normalized. On the other hand, with Eq. 4, employed by Marshall et al. (e.g., 1983), one can treat the normalization as a fitting parameter.
We use the MINUIT package (James 1994) distributed as a part of the CERN program library for the minimization procedure and parameter error search. The MINUIT command "MINOS" gives errors taking parameter correlations into account. In cases of global expressions (Sect. 5), "MINOS" fails to give errors for many parameters if we use Eq. 4. In this case, we use Eq. 5 and the normalization and its errors are calculated independent of the fitting process as follows. The model normalization A (as a generic symbol of various normalizations that appears in the following subsections) is determined, such that the total number of expected objects is equal to the number of AGNs in the sample (N obs ). The estimated 1σ confidence error for A is taken to be A/ √ N obs and does not include correlations of errors with other parameters. In our results, we indicate which method has been used for each fit.

Photometric redshift full probability density distribution
For the CLASXS, CLANS, COSMOS, LH, and Chandra Deep fields North and South datasets, significant fractions of redshifts rely on photometric redshift determinations. Especially for fainter optical sources, which are systematically selected against for spectroscopic measurements, the error in the primary peak and catastrophic failure rate increase rapidly. Thus we include the probability density distribution (PDF) in redshift space for photometric redshifts, whenever available.
Let the PDF of the i-th object be p i (z) (normal-ized with p(z)dz = 1). The X-ray luminosity of the i − th object L x,i (z) is now a function of z. The PDFs can be incorporated into Eqs. 4 & 5 by replacing N (log L xi , z i ) by: We use the PDFs derived by Salvato et al. (2011) andFotopoulou et al. (2012) for the COS-MOS and LH samples respectively. They give the PDFs in bins of ∆z = 0.01 in the 0 < z < 7 range. To reduce the computational time and at the same time not to sacrifice the accuracy, we rebinned the PDFs into ∆z/(1 + z) ≈ 0.03 bins and then neglected the new bins that have a probability of less than 2%. For the CLASXS, CLANS, and CDF samples, no PDF information was available in the public domain at the time of our analysis. However, for the CDF-S, Xue et al. (2011) show the first and second peaks of the photometric redshift PDF derived by Luo et al. (2010) for some of their photometric redshifts. For those, we assumed that the PDF is the sum of two delta functions centered at these peaks, each of which has a weight of 0.5. The effects of considering the PDFs are discussed in Sect. 7.1.

Using Realistic AGN Spectra
Our sample contains AGNs at various absorption levels. Thus we have considered realistic absorbed spectra in our XLF calculations, following the approach by U14 (see also Miyaji et al. 2000b;Ueda et al. 2003;Ebrero et al. 2009). In this work, we do not estimate the N H on individual AGNs nor derive the N H function from our own dataset. Instead, we include in our model fitting the latest luminosity and redshift dependent N H function, f (L X , z; log N H ), from the literature and AGN template spectrum. Ueda et al. (2014) derived the refined N H function and its evolution, where the evolution of absorbed AGN fraction is based on that of Hasinger (2008). Writing down the full form of the N H function is beyond the scope of this paper and the reader is directed to citetueda14. The N H function is normalized as 24 20 f (L X , z; log N H )d log N H = 1 (unabsorbed AGNs with log N H < 20.0 are included in the 20 ≤ log N H < 21 bin).
At z = 0, the N H function of U14, which is based on the 9-month Swift BAT, AMSS and SXDS surveys, is in good agreement with that of Burlon et al. (2011), which is based on the three-year Swift BAT survey.
Also U14 generated a set of sophisticated template X-ray spectra of AGNs, which include (i) an underlying power-law continuum with a high energy cutoff at E c = 300 keV, (ii) a reflection component, (iii) a scattering component from surrounding gas outside the torus, corresponding to 1% with a torus opening solid angle is 2π, and an intrinsic absorption with Compton-scattering derived from Monte-Carlo simulations. For our purpose, we take the spectral energy distribution (SED) templates in Fig. 7. of U14 for log N H [cm −2 ] = 20.5, 21.5, 22.5, and 23.5, but ignored the luminosity and redshift dependence of the SED. In these template spectra, Γ = 1.94 and Γ = 1.84 are assumed for the primary power-law continuum for log N H [cm −2 ] < 22 and > 22 AGNs respectively.
In this section, we refine our definitions. The symbol L X represents an intrinsic rest-frame 2-10 keV luminosity of the primary power-law, which does not include the reflection component. This is because including the reflection component double counts the intrinsic isotropic luminosity of the AGN. The reflection component is typically ∼ 15% of the primary power-law in the 2-10 keV band. The symbol dΦ/d log L X is the comoving number density of AGNs per dex in L X for all Compton-thin (CTN) AGNs (log N H < 24.0).
To feed realistic spectra into the maximum likelihood procedure, the N function, Eq. 6, should be a function of observable quantities for each sample. We define the "observed" 2-10 keV luminosity defined in Eq. 1.
In this case, clearly distinguishing between observed and intrinsic X-ray luminosities, Eq. 6 can be rewritten as: where S X,obs is now a function of N H in addition to L X and z. The integral is over log N H and the intrinsic luminosity L X is a unique function of L X,obs , N H , and z. The ratio of d log L X,obs /d log L X becomes non-unity only if the spectrum changes with L X for a given N H , which is the case for the Swift BAT sample (see below).
In each template spectrum, with an intrinsic 2-10 keV luminosity of the power-law continuum L X , we can calculate the observed luminosity at 2(1 + z) − 10(1 + z) keV (L obs X,tmp (N H , z)). The true 2-10 keV flux from this object, under this template spectrum, is now: The difference between S tru X,obs to S Γ=1.7 X,obs (or the catalogued flux with an ECF calculated with Γ = 1.7 or our assumed spectrum for the ECF) is often neglected in the literature. In reality, it becomes important for highly absorbed AGNs (log N H 23) at low redshifts (z 1) in "2-10 keV" surveys, because the effective areas of the instruments used in generating our sample other than Swift BAT, drop rapidly as we go from 2 keV to 10 keV. Thus the ECF increases rapidly as the spectrum gets harder.
An ideal approach would be to use separate response curve and countrate-survey area relation for each survey, as U14 did, since each survey use different instruments and/or different observational epochs, for which the calibration may vary. However, a limitation of the software that we use for this work is that it can only accept one fluxlimit -solid angle curve, as the software is a legacy from a series of single-band soft X-ray luminosity work Miyaji et al. (2000aMiyaji et al. ( , 2001; Hasinger et al. (2005). While this reduces the accuracy of the estimate of S tru x,obs , the calculation speed gained by this simplification has allowed us to explore analytical forms for various sub-sets and cases (see below). We find that the relative ECFs between those for absorbed AGN spectra and for Γ = 1.7 is not very sensitive to which instrument we use. For example, the ratios of the ECFs for an extremely hard object (Γ = 1.7 with log N NH = 23.5 at z=0.2) to that for an unabsorbed Γ = 1.7 spectrum are 2.0, 2.1, and 1.7 for the XMM-Newton PN+2MOS (countrate at 2-8 kV), the Chandra ACIS-I Cycle 8 (countrate at 2-7 keV) and the ASCA GIS (CR at 2-10 keV) respectively. Here we use the response for ACIS-I for the conversion between S tru x,obs to S Γ=1.7 x,obs . A caveat is that this causes some inaccuracy for the HBSS (CR is at 4.5-7.5 keV) and XMM-Newton LH surveys (CR is at 2-4.5 keV). However, the impact is relatively minor, since these latter surveys contain relatively small fraction of sources in the relevant flux range, i.e., the HBSS is overwhelmed by the ASCA MSS+LSS and the LH sample by the COSMOS sample in their respective flux coverages.
For the Swift BAT sample, the selections are made using the 15-55 keV range, which is very little affected by the absorption up to log N H = 24. Despite the limitation of our software described above, the Swift BAT sample can be treated separately because the flux regime covered by this sample does not overlap with any other. We have made the conversion between the 2-10 keV and the 15-55 keV bands as follows.
In the comparison between Swift 14-194 keV and de-absorbed MAXI 2-10 keV luminosities, Ueda et al. (2011) pointed out that the effective power-law index between these bands varies with luminosity from Γ l = 1.7 to Γ u = 2.0. We have confirmed this by comparing the 14-55 keV luminosities from the (Ajello et al. 2012) catalog with 4-10 keV luminosities in the latest MAXI catalog of Hiroi et al. (2013). To approximately represent the change, we use a smoothly varying luminositydependent effective spectral index Γ eff = Γ l + 0.5 erf[(log L X − 44.0)/2.0](Γ u − Γ l ) (11) for the conversion between the 2-10 keV and 15-55 keV fluxes, where erf is the error function. U14 hypothesize that the main cause of this slope change is the difference between type 1 and type 2 AGNs in the strength of component reflected by the torus (i.e., the torus tend to be thicker in type 2 AGNs) in combination with the relation that the type 2 fraction decreases with luminosity. Because of the definition of L X in this section, the 2-10 keV flux is calculated based on (1 + 0.15)L X following the assumption that 15% of intrinsic 2-10 keV luminosity is added by the reflection component. A small K-correction has been made based on Γ eff . In order to use our flux limit -survey area curve, we have back-converted from the 15-55 keV flux to observed flux S Γ=1.7 X,obs using our reference slope Γ = 1.7.
As a technical note, we comment that U14 chose to use the count rate as an independent variable to the N-function instead of log L X,obs .

The N obs /N mdl method for the Binned XLF
A less biased way of presenting the binned XLF than the classical 1/V a (Avini & Bahcall 1980) is the N obs /N mdl estimator described in Miyaji et al. (2001) (see also La Franca & Cristiani 1997). The basic procedure is: 1. Divide the combined sample into several redshift shells. For each redshift shell, fit the AGN XLF with a smooth analytical function obtained with the maximum-likelihood fit as described above.
2. For each redshift shell, count the number of objects that fall into each luminosity bin to obtain the observed number of objects (N obs ).
3. For each luminosity bin, evaluate the analytical fit at the central luminosity/redshift (dΦ mdl /dlogL X ).
4. Calculate the predicted number of AGNs in the bin (N mdl ).

The final result is:
We note that, in the case that the XLF model is constant within the bin, the N obs /N mdl estimator is mathematically equivalent to that proposed by Page & Carrera (2000). In the case of taking absorption into account with the method described above in Sect. 4.3, and binning is by log L X defined there, deriving N obs is not straightforward, since we cannot determine log L X individually from single-band fluxes. The N H values derived from hardness ratios are subject to large errors and not necessarily available for all objects in the sample. On the other hand, it is straightforward to calculate from the model probability distribution of N H for given redshift and observed flux: p(z, S x,obs ; log N H ), which can be normalized as: This is different from the N H function f (L X , z; log N H ), which is the distribution of N H for given redshift and intrinsic luminosity. For our purpose, for the i-th object with a given (z i , S x,obs,i ), we use the probability distribution of the intrinsic luminosity: If we would like to calculate N obs in the redshift-intrinsic luminosity bin (z 1 , z 2 ; log L 1 , log L 2 ): where i is over objects with redshifts between z 1 and z 2 . In practice, we calculate N obs as weighted sums over log N H = 20.5, 21.5, 22.5 and 23.5 templates. Similarly, we include the PDFs of photometric redshifts into N obs estimations by adding the portion of the PDF of each objects that falls into the redshift bin, as has been made also by Aird et al. (2010). One limitation is that the errors assuming Poisson statistics are not valid if N obs is a sum of weights, unlike the cases of our previous works (Miyaji et al. 2001;Hasinger et al. 2005, U14), in which N obs were always integers. However, we use the approximate 1σ Poisson errors from the Gehrels (1986) Eqs. (7) and (12) (S = 1) for the upper and lower 1σ errors respectively. Thus, the errors are only approximate. Also the errors in different bins are not independent of one another.

The Luminosity-dependent Density Evolution Model
For the analytical expression of the XLF at z=0, we use the smoothed two power-law formula.
The normalization A z=0 44 is the XLF value at (log L X , z) = (44, 0). There is also a convention to use A * (2nd line in Eq. 16) as a normalization parameter. We primarily use A z=0 44 , because A * is strongly coupled with the break luminosity L X, * , and therefore it is not possible to estimate its error accurately with our maximum-likelihood fitting procedure (see below). Following the most successful analytical form of the evolution in the literature, we use the luminosity-dependent density evolution: where e d (z, L X ) is the density evolution factor normalized to z=0. We describe it with a threesegment power-law model: As in the previous works, the first break positions z p1 may depend on luminosity at the low luminosity regime: We also include the luminosity dependence of the evolution indices p1 and p2: p 1 (L X ) = p 1,44 + β 1 (log L X − 44) (20) p 2 (L X ) = p 2,44 + β 2 (log L X − 44). (21)

Global LDDE Results
The best-fit parameters of our LDDE model and 1σ errors, corresponding to ∆L = 1, are shown in Table 2. For the global expression, the fits using the likelihood function Eq. 5 give much more stable error search results than those using Eq. 4. As stated above, minimizing L in Eq. 5 cannot determine the normalization. Thus the A z=0 44 value and its error are determined by scaling it such that it gives the total number of objects (N tot ) for best-fit values of other parameters. The fractional error of 1/ √ N tot is used for its error. Thus the errors of A z=0 44 in Table 2 do not contain correlation of errors with other parameters. The traditional normalization A z=0 * (see Eq. 16) is also shown without an error for reference.
The parameters for the fit in which the realistic AGN spectra, N H function (Sect. 4.3), and the PDFs of photometric redshifts (Sect. 4.2) are fully taken into account are listed under "Full" in Table 2. For our final results, we show the case where z b1,44 is fixed to 1.1 and z b1,44 = 2.7, since fitting with all parameters free causes failures in the minimization processes. The values of these fixed parameters have been decided based on the results of the fits in luminosity-class divided samples, which is explained later in Sect. 6.2. Since the fits in each luminosity class also use information on the global expression, thus the process has been iterative, i.e., we decide on fixed parameters of the global expression referring to the results of the luminosity-class divided sample based on the global expression in the previous step.

Adaptive binning
In this section, we present the binned XLF, as well as analytical formula of the XLF divided into redshift shells. We also present the number density curve divided in luminosity classes. The binning has been decided on a adaptive-basis, so that the sizes of error bars and bin-sizes are optimal for the plots to be visible and informative.
First, we divide the sample into 11 redshift shells in the same way as those used by U14, except that the two highest redshift shells by U14 (z > 3.0) are combined into a single one for obtaining a sufficient number of objects for a separate two power-law fit. For each redshift shell, we construct luminosity bins such that each contains 20 objects by default, with minimum and maximum bin sizes of ∆ log L X = 0.125 and 0.5 respectively. In case the most luminous bin contains less than 20 objects, the border between the last two bins is adjusted to contain approximately equal numbers of objects within the bin size constraints.
For the luminosity division, the range 42.0 ≤ log L X ≤ 46.0 is divided into four ∆ log L X = 1.0 classes. Each luminosity class is divided into redshift bins such that each bin has 20 objects (using the best fit photometric redshifts for those without spectroscopic redshift) with forced minimum and maximum sizes of ∆ ln(1 + z) = 0.07 and 0.5 respectively. If the highest redshift bin has less than 20 objects, the border between the highest two redshift bins are adjusted in the same way as the case of luminosity binning. Due to the small number of objects, the default number of objects per redshift bin for the 45.0 ≤ log L X ≤ 46.0 class has been set to 7, with the same minimum and maximum bin size constraints.

Analytical form for each redshift shell
For the analytical expression of the XLF in each redshift shell, we use the smoothed two power law formula. Because the redshift shells have finite widths, the fit results depend on the evolution of the SXLF within them: where e d (z, L X ) is the density evolution factor. While the final results are insensitive to the detailed behavior of e d (z, L X ) within the shell at most locations in the (L X , z) space, we have taken our best-estimate by using the luminositydependent density evolution (LDDE) model derived above in Sect. 5. The luminosity range of the fit is from log L X,obs =41.5 to 46.0. The best fit parameters of Eq. 22 for each redshift shell are shown in Table 3. The normalization is defined by: where z c is the central redshift of the shell, which is defined by z c + 1 = (1 + z min )(1 + z max ), where z min and z max are the minimum and maximum redshifts of the bin. One sigma parameter errors are calculated in the same way as in Sect. 5.2. In this section, we have been able to obtain 1σ errors of most parameters with Eq. 4, and unlike the case of Sect. 5.2, the errors of A z=zc 44 contain the correlations among parameter errors. We also show the traditional normalization A z=zc * without errors. Table 4 tabulates the binned N obs /N mdlestimated XLF for the full treatment case, along with the observed number of AGNs (N obs ) and the final estimated values XLF value at the center of each bin. The full XLFs in the 11 redshift is for all other parameters fixed.
d The traditional normalization defined as the twice of the XLF value at the break luminosity at z = 0 (Lx, * ). shells are plotted in Fig. 4 in separate panels. The best-fit smoothed two power-law model (2PL) and the best-fit LDDE model evaluated at z c are overplotted to assess the goodness of the models. The best-fit two power law models at z c = 0.1 (lowest redshift bin) and at z c = 1.8, where the A 44 is near the peak, are also shown in all panels for reference.      Figure 5 shows the variation of the smoothed two-power law parameters as a function of redshift. The normalization is the XLF at log L X = 44 and we see the familiar feature of growth (with increasing redshift) between z = 0 and z ∼ 1, followed by a plateau, then a decline at z > 3. A notable feature is a sudden drop of the low luminosity slope γ 1 from z < 0.6 (γ 1 ∼ 1), to z ≥ 0.6 (γ ∼ 0.5). No significant change of slope is observed in each side of this redshift. An F-test comparison between the best-fit models for the case where γ 1 is constant over redshift versus the case where γ 1 changes at z = 0.6 shows that the probability that the former is accepted is only 7 × 10 −5 . The high luminosity slope is consistent with being constant at γ 2 ∼ 2.7, except for two data points at z ∼ 0.7 and z ∼ 1.1, where γ 2 is significantly smaller, although it may be caused by an observational bias as discussed in Sect. 7.

Evolution of the Space Density in Luminosity Classes
In this section, we investigate the evolution of the AGN space density in different luminosity classes as a function of redshift. We have made the fit to an analytical expression at the center of each of the four luminosity classes: 42.0-43.0, 43.0-44.0, 44.0-45.0, and 45.0-46.0 in log L rmx . The XLF at the central luminosity of the bin log L c , which is defined as the mean of the minimum and maximum bounds of log L x bin can be expressed as: where the normalization A z=0 Lc is the XLF value at the central luminosity of the class at z = 0.
For each luminosity class, we have made a maximum-likelihood fit to the the redshiftdependent evolution parameters in the evolution factor expressed as a three-segment power-law model Eq. 18, evaluated at L X = L c .
During the fitting process, the luminosity dependence of the XLF within the shell is fixed to the best-fit LDDE case. In each class, there are six fitting parameters: A z=0 Lc , p 1 , z b1 , p 2 , z b2 and p 3 . As in the case of redshift shell divided samples, maximum-likelihood fits have been made for each luminosity class using Eq. 4, where the normalization is also a fitting parameter. Thus the  reported normalization errors contain the effects of correlations with other parameters. The resulting best-fit parameters and errors are reported in Table 5. One limitation of dividing the sample into luminosity classes is that we do not know the intrinsic luminosity L X of each object. We only know the observed flux and redshift. For the fit of the luminosity class (log L X,min , log L X,max ), we select objects that fall into the log L X,min ≤ log L X,obs ≤ log L X,max range. Thus the fitting in each luminosity class uses some objects that are outside of the log L X,min ≤ log L X,obs ≤ log L X,max range, while others that fall into this range are not used. The fitting process using Eq. 9 properly takes care of the expected number of AGNs in the observed luminosity range calculated based on the intrinsic luminosity XLF model. Figure 6 shows the evolution of the AGN number density and emissivity in each class for the full treatment. The binned number densities and emissivities, calculated with the N obs /N mdl method using the luminosity-class divided models, are shown with 68% Poisson errors. The total emissivities of the 42 ≤ log L X [h 70 erg s −1 ] < 46 range are also plotted as a function of redshift. The model number densities and emissivities are also shown for the best-fit LDDE model in Table 2 as well as class-divided fits. The evolution of total emissivity in 42 < log L X < 46 is also plotted. Table 6 tabulates the binned N obs /N mdl -estimated values evaluated at the central luminosities L c (instead of number densities integrated over the luminosity class, which are plotted in Fig. 6), as well as emissivities with best fit models.
The results in Table 5 show that the first break redshift z b1 increases with luminosity with z b1 ∼ 1.1 at log L X ≈ 44[h 70 erg s −1 ]. The second break redshift does not vary significantly and is consistent with z b2 ∼ 2.7 at log L X > 43 [h 70 erg s −1 ]. The slope after the second break is consistent with p 3 ≈ −5.6. These have been fed back as fixed parameters for a refinement of the global expression in Sect. 5. For the 42 ≤ log L X < 43 class, the values of z b2 and p 3 could not be determined, because the sample does not extend to the second break z b2 ∼ 2.7. For the 45 ≤ log L X < 46 class, the low number of objects has limited our analysis. While we obtain best fit values of z b1 ≈ 1.5 and z b2 ≈ 3, their 1σ errors are essentially unconstrained. Thus we fixed the z b1 = 1.2 and z b2 = 2.6 based on the results of the fit to the 44 ≤ log L X < 45 class. We have also attempted a two-segment law fit (one break redshift) for this class, but the three segment form with the two break redshifts fixed at these values gives a better fit.
While the value of the parameter p 1 increases with luminosity, the measurements are over different redshift ranges for different luminosity classes and should not be compared with one another. The p 2 parameters are also measured over somewhat different redshift ranges. However, it is clear from Fig. 6 that the slope between z b1 and z b2 goes from negative to positive (as a function of log(1 + z)) with luminosity. The full binned XLF divided in luminosity classes is shown in Table 6. Table 5 Best-fit evolution parameters for each luminosity class log L X -range a log Lc a A z=0  Table 3 apply.
b Fixed parameters are indicated by an asterisk ('*').   Figure 6(a) clearly shows a shift of the number density peak with luminosity up to log L X ≈ 44.5 in the sense that more luminous AGNs (QSOs) peak earlier in the history of the universe, while the low luminosity ones arise later. Also, there is a clear decline of the derived space densities at all luminosity classes. A notable feature is a clear plateau structure at intermediate redshifts at intermediate luminosities, which we have been able to trace in detail thanks to the addition of the COSMOS dataset. This has been reflected in the total emissivity curve, which stays almost constant between 1 < z < 2.7. The slope between the first break z b1 and the second break z b2 changes from negative to positive.

Sensitivity to pz-PDZs and Spectral Models
For comparison, we show the parameters for the case where we only take the best-fit single photometric redshift (PZ) per object, instead of full pz-PDFs, under "No PDZ" for the global expression ( Table 2). The fitted parameters are not significantly different between the two cases except p2 44 , which is the slope between the first break and the second break redshifts at log L X = 44. The difference in the normalizations are more than the errors given, but these errors are underestimates, since, for the global expression, we use Eq. 5 and correlations with other parameters have not been taken into account, as explained in Sect. 4.1 and 5.2. To observe the effects more closely, we plot the ratio of the XLF estimates based on the full PDF and no PDF samples for each redshift shell in Fig. 7. The same are plotted for the number density evolution in each luminosity class in Fig. 8. Upon the examining the difference between the best-fit models in these figures, the effects of the pz-PDF seem to be as large as ∼ 30% at some locations (blue dashed lines). However, comparing the data points from the N obs /N mdl (black filled hexagonal versus blue open triangle data points in these figures), we see that there is no significant difference and the apparent discrepancy of the models are at locations that are not well-constrained by the data. Therefore the apparent deviation is not caused by the pz-PDFs, but by the different convergence of of the model parameters at the level of a fraction of statistical errors. This is an optimistic case, because in the COSMOS survey, where the majority of the photometric redshifts with the pz-PDZ came from, has exceptionally high quality photometric data, achieving an accuracy of ∆z/(1 + z) ≈ 0.015 with ∼ 5% outliers (Salvato et al. 2011). At the time of our analysis, full pz-PDFs were not available for the CLANS, CLASXS, CDF-N and CDF-S data in our disposal, in which the spectroscopic redshift completeness ranges from 60-80%. These fractions are similar to those of the LH and COS-MOS dataset. For these, we take best-fit photometric redshift values, or for CDF-S sources, if the information about the second peak is available, (Luo et al. 2010;Xue et al. 2011), we give a weight of 0.5 in each peak for the full analysis. The photometric redshifts of the CDF-S have a normalized median absolute deviation (NAMD) of σ NMAD = 0.059 with outliers of 8.6% (Luo et al. 2010), indicating a similar quality to those of COS-MOS. Those of the CLANS, CLASXS and CDF-N surveys use fewer photometric bands and the uncertainties are as large as σ z /(1 + z) ∼ 0.16, but with ∼ 6% outliers (Trouille et al. 2008). The uncertainties of these PZ have not been taken into account in our full analysis. For a further test, we create fake PDZ's for objects with photometric redshifts (except those with a second peak) in these surveys with six delta functions. We assume a Gaussian profile centered at the cataloged PZ value z p for the main pz-PDF peak with an appropriate width corresponding to the σ NAMD or σ values for the each survey. In our fake PDZ, the main peak is represented by three delta functions, one at the center (representing the central 1/2 of the PDF main peak), the others at the 12.5% and 87.5% percentiles of the profile (representing the upper and lower 1/4 of the PDF main peak). The central one has a weight that is half of 0.5(1 − outlierfraction) and the other two each has a weight that is half of the central one. Also added are three delta function to represent outliers. Each of the three has a weight of 1/3 of the outlier fraction and the redshift is taken randomly from spectroscopic redshifts of X-ray sources that have fluxes within a factor of three of that the object. We keep using the real PDZs for the COS-MOS and LH samples. We make global fit the sample and all the fitted parameters are within a fraction of 1σ errors from the original full analysis. Thus we conclude that the uncertainties of the PZ's have minimal net effect to our XLF analysis.
In the literature, it has been often assumed for the "hard" (2 ≤ E[keV] ≤ 7 − 10) surveys that absorption is negligible and the HXLF is calculated without taking absorption into account, simply by assuming a power-law spectra (e.g., Silverman et al. 2008;Yencho et al. 2009;Aird et al. 2010). For comparison, we also compute the case where all the AGNs have simple a Γ = 1.8 power-law. The results of this case, with PDFs of the photometric redshifts are taken into account, are listed under "Γ = 1.8" in Table 2 and the comparisons with the "Full" case are overplotted in Figs. 7 and 8. As seen in these figures, the XLFs based on the simple spectral assumption of Γ = 1.8 without absorption correction is subject to underestimates of the XLF at low luminosities (log L X 44) at low redshifts (z 1), by up to 50%. If one uses a sample selected at E 5 keV, however, this systematic error is expected to be reduced dramatically (S. Fotopoulou et al. 2015, submitted to A&A).
In principle, the results depend on the N H function and its evolution. Using the 14-195 keV selected Swift BAT AGNs with detailed X-ray spectroscopy at E < 10 keV and detailed modeling of selection effects, the N H function derived by U14 is robust, especially at low redshifts, and in agree-ment with that of Burlon et al. (2011) derived using a larger sample of Swift BAT AGNs. The fraction of absorbed AGNs among CTN AGNs is observed to increase with redshift (e.g. Hasinger 2008). However, it is not clear whether the increase of absorbed fraction saturates at z ∼ 2 or keeps increasing to higher redshifts, and Hasinger (2008) commented that a model where the absorption fraction keeps increasing at z > 2 is also marginally acceptable. In our analysis above, we have used the U14 result, where the fraction of 22 < log N H < 24 AGNs among CTN AGNs at log L X = 43.75 grows as ∝ (1 + z) 0.45 at (z < 2.0) and stays constant at higher redshift. For a sensitivity check, we assume that the absorbed fraction grows as ∝ (1+z) 0.32 up to our sample limit, where the index has been chosen to have the same absorbed fraction at z = 3. The results of the global fit under this assumption have not changed significantly, and all the parameters are consistent with our "Full" case within a fraction of 1σ. This can be understood because, at high redshift, the 0.5-10 keV band, which is usually used for the spectral analysis or hardness ratio analysis to measure N H , corresponds to rest-frame E > 1.5 keV at z > 2. Thus it is already not sensitive to measuring absorption. For the same reason, from the 2-10 keV flux, we measure rest frame higher energy luminosity, which has not been affected by absorption. Since we want to derive the intrinsic XLF, the measured X-ray luminosity is insensitive Fig. 7.-The comparison of best fit models in each redshift shell and the N obs /N mdl estimates. The best fit two power-law models of the case of the no pz-PDFs and Γ = 1.8 are divided by those of the full case. The N obs /N mdl XLFs are also divided by the full model. See labels in the lowest rightmost panel for the meanings of line styles, symbols, and colors (electronic version). The data points of the no PDF and Γ = 1.8 cases are shifted horizontally by +0.04 and −0.04 respectively in log L X for display purposes. (Color version of this figure is available in the electronic version.) to the details of the N H function at high redshifts.

Detailed behavior of the XLF
Thanks to the addition of the COSMOS dataset, we have been able to trace the detailed behavior of the XLF in the intermediate redshift/intermediate luminosity range. We describe a number of notable features.
• The low luminosity slope γ 1 flattens suddenly at z ∼ 0.6. The slopes stay almost constant within each of the z < 0.6 and z ≥ 0.6 regimes, having γ 1 ≈ 1 and ≈ 0.5 respectively. This sudden flattening of the slope practically excludes ILDE/LADE models, where the shape of the XLF stays unchanged over redshift, as a global expression. The high redshift faint end slope is consistent with those of optically selected QSOs at 1 z 3.5 (γ 1 ∼ 0.45) measured by Bongiorno et al. (2007) and z ∼ 4 QSOs by Ikeda et al. (2011) (γ 1 = 0.7±0.2). Note, however, that the faint end slope measured by Glikman et al. (2010) is significantly steeper (γ 1 = 1.3 ± 0.2 at z ∼ 4).
• The high luminosity slope stays almost constant (γ 2 ≈ 2.7), with the exception of two redshift shells at 0.6 < z < 0.8 and 1.0 < z < 1.2, which have flatter slopes γ 2 ≈ 2.0. For the shell between these (0.8 < z < 1.0), γ 2 is consistent with those for the rest of the shells. We note that the high luminosity end of the XLF is poorly sampled at these redshifts and the apparently low γ 2 values may be caused by the higher weights at luminosities just above the break and less weight at highest luminosities. Thus there is no evidence for systematic change of the high luminosity slope with redshift. The upcoming eROSITA mission will provide better sampling of this regime, which will enable us to see if the high end slope is constant over all redshifts.
• The addition of the new data including COSMOS revealed new details about AGN downsizing. The number density/emissivity growth curves have two distinctive breaks, one at z b1 ≈ 1 and the other at z b2 ≈ 2.7, with a "plateau" between them. In any luminosity class, a rapid rise of number density with cosmic time is observed at z b2 (except the lowest luminosity bin) and a rapid decline is observed after z b1 . At low luminosities, (log L X 44), the peak is at z b1 and high luminosities, (log L X 44), z b2 becomes the peak. Also the position of the first break moves from lower to higher redshifts with luminosity up from z ∼ 0.6 to z ∼ 1.4 up to log L X 44.5. No evidence for a change in z b2 with luminosity is observed.
Our results are in excellent agreement with those of U14 as seen in Fig. 9(a). Small excesses at the level of 0.2 dex at a few data points can be well explained by cosmic variance. Especially our a few × 10 15 S X,obs 10 −14 erg s −1 cm −2 range is dominated by the COSMOS data, which may be affected by large scale structures seen as redshift spikes Brusa et al. 2010;Civano et al. 2012). Other possible differences at the same level may be explained in the treatment of the probability density distribution of the photometric redshifts.
Our results at the highest redshift shell z > 3 as well as the redshift evolution factor (p 3 ) are consistent with the recent results by Vito et al. (2014).

Comparison with Theoretical Models
There are a number of attempts to reproduce/understand the down-sizing behavior of AGN evolution within cosmological simulations and semi-analytical models. In order to understand the implication of our results to the physical scenarios on accretion processes by comparing observations and theory, we should consider which features of the theoretical predictions are a consequence of the assumed physical picture and which features are a consequence of adjusting free parameters to match the AGN luminosity functions in various redshifts.
In their attempt to reproduce the faint end AGN LF, Degraf et al. (2010) considered a BH growth by infall of surrounding gas via Bondi, H. (1952) accretion, which some authors call "radio mode" or "hot halo mode", and merging with other BHs in their simulations including hydrodynamic gas component. They compared their results with luminosity-class divided number density evolution curves of soft  and hard X-ray . They recognize that their model overproduces number densities of X-ray selected AGNs at z > 1. We compare our number density evolution with two semi-analytical models of Marulli et al. (2008) and Fanidakis et al. (2012), where we have been able to convert their published results to the density curves in our 2-10 keV luminosity classes in Fig. 9(a). This figure shows that the overproduction is also seen in these semi-analytical models (SAMs).
Since these authors did not provide the number density curves in the 2-10 keV range, we explain the derivation of the luminosity-class divided 2-10 keV number density curves from their published figures. For Marulli et al. (2008), we use their Fig. 7 and take the "best" bolometric luminosity function model. The bolometric LF has been converted to 2-10 keV XLF model using the bolometric correction by Hopkins et al. (2007) and integrated over our luminosity classes. For Fanidakis et al. (2012), we use their soft (0.5-2 keV) X-ray "before obscuration (total)" curves provided in their Fig. 18 (b) instead of their "after obscuration (visible)" to compare with our 2-10 keV density evolution curves. We convert their "before obscuration" to those for 2-10 keV assuming a photon index of Γ = 1.9, which is representative of unobscured AGNs. Since they use slightly different cosmological parameters, we also convert the density curves to our adopted cosmology to overplot in Fig. 9(a).
In both cases, we see overpredictions of AGN number density at highest redshift at low luminosity bins. As AGN model components, the Marulli et al. (2008) model took into account the "radio" mode accretion and merger driven "quasar" mode. The Fanidakis et al. (2012) model took into account the "radio" and the "starburst" modes, where the "starburst mode" includes merger-driven "quasar mode" and secular (disk instability) components. Marulli et al. (2008) considered a number of AGN light curve models and searched for the one that fits well with the bolometric LF of Hopkins et al. (2007), thus had freedom in adjusting to observations. On the other hand, Fanidakis et al. (2012) had four free parameters to fit to the LF. Except for the duty cycle, which is adjusted to the global normalization of the LF, the parameters are physically motivated. The Fanidakis et al. (2012) model overpredicts the number densities at high redshift at all luminosity classes. In their comparison with the 0.5-2 keV XLF, they took into account the effects of luminosity dependence and redshift evolution of obscured/type 2 AGN fraction f 2 by Hasinger (2008). For the redshift evolution of f 2 , Hasinger (2008) proposed a model with f 2 ∝ (1 + z) 0.62 up to z = 2.06 and which stays constant at higher redshifts as preferred, while Fanidakis et al. (2012) used a marginally acceptable model with f 2 ∝ (1 + z) 0.48 without saturation, and thus their obscured fraction keep increasing to higher redshifts. The latter form of the obscured fraction evolution was needed in the model of Fanidakis et al. (2012) for satisfactory consistencies with the optical QSO and soft XLF to their model. On the other hand, the high-z overprediction of the 2-10 keV number densities show that they indeed over predicted high redshift AGN population and over-corrected for the effects of obscuration to match with the 0.5-2 keV number density/optical QSOs. The tendency of overprediction of high redshift low luminosity/mass population is also a problem with predicting stellar mass function of galaxies (e.g. Guo et al. 2011), the reason for which is still unclear. On the other hand, if their basic prediction on high redshift accretion is correct, the discrepancy may be an indication of a large population of highly Compton-thick AGNs at high redshifts. The high redshift over-prediction problem is also present in the recent model by Sijacki et al. (2014) at all luminosities.
Another interesting model is by Draper & Ballantyne (2012), where they considered merger and secular modes of AGNs triggering. Their work was not directly based on cosmological simulations. Instead they took simple formulae for the evolution of galaxy number density with stellar mass above some threshold value, gas fraction, and merger rate from the literature. They then considered secular and merger driven AGN triggering followed by a parameterized model of luminosity decay. Their model 2-10 keV number density evolutions are overplotted with our data in Fig. 9(b), which has been derived from their Fig. 7. Since their Fig.  7 shows curves for log L X > 42, > 43, and > 44, we draw their model curves of log L X = 42 − 43 and 43 − 44 curves by subtracting their > 43 curve from > 42 one and > 44 curve from > 43 one respectively. For the 44−45 curve, we use their > 44 curve, since the contribution of the log L X > 45 AGNs to the number density is only ∼ 0.2% of > 44 AGNs for γ 2 = 2.7. We see that their "sum" of the contributions of the two processes is in good agreement with our data. Since they have three free parameters to adjust in their AGN light curves, which are the timescale, slope of the decay, and the peak Eddington ratio, they had much freedom in adjusting to observed XLFs. Thus it is no surprise that their theoretical curves agree well with our observation. However, they had difficulties in reproducing the hard X-ray number density curve with only one mode, i.e., merger only or secular only. Their summed curves show the flat-top (two break) structures, where the slope in the z b1 < z < z b2 range increases from negative to positive with luminosity.

Conclusions
We have investigated the detailed behavior of the intrinsic 2-10 keV X-ray luminosity function (XLF) of active galactic nuclei (AGNs) from a combination of samples in various depths/width. We summarize our conclusions.
1. A total of ∼ 3200 unique X-ray AGNs from a combination of wide to deep samples covering the redshifts of 0.015 < z < 5.8 and six orders of magnitude in flux have been used to investigate the behavior of the XLF. In particular, we present here the XLF from X-ray sources from both of the XMM and Chandra COSMOS surveys, which comprise about 43% of our global sample.
2. Recent N H functions and realistic spectral templates are fully incorporated in the maximum-likelihood fitting to analytic expressions as well as the N obs /N mdl -estimated binned intrinsic XLF. The photometric redshift probability density distributions (pz-PDFs) of the COSMOS and Lockman Hole samples as well as secondary photometric redshifts of CDF-S are also considered.
3. We present a full parametric luminositydependent density evolution (LDDE) model, and redshift-shell and luminosity class sep- arated models, as well as the binned XLF using the N obs /N mdl estimator.
4. The low luminosity end slope of the XLF flattens suddenly at z > 0.6. This behavior practically excludes any global expression of XLF evolution that assumes an unchanged XLF shape over cosmic time.
5. We investigate the net effects of ignoring probability density distribution (pz-PDF) of photometric redshifts. The effect of including pz-PDF into the analysis does not alter the final XLF results significantly.
6. We investigate the effects of assuming a simple AGN spectrum of Γ = 1.8. Under this simple assumption, the XLF is subject to underestimates by up to ∼ 50% in the low redshift low luminosity regime.
7. The detailed behaviors of AGN down-sizing have been revealed. A clear two-break structures has been revealed in the number density evolution curves at every luminosity class above log L X > 43. This behavior is qualitatively consistent with a two-mode AGN evolution involving major merger and secular processes.
8. Most current semi-analytical models of AGNs over-produce AGN number densities in the high redshift, low luminosity regime. Their semi-analytical treatment of accretion processes may have to be revised to reproduce our results. Alternatively, if their accretion scenario is correct, this may suggest a population of heavily Compton-thick AGNs in this regime. . We gratefully acknowledge the contributions of the entire COSMOS collaboration consisting of more than 100 scientists. More information about the COSMOS survey is available at http://www.astro.caltech.edu/∼cosmos. The authors thank Len Cowie for allowing us to use their unpublished spectroscopic redshifts and Murray Brightman for his advice on CTK AGNs.