Constraining the epoch of reionization from the observed properties of the high-z Universe

We combine observational data on a dozen independent cosmic properties at high-$z$ with the information on reionization drawn from the spectra of distant luminous sources and the cosmic microwave background (CMB) to constrain the interconnected evolution of galaxies and the intergalactic medium since the dark ages. The only acceptable solutions are concentrated in two narrow sets. In one of them reionization proceeds in two phases: a first one driven by Population III stars, completed at $z\sim 10$, and after a short recombination period a second one driven by normal galaxies, completed at $z\sim 6$. In the other set both kinds of sources work in parallel until full reionization at $z\sim 6$. The best solution with double reionization gives excellent fits to all the observed cosmic histories, but the CMB optical depth is 3-$\sigma$ larger than the recent estimate from the Planck data. Alternatively, the best solution with single reionization gives less good fits to the observed star formation rate density and cold gas mass density histories, but the CMB optical depth is consistent with that estimate. We make several predictions, testable with future observations, that should discriminate between the two reionization scenarios. As a byproduct our models provide a natural explanation to some characteristic features of the cosmic properties at high-$z$, as well as to the origin of globular clusters.


INTRODUCTION
Reionization of the cosmic gas after it became essentially neutral at a redshift z ∼ 1100 plays a crucial role in the galaxy formation process. Unfortunately, the details of the epoch of reionization (EoR) are to a large extent unknown.
The absence of global absorption shortwards of the rest-frame Lyman-α (Lyα) line in the spectra of quasars at z < 3 made Gunn & Peterson (1965) realize that the hydrogen present in the nearby intergalactic medium (IGM) was ionized (x HI 10 −4 ) except for small intervening systems yielding discrete absorption lines, the so-called Lyα forest. The Gunn-Peterson trough caused by neutral intergalactic hydrogen was finally found by Becker et al. (2001) and Djorgovski et al. (2001) in the spectra of quasars at z ∼ 6, suggesting that value for the redshift z ion,H of full hydrogen ionization.
Several subsequent studies using: (a) the mean opacity of the Lyα forest (Fan et al. 2006); (b) the size of the proximity zone around quasars (Wyithe et al. 2005;Fan et al. 2006;Bolton & Haehnelt 2007a;Lidz et al. 2007;Maselli et al. 2007); (c) the detection of damping wing absorption by neutral IGM in quasar spectra (Mesinger & Haiman 2004, 2007Mortlock et al. 2011b); (d) its non-detection in the spectra of gamma ray bursts (Totani et al. 2006;McQuinn et al. 2007); (e) the abune.salvador@ub.edu dance of Lyα emitters (LAEs) (Malhotra & Rhoads 2004;Haiman & Cen 2005;Furlanetto et al. 2006;McQuinn et al. 2007;Mesinger & Furlanetto 2008;Ouchi et al. 2010;Jensen et al. 2013); and (f) the covering fraction of dark pixels in the Lyα and Lyβ forests (McGreer, Mesinger, & Fan 2011) confirmed a value of z ion,H of 6.0 +0.3 −0.5 . Such a value is not inconsistent with the detection of LAEs beyond z = 6. Not only is the frequency of Lyα photons shifted in origin due to outflow velocities of the emitting gas clouds, but they are also sufficiently redshifted when leaving the ionized bubbles around galaxies for not being absorbed by the neutral hydrogen present in the IGM. In any event, the LAE abundance rapidly drops beyond z = 6 (Hayes et al. 2011), the most distant spectroscopically confirmed LAE lying at z = 6.96 (Ota et al. 2008; but see Zitrin et al. 2015). This would imply that the volume filling factor of ionized regions, Q HII (z), is decreasing at z = 7.
An alternative way to estimate z ion,H is from the analysis of the large-scale cosmic microwave background (CMB) anisotropies. Assuming instantaneous reionization, one can determine z ion,H and the optical depth to Thomson scattering of free electrons by CMB photons, τ (hereafter simply the CMB optical depth). The 9-year data gathered from the Wilkinson Microwave Anisotropy Probe (WMAP9, Hinshaw et al. 2013) led to z ion,H = 10.6 ± 1.2 and τ = 0.089 ± 0.014. The estimates drawn from the temperature three-years Planck data combined with the polarization WMAP data were z ion,H = 9.9 +1.8 −1.6 and τ = 0.078±0.014 (Planck Collaboration et al. 2015), and the most recent ones inferred taking into account the own Planck data on the large-scale polarization anisotropies have decreased to z ion,H = 8.8 +0.9 −0.8 and τ = 0.058 ± 0.012 (Planck Collaboration et al. 2016).
The marked difference between the latter z ion,H estimates and the previous ones drawn from the spectra of distant luminous sources seems to indicate that reionization, far from being instantaneous, was extended, and possibly non-monotonous (Wyithe & Loeb 2004;Cen 2003;Chiu, Fan, & Ostriker 2003;Hui & Haiman 2003;Sokasian et al. 2004;Naselsky & Chiang 2004;Gnedin 2004;Furlanetto & Loeb 2005;Wyithe & Cen 2007;Bolton & Haehnelt 2007b). In this sense, the larger values of z ion,H inferred assuming instantaneous reionization would just indicate that the ionized fraction was quite low by z ∼ 10. Of course, the uncertainty in the associated τ value is also large. However, provided the ionized fraction were small enough at z 15, the value of τ should not be too biased ).
The small-scale CMB anisotropies also inform on the duration ∆z of reionization (Zahn et al. 2012;Planck Collaboration et al. 2016). However, its derivation depends on the uncertain correction for dust emission, the dominant signal at the relevant spatial frequencies, and is very sensitive to the assumed morphological model of ionized bubbles.
Another integral quantity that constrains the EoR, similar to τ but more sensitive to low redshifts, is the weak comptonization distortion of the CMB radiation spectrum. Although its measured upper limit of y < 1.5 × 10 −5 (Mather et al. 1990) is a rather loose constraint, it has the advantage with respect to τ of being inferred with no modeling.
But the only way to probe the EoR at every z is through 21 cm hyperfine line observations (see Morales & Wyithe 2010 for a comprehensive review). Unfortunately, the cosmological signal is 5 orders of magnitude smaller than that produced in galaxies, which greatly complicates the technique, and makes it necessary to have a previous good knowledge of the galaxy formation process (Ichikawa et al. 2010;Iliev et al. 2008). Presently, this approach has only allowed one to put a lower limit to ∆z of 0.06 (Bowman & Rogers 2010).
Among the various mechanisms that might contribute to reionization (see e.g. Chen 2007), the most simple one is photo-ionization by luminous sources, namely active galactic nuclei (AGN), normal galaxies with ordinary Population II (Pop II) stars, and first generation metal-free Population III (Pop III) stars. The X-ray background demonstrates that AGN emit one order of magnitude less UV photons than needed at z ∼ 6 (Worsley et al. 2005;Cowie et al. 2009;McQuinn 2012). This implies that either massive black holes (MBHs) are too scarce (Willott et al. 2010a), or their associated AGN are too obscured by dust (Treister et al. 2011), while more abundant mini-quasars associated with stellar black holes would have too short duty cycles (Alvarez, Wise, & Abel 2009;Milosavljević et al. 2009).
Likewise, the star formation rate densities derived from the rest-frame UV luminosity functions (LFs) of normal bright galaxies (McLure et al. 2010;Lorenzoni et al. 2011;Bouwens et al. 2011Bouwens et al. , 2012a show that these objects are insufficient to ionize the IGM by z ∼ 6. Therefore, ionization would be achieved either as a result of a substantial abundance of faint starforming galaxies (Robertson et al. 2015;Ishigaki et al. 2015;Bouwens et al. 2015) or of the still undetected first generation Pop III stars (Sokasian et al. 2004;Shin, Trac, & Cen 2008).
The situation is even more uncertain regarding He II reionization. The He II mean opacity inferred from the Lyα forest suggests that the redshift of complete He II reionization, z ion,He , should be less than or approximately equal to 3 (Ricotti et al. 2000;Schaye et al. 2000;Dixon & Furlanetto 2009;McQuinn 2009;Becker et al. 2011), and probably no smaller than 2 as the AGN emission begins to decline at that redshift (Cowie et al. 2009). AGN are indeed the most plausible ionizing sources responsible of He II reionization as they emit more energetic photons than normal galaxies, and Pop III stars, the other sources of energetic photons, do not form at z < 6.
It is thus clear that to fully determine the EoR we must monitor the formation and evolution of the three kinds of luminous sources and their feedback since the dark ages. The problem with this approach is that some aspects of the physics involved are poorly known and must be treated through numerous parameters. This affects not only semi-analytic models (SAMs) but also cosmological hydrodynamic simulations. Indeed, simulations accurately deal with dark matter (DM), but most of the baryon physics is at sub-resolution scale. Consequently, they use essentially the same recipes and parameters as SAMs in this regard. This problem is particularly hard for simulations because the huge CPU time they involve severely limits the coverage of the whole parameter space. Thus SAMs are the best option.
The observational data on the Universe at high-z are rapidly growing, and it might nowadays be possible to address successfully this problem. That is, it might be possible to adjust all the free parameters of a SAMlike model such as the Analytic Model of Intergalacticmedium and GAlaxies (AMIGA), specifically designed by  to deal with the interconnected process of galaxy formation and reionization. In the present paper we show that, combining the above mentioned information on reionization with the current data on the evolution of the Universe at z > 2, we can constrain, indeed, the freedom of this model, and have an unprecedented detailed picture of the EoR.
The paper is organized as follows. In Section 2 we recall the basic characteristics of AMIGA, paying special attention to its free parameters. The observational data used to adjust those parameters are given in Section 3, and the fitting strategy used is described in Section 4. Our results are presented in Section 5, and discussed and summarized in Section 6.
AMIGA is a very complete and detailed, self-consistent model of galaxy formation particularly well-suited to monitor the coupled evolution of luminous sources and IGM. All processes entering the problem are treated as accurately as possible according to our present knowledge. In particular, AMIGA is the only model of galaxy formation developed so far that includes Pop III stars, normal galaxies, and AGN, together with an inhomogeneous multiphase IGM.
AMIGA deals with the properties of halos and their baryon content by interpolating them in grids of masses and redshifts that are progressively built, as halos merge and accrete, from trivial boundary conditions. This procedure is less computer memory-demanding than the usual method based on Monte-Carlo or N -body realizations of halo merger trees, and enables us to reach very high redshifts and very low masses while maintaining a good sampling over the entire range in redshift and mass, respectively. In this way we can readily integrate at every z the feedback of luminous sources for the halo mass function (MF), and accurately evolve the IGM.
The total number of free parameters in AMIGA is notably reduced in comparison to ordinary SAMs thanks to causally linking physical processes that are usually treated as independent from each other. This greatly diminishes the degrees of freedom and renders the model predictions more reliable. The reader is referred to  for a full description of AMIGA. Here we provide a brief overview focusing on the model parameters.
2.1. IGM Ionizing photons have a very short mean free path in a neutral hydrogenic gas (hereafter the gas outside halos is referred to as the "IGM"). For this reason, luminous sources ionize (and reheat) the IGM in small bubbles around them, which progressively grow and percolate. H II bubbles in turn harbor hotter He III sub-bubbles due to the most energetic fraction of UV photons emitted. Should the ionizing flux become less intense, the singly or doubly ionized bubbles stop stretching, and the surrounding H I or He II regions begin to grow again by recombination (or new dense H I or He II nodes begin to develop after complete reionization). In other words, the neutral, singly, and doubly ionized IGM phases remain well-separated except for small residual fractions.
X-ray photons produced in supernovae (SNe) and at a lesser extent in Pop III stars, ordinary stars in binary systems, and AGN have instead a large mean free path, and form a uniform background that reheats the IGM by Compton scattering. 1 On the contrary, adiabatic cooling due to cosmic expansion, collisional cooling in hot neutral regions, and Compton cooling by CMB photons oppose themselves to the heating. The latter mechanism is particularly effective at high-z where CMB photons are more abundant.
The intergalactic gas falling into halos is shock-heated to their virial temperature. Thus the only halos able to trap gas at any given moment are those with virial temperatures higher than the time-varying temperature of the surrounding IGM. The accretion of gas onto galaxies is calculated in AMIGA as usual in SAMs: taking into account the cooling of the hot intrahalo gas and its freefall to the halo center. (Hereafter the gas trapped within halos is simply referred to as "the hot gas", and the gas inside galaxies as "the cold gas".) However, such a classical scenario strictly holds only in the case of massive halos where the gas infall rate is limited by the cooling rate. In low-mass halos it is limited by the rate at which the halo accretes gas from outside. Since such an accretion is calculated in AMIGA independently of the symmetry of the problem, our treatment thus implicitly deals with accretion of cold gas onto galaxies according to the filament scenario (Dekel et al. 2009).
Not only do galaxies accrete gas, but they also lose it through SN-and AGN-driven winds. Thus the hot gas continuously incorporates metals through three channels: i) via halo mergers, through the gas within the merging halos; ii) via accretion, through the gas (if any) in small accreted halos and directly in the form of diffuse gas (inflows); and iii) via mass losses from the galaxies in the own halo (outflows). Halos themselves may lose the hot gas if it becomes unbound. AMIGA accurately takes into account all these mass exchanges between the different gaseous phases within ionized regions.
In neutral regions IGM remains unpolluted except for possible recombination periods after complete reionization. Therefore, Pop III stars with metallicity below some threshold value Z c for atomic cooling, in the range between 10 −5 and 10 −3 solar units, can only form through H 2 -molecular cooling in such neutral, pristine regions, which are immediately photo-dissociated and photo-ionized around those stars.
The evolving extent of the different IGM phases is described by means of the volume filling factors of singly and doubly ionized regions, Q HII ≡ n HII / n H and Q HeIII ≡ n HeIII / n He , which result from solving the improved differential equations where a dot means time-derivative, with trivial initial conditions in the dark ages. In equation (1) t is the cosmic time, subscripts S, SI, and SII mean either H, H I, and H II or He, He II, and He III, angle brackets stand for the average over the region indicated by the subscript (or over the whole IGM in the absence of subscript), n is the comoving IGM baryon density, a is the cosmic scale factor, and µ e is the electronic contribution to the mean molecular weight, i.e. (µ e ) −1 i is equal (neglecting metals) to zero, X + Y /4, or X + Y /2, in neutral, singly, or doubly ionized regions, respectively, X and Y being the usual hydrogen and helium mass fractions.Ṅ SII is the comoving metagalactic ionizing emissivity due to both luminous sources and recombinations, calculated according to Meiksin (2009), 2 α SI is the temperaturedependent recombination coefficient to the SI species, and C ≡ n 2 HII / n HII 2 , is the so-called clumping factor giving the ratio between the real recombination factor and that in an ideal homogeneous universe at a given z (see below). The term in claudators on the right of equation (1) takes into account inflows/outflows, and the volume average in region SII of the function of temperature f (T ) ≡ α SI (T ) is calculated using the Taylor expansion f (T i ) + (d 2 f /dT 2 i )σ 2 Ti /2 from the mean temperatures and variances, T i and σ 2 Ti , in phases i = I, II, and III for neutral, singly, and doubly ionized regions, respectively, solutions of the differential equations  d In equations (2) and (3) ε i , n i , and µ i are respectively the phase-i comoving energy density, comoving baryon density, and volume average mean molecular weight, i.e. µ −1 II + Y /(4X) Q HeIII (neglecting metals). The first term equal to 2 on the right of equation (2) arises from adiabatic cooling, while the second term includes the effects of Compton heating/cooling by CMB, Compton heating by X-rays, heating/cooling by ionization/recombination of the various hydrogen and helium species, cooling by collisional ionization and excitation, achievement of energy equipartition of newly ionized/recombined material, and inflows/outflows from halos. 3 Thus the evolution of luminous sources determines through equations (1)-(3) the evolution of the IGM, and conversely, the evolution of the IGM determines in the way shown next the evolution of luminous sources, as well as the halo MF in regions with different ionization state . Once the evolution of all these components has been determined we can calculate the CMB optical depth and the Compton distortion y-parameter In equations (4) and (5) H(z) is the Hubble parameter, σ T is the photon Thomson scattering cross-section, c is the speed of light, M hg (M, z) and T hg (M, z) are the average mass and temperature of the hot gas in halos with mass M , N i (M, z) ≡ N ion (M, z)Q i (z)/(Q II (z) + Q III (z)) is the halo comoving abundance in region i, being N ion (M, z) the halo MF in ionized regions, and M min,i (z) is the minimum halo mass in region i able to trap gas, all of them dependent on z.
2.2. Luminous Sources Halos grow through major mergers and smooth accretion.
AMIGA follows this growth in a well-contrasted analytic manner with no free parameters (Raig, González-Casado, & Salvador-Solé 2001; Salvador-Solé, Solanes, & Manrique 1998;Salvador-Solé et al. 2007). This allows us to accurately calculate, at any given moment, their abundance (Juan et al. 2014a,b) and inner structure and kinematics (Salvador-Solé, Viñas, & Manrique 2012;Salvador-Solé et al. 2012), which in turn sets the structure and temperature of their hot gas (Solanes et al. 2005). The amount and metallicity of such a hot gas and the properties of the central galaxy and its satellites are monitored in AMIGA from the individual halo aggregation history.
As long as the hot gas has a metallicity below Z c , molecular cooling takes place (see . When the temperature reaches a minimum value that cannot be surpassed, the cold gas accumulates in a central disk until the Bonnor-Ebert mass is reached. The cloud then collapses and fragments, giving rise to a small Pop III star cluster of about 1000 M ⊙ . 4 The initial mass function (IMF) of Pop III stars is poorly known, but at the present stage we do not need it. We only need the mass fractions f 2 and f 3 of Pop III stars with masses in the ranges 130 ≤ M ⋆ /M ⊙ < 260, and M ⋆ < 260 M ⊙ , respectively. Indeed, stars in the former range explode in pair-instability SN (PISN), release about half their mass in metals, and do not leave any black hole as remnant, while those in the most massive range essentially collapse into a black hole without producing metals (Heger & Woosley 2002). Thus, the yield of massive Pop III stars, p III , and the mass fraction of a Pop III star cluster ending up in a coalesced mini-MBH, β III , satisfy the relations 5 p III = 0.5 × f 2 (6) On the other hand, the H I-and He II-ionizing photon emissivities of massive Pop III stars depend on the mass fraction f 2 + f 3 in a well-known way (Schaerer 2002). 6 Lastly, Pop III stars with masses below 130 M ⊙ contributing to the mass fraction f 1 = 1 − f 2 − f 3 behave like ordinary Pop II stars, except for their low metallicity. Therefore, the feedback of Pop III stars is fully determined by the mass fractions f 2 and f 3 . When the metallicity of the hot gas exceeds Z c , it undergoes atomic cooling. The gas then contracts, keep-ing the initial angular momentum, and settles in a rotationally supported disk. The structure of disks is calculated according to Mo, Mao, & White (1998) from the specific angular momentum of the gas at the cooling radius, which equals that of dark matter as provided by numerical simulations. This leads to their effective (half-mass) radii r D and central surface density where M D is the disk mass. 7 If the disk is unstable, the cold gas is collected in a spheroid. Spheroids also form in galaxy mergers. Indeed, when a halo is captured by another one, it is truncated around the central galaxy which becomes a satellite orbiting within the new halo. Satellites undergo dynamical friction and orbital decay, being eventually captured by the more massive central galaxy or merging with it if the mass ratio between the two objects is greater than 1:3. Satellites also interact between themselves and with the hot intrahalo gas. However, except for ram-pressure stripping which is self-consistently calculated in AMIGA, these interactions (tidal stripping, gas starvation at the shock front of the hot gas, tidally induced disk-to-bulge mass transfer,...) are neglected in the present work as they should play no significant role at high-z when groups and clusters of galaxies are little developed and populated (but see Sec. 5.4).
Spheroids resulting from unstable disks or galaxy mergers have the Hernquist (1990) profile whose projection fits the r 1/4 law. Collisions between gas clouds yield the dissipative contraction of recently assembled spheroids with mass M B as stars form, according to the physically motivated differential equation for the effective radius r B ) from the initial value given by the r B − M B relation of non-contracted spheroids (see below) until the gas is exhausted or it reaches the density of molecular clouds (∼ 10 6 particles cm 3 ). In equation (8) M cg,B (t) and Z cg,B (t) are the mass and metallicity, respectively, of the cold gas in the spheroid at t,τ acc is the time elapsed from the formation of the spheroid to the quenching of star formation by the enlightened central AGN (see below), and ρ dis is a characteristic gas density for dissipative contraction of spheroids. When the cold metal-rich gas falls into a galactic component C, i.e. a bulge (C=B) or a disk (C=D), new Pop II stars form. The star formation rate (SFR) satisfies the usual Schmidt-Kennicutt laẇ where M cg,C (t) is the mass of cold gas available, τ dyn,C is the dynamical timescale at the half-mass radius of the galactic component, and α G is the star formation efficiency. The amount of interstellar gas heated by type II SNe in Pop II star formation episodes in the galactic component 7 Caution with the notation used in the present paper for the effective radius of a galaxy component. It coincides with that used in (Manrique & Salvador-Solé 2015) for its length scale.

C is
where V C is the circular velocity at r C , V hg is the typical thermal velocity of the hot gas in the halo where the heated gas is deposited, E SN = 10 51 erg is the typical energy liberated by one typical SN explosion, η SN = 0.014 M ⊙ −1 is the number of such explosions per unit stellar mass in a typical ∼ 0.20 Gyr duration starburst for the adopted IMF, and ǫ C is the SN heating efficiency of the component. Due to the different SFR, geometry and porosity of the gas in spheroids and disks, ǫ C is let to take different values in bulges (ǫ B ) and disks (ǫ D ). Disks smaller than the corresponding spheroids are however assumed to be oblate pseudo-bulges with ǫ B .
The amount of intrahalo gas heated through PISN during the explosion of very massive Pop III stars and leaving the halo is also calculated according to equation (10), but with ǫ C = ǫ B , E SN = E PISN = 10 53 erg, η SN = η PISN = 0.0015 M ⊙ −1 (Heger & Woosley 2002;Schaerer 2002), V C equal to the halo circular velocity at half-mass radius, and V hg = 0.
The instantaneous emission at all the relevant wavelengths of normal galaxies is calculated according to their individual stellar formation and metallicity histories using the stellar population synthesis model (SPSM) by Bruzual & Charlot (2003) for the adopted IMF (see below) and taking into account the typical yield of Pop II stars (Vicenzo et al. 2016) and all the mass and metallicity exchanges between the different galaxy and halo components.
In particular, as long as halos accrete, the metalenriched gas ejected from the central galaxy remains, due to viscosity, at the cooling radius where it is the next to cool. AMIGA assumes that, at the next major merger, all the hot gas of the progenitors is mixed up, except for their inner h mix mass fraction which retains all the reheated gas remaining in those progenitors. This means that metals ejected from central galaxies take at least two consecutive major mergers to get well mixed within the halo. As a consequence, the hot gas in the the inner part of the halo is always somewhat more metal-rich than in the outer part.
MBHs originate as remnants of the most massive Pop III stars that, as mentioned before, are supposed to coalesce in one mini-MBH per star cluster. When Pop III star clusters (or their remnants) are subsequently captured by normal galaxies and reach their spheroids, the corresponding mini-MBHs migrate by dynamical friction to their centers where they merge with the MBH lying there. MBHs at the center of spheroids progressively grow not only through mergers, but mainly through accretion of gas (and to a much lesser extent of stars). AMIGA does not assume any specific fraction of the gas fallen into spheroids to fuel the central MBH. Instead, it assumes that, each time new gas falls into a spheroid, part of it forms stars, and part fuels the MBH. The starburst lasts until it is quenched by the enlightened AGN that heats mechanically an amount of gas equal to and expels it back into the halo through super-winds. In equation (11) ǫ AGN is the so-called quasar-mode AGN heating efficiency, and L is the AGN bolometric luminosity L(t), averaged over the typical duration (τ acc = 0.15 Gyr; Kelly et al. 2010) of one such enlightening episodes. The bolometric light curve is modeled according to Hatziminaoglou et al. (2003) through the expression where V last is the Keplerian velocity at the last marginally stable orbit (at 9.2 times the Schwarschild radius) around the MBH, andṀ BH is the time-derivative of the MBH accretion curve, M BH (t), assumed with a universal dimensionless bell-shaped form, whose maximum at τ acc sets the quenching of the starburst going on in the spheroid. The fraction of H I-(and He II-) ionizing photons escaping from halos with a virial temperature below 10 4 K is self-consistently calculated by subtracting those photons ionizing the neutral intrahalo gas, taking into account recombinations. Above 10 4 K, a fixed value, f esc , is adopted, distinguishing between galaxies, f esc,G , and AGN, f esc,AGN .
The fraction of the SN explosion energy converted to X-ray photons through free-free emission of the SN remnant and inverse Compton scattering of CMB photons by relativistic electrons is about 1% (Oh & Haiman 2003), and the fraction of the AGN bolometric luminosity emitted in soft X-rays is 4% (Vasudevan & Fabian 2007). About half of those soft X-ray photons Compton heat the IGM in neutral regions due to the residual free electrons. The other half produce secondary ionizations and excitations (Oh & Haiman 2003), neglected in AMIGA.

Parameters
All quantities entering the previous equations are selfconsistently evolved in AMIGA from trivial initial conditions. Only a few of them must be supplied as external inputs.
Some of those external inputs correspond to reasonably well known functions of redshift, stellar mass, or galaxy mass. These are: -The clumping factor: Hydrodynamic simulations show that C(z) is equal to about 3 at z ∼ 6 (Finlator et al. 2012 and references therein), and evolves essentially as Strictly speaking, this fitting expression was derived by Finlator et al. (2012) from their simulations where reionization started at z ∼ 13, and ended up at z ion,H ∼ 8. In the real Universe, z ion,H is instead closer to z ∼ 6, and starts at a redshift z start to be determined. We thus adopt the following generalization of expression (13), with identical linear behavior but satisfying C(z start ) ≡ C F (13) = 1 and C(6) ≡ C F (8) = 3.
-The IMF of ordinary Pop II stars: Observations show that it can be approximated by a Salpeter IMF. In AMIGA we adopt the Salpeter slope, −2.35, for large masses up to 130 M ⊙ , and the slope −1 for small masses in the range 0.1 < M ⋆ /M ⊙ < 0.5. Such an IMF is consistent with the observed local one (Wilkins, Trentham, & Hopkins 2008), and similarly top-heavy as the Chabrier (2003) IMF.
-The r B −M B relation for non dissipatively contracted spheroids: Assuming for simplicity that it is a power-law, r B = AM γ B , and taking into account that spheroids at the two mass ends of the observed relations for nearby objects (Shen et al. 2003(Shen et al. , 2007Phillips 2003) have not suffered dissipative contraction, 8 one is led to A ≈ 0.049 Kpc M ⊙ −γ , and γ ≈ 0.20. Even though there is some uncertainty in the preceding expressions, the results of the present study have been checked to be very robust against reasonable variations in the coefficients.
But most of the external inputs are constant or can be considered as such in a first approximation. 9 Many of them are either well determined, as in the case of the threshold mass ratio for galaxy captures with destruction, the yield of Pop II stars, and the energy and frequency of SNe, or they play an insignificant role in the results, as the fraction of the energies released by SNe and AGN that are converted into X-rays, so they can safely be taken with the fixed values specified in Section 2.2. Others can also be taken with any reasonable fixed value because their uncertainty is absorbed in other parameters. This is the case of the factors 0.5 and 1.0 in equations (6)- (7), the constantτ acc in equation (8), and the product V 2 last τ acc times the factor arising from the universal dimensionless MBH accretion curve in equation (11).
However, all the remaining external inputs must be left free since they are poorly determined, and have a significant effect on the results. These are thus the real parameters of the model to be adjusted through the fit to the observational data. Their complete list is the following: Z c : threshold metallicity for atomic cooling, f 2 : mass fraction in intermediate Pop III stars, f 3 : mass fraction in massive Pop III stars, ρ dis : characteristic density for dissipative contraction, α G : Pop I & II star formation efficiency, ǫ B : SN heating efficiency in spheroids, ǫ D : SN heating efficiency in disks, ǫ AGN : AGN quasar mode heating efficiency, h mix : mixing hot gas mass fraction, f esc,G : escape fraction from galaxies, and f esc,AGN : escape fraction from AGN. As we will see in Section 4, there is some degeneracy between these 11 parameters. In addition, the unknown Pop III star IMF must satisfy some conditions that have not yet been enforced. As a consequence we will end up with only 9 degrees of freedom. Although this may still seem quite a large number, it is similar to that found in all recent analytic models of the EoR 10 that find 8 The most massive ellipticals have formed through dry major mergers, and spheroids with low enough mass reach the maximum density of molecular clouds. 9 For instance, the escape fractions from galaxies and AGN could vary with z. However, neither from a theoretical nor an observational point of view is that dependence clear. Thus we assume them with fixed typical values, as usual.
10 In all of them the clumping factor and the IMF of ordinary z ion,H = 6 from the fit to the observed SFR density of galaxies at z > 6 ( Robertson et al. 2015;Ishigaki et al. 2015;Bouwens et al. 2015). We remark, however, that our model is expected to recover not only the observed redshift of full hydrogen ionization, but the whole evolution of the Universe from the fit to all the observed cosmic histories at z > 2. Such an achievement with so few parameters is possible thanks to the fact that our model is self-consistent so that many quantities can be calculated and need not to be considered as free parameters.

THE OBSERVED HIGH-Z UNIVERSE
All the observational data currently available on the evolution of the Universe at z > 2 (see Fig. 1) refer to: (a) the cold gas mass density history (CGH), (b) the stellar mass density history (STH), (c) the MBH mass density history (MBHH), (d) the hot gas metallicity history (HGMH), (e) the cold gas metallicity history (CGMH), (f) the stellar metallicity history (STMH), (g) the IGM metallicity history (IGMMH), (h) the galaxy morphology history (GAMH), (i) the galaxy size history (GASH), (j) the SFR density history (SFH), (k) the H I-ionizing emissivity history (IEH), and (l) the IGM temperature history (IGMTH). Taking into account that there are two independent IEHs, one for galaxies and the other for AGN, this represents 13 independent global, i.e. averaged or integrated, cosmic histories. 11 In addition, the following differential properties are also available at a few discrete redshifts (see e.g. Figs. 8 and 11): (m) the galaxy stellar MFs (or UV LFs), and (n) the MBH MFs (or AGN optical and X-ray LFs). The mass density and emissivity estimates at different redshifts are derived from these MFs or LFs at the corresponding z, extrapolated beyond the observational limits. 12 Therefore, these MFs or LFs harbor more information than the corresponding global properties above. However, they are harder to handle. Thus we will concentrate in trying to fit the former cosmic histories, and then check whether the galaxy and MBH MFs at a few representative redshifts are also well recovered.
We briefly discuss next the techniques, approximations, and models used to infer all those data. All quantities given below refer to the Salpeterlike IMF adopted in AMIGA (applying the conversion factor from a strict Salpeter IMF provided by Wilkins, Trentham, & Hopkins 2008) and to a Hubble constant of 67.3 km s −1 Mpc −1 .
(a) Cold gas mass densities are computed from the distribution function of H I column densities per unit redshift in damped Lyα systems (DLAs) in the spectra of quasars up to z ion,H , integrated from zero to infinity.
The estimates shown in panel (a) of stars are also fixed, as well as the uncertain z-dependent dust attenuation. 11 The STH is the time-integral of the SFH, but the two observables are inferred independently. In fact, as we will see, they turn out to be inconsistent with each other.
12 Given the asymptotic logarithmic slopes of those MFs or LFs, such extrapolations are not expected to have any dramatic effect (but see Sec. 5.3).  Péroux et al. (2003) and Prochaska, Herbert-Fort, & Wolfe (2005). Note that some amount of H 2 (and other molecules), not accounted for in these estimates, must also be present. However, the real gas densities are not expected to be substantially greater.
(b) Stellar mass densities are inferred by fitting the SED of high-z galaxies, obtained from several optical and infrared rest-frame bands, to template SEDs generated by SPSMs assuming an IMF and a dust absorption correction. 13 The resulting galaxy stellar MFs are extrapolated down to 10 8 M ⊙ .
In panel (b) we show one illustrative sample of datasets obtained by several authors (Reddy & Steidel 2009;Stark et al. 2009;González et al. 2010;Labbé et al. 2010a,b;González et al. 2011;Mortlock et al. 2011a;Caputi et al. 2011). The large scatter in the data is mostly due to the very uncertain correction for dust attenuation.
(c) MBH masses at high redshifts are estimated by the method of reverberation mapping that assumes the broad line emitting clouds around AGN are virialized. The largest uncertainty in those masses arises from the unknown geometry and inclination of the system, and is encapsulated in a virial factor (f σ or f FWHM , depending on whether the potential well is traced through the cloud velocity dispersion or the emission-line width, respectively) of order one. The MBH mass density estimates at z < 5 and z = 6.1 plotted as open circles in panel (c) have been computed using the MBH MFs derived by Kelly et al. (2010) and Willott et al. (2010b), respectively, from AGN LFs in rest-frame UV. These MFs are, however, greatly affected by AGN obscuration and inactivity, and are substantially lower than the filled circle at z = 6.5 inferred by Treister et al. (2011) from the extrapolated AGN LF in X-rays, less affected by these effects. 14 On the contrary, this latter estimate is fairly aligned with the filled circles at the previous redshifts, derived assuming a constant MBH-to-spheroid mass ratio, µ, equal to 8 × 10 −3 (Kormendy & Ho 2013, taking into account the change in IMF), 15 and the upper envelope of the stellar mass densities given in item (b). 16 (d) Absorption lines of metal ions, usually C IV, are observed outside the Lyα forest in the spectra of distant quasars. The absorbers have column densities N according to a distribution function at any given  (i) median effective radii of spheroid-dominated galaxies (Buitrago et al. 2008); (j) SFR densities (Reddy & Steidel 2009;Lorenzoni et al. 2011;Bouwens et al. 2012b;Oesch et al. 2013;Coe et al. 2013;Ellis et al. 2013); (k) total ionizing emissivities (Becker & Bolton 2013) and the contribution from AGN (Cowie et al. 2009); and (l) IGM temperatures (Lidz et al. 2010;Bolton et al. 2010Bolton et al. , 2012. Error bars are 1σ deviations, except for panel (g) where they give upper and lower limits. The density parameters are for the current value of the critical cosmic density, ρc(z = 0), and the scaling factors in the right column are: r 0 = re(z = 0),ρ 0 = 1 M ⊙ yr −1 Mpc −3 ,Ṅ 0 = 10 51 photons s −1 Mpc −3 , and T 0 = 10 4 K. Other cosmic histories with no data also included for the densities and metallicities in the left and middle columns to cover all baryonic components: (x) hot gas mass densities, (y) IGM mass densities, and (z) metallicities of matter fallen into MBHs.
widest column density range available (10 12 −10 15 cm −2 ), converted to C abundances for a [C IV/C] abundance ratio of 0.05 appropriate to high column densities (Simcoe 2006), and then to metallicities assuming the C mass fraction of ejecta from ordinary Pop I & II stars given by (Ryan-Weber et al. 2009). Indeed, the absorbers are believed to be neutral polluted gas regions around intervening halos, with density contrasts spanning from 10 to 100 in redshifts going from 6 to 2, respectively. Consequently, these metallicities should coincide with that of the hot gas in low-mass halos having accreted gas for the first time (i.e. with masses between the minimum trapping mass at that z and, say, 1 dex larger).
(e) Gas-phase metallicities of high-z star-forming galaxies are derived from emission lines in their composite spectra. The estimates obtained for different metal ions are calibrated against each other, and converted to the [O/H] abundance ratio, then to the global metallicity, using the relation provided by Sommariva et al. (2012). In panel (e) we plot typical values obtained for small samples of galaxies with stellar masses above ∼ 5 × 10 9 M ⊙ at z around 3.5 and 2.2, respectively inferred by Maiolino et al. (2008) and Cullen et al. (2014), and several individual values for a few galaxies with similar masses at the remaining redshifts inferred by Sommariva et al. (2012).
(f) Stellar metallicities of high-z star-forming galaxies are derived similarly to the gas-phase metallicities in item (e) but from ion absorption lines. In panel (f) we plot the values obtained for a small sample of galaxies with stellar masses above ∼ 5 × 10 9 M ⊙ at z around 2.0 inferred by Halliday et al. (2008) and several individual values corresponding to a few galaxies with similar masses inferred by Sommariva et al. (2012).
(g) There are no direct measures of the average IGM metallicities. However, we can put in panel (g) some upper and lower limits. Indeed, as mentioned in item (d), metal absorption lines outside the Lyα forest in the spectra of distant quasars would be produced by the polluted IGM around halos. Since the metallicity of those systems decreases with increasing z due to the increasing effect of galaxy ejecta, their values for the most distant systems (5 ≤ z < 6) are the best upper limits to the average "primordial" (i.e. enriched by Pop III stars only) IGM metallicity we can have. For the present purposes, these metallicities are derived assuming a [C IV/C] abundance ratio of 0.5 (Ryan-Weber et al. 2009), appropriate to low column densities, and a carbon mass fraction of Pop III stars (model C in Schaerer 2002) as it corresponds to IGM polluted by Pop III stars only. On the other hand, the IGM metallicity can become lower than Z c at low-z due to the lower metallicity-to-ionizing photon contribution of normal galaxies into the IGM compared to that of Pop III stars that dominate at high-z. A lower limit of 0.5 × 10 −4 Z ⊙ is thus an educated guess.
(h) Morphological fractions are inferred from fitting a Sérsic (1968) law to galaxy surface-brightness profiles using circular apertures. This yields the Sérsic index, n, and the effective (half-light) radius, r e , of the galaxy. Alternatively one can decompose the surfacebrightness profile in a bulge and a disk. The fractions of spheroid-dominated star-forming galaxies with stellar masses ∼ 10 11 M ⊙ plotted in panel (h) were obtained by Buitrago et al. (2008) for a threshold Sérsic index value of n > 2, and by , at z = 2.5, for starforming galaxies with bulge-to-total light ratios greater than 0.5. Both methods yield similar results.
(i) Galaxy sizes depend on the galaxy morphology. The median effective radii of spheroid-dominated (n > 2) and disk-dominated (n ≤ 2) objects scaled to the corresponding local values (Shen et al. 2003) for galaxies with stellar masses ∼ 10 11 M ⊙ plotted in panel (i) as filled circles and open circles, respectively, were obtained by Buitrago et al. (2008) for the same galaxy sample as used to draw the morphological fraction in item (h). The estimates inferred by  directly using the effective radii of bulges and disks yield substantially different results and have not been considered. On the other hand, we will concentrate, from now on, in the effective radii of spheroid-dominated galaxies, as those of disk-dominated galaxies depend on the relative spatial distribution of stars and gas in the disk which is not modeled in AMIGA.
(j) SFR densities,ρ * , are derived from the luminosity of observed galaxies in broadband or narrowband filters that traces active star formation, correcting for incompleteness from the extrapolated LFs and for dust absorption, applying IMF-dependent SFR calibrations. In panel (j) we plot the estimates inferred by Reddy & Steidel (2009) (k) Similarly to the escape fraction of ionizing (or Lyman-continuum, LyC) photons, which has two different versions, one for galaxies and the other for AGN, we will distinguish between the H I-ionizing emissivity from normal galaxies and AGN. At high-z, the LyC emissivity of normal galaxies is usually estimated from their SFR densities assuming a given escape fraction of ionizing photons Richard et al. 2006;Stark et al. 2007). These estimates are not considered here as they are equivalent to the SFR densities plotted in panel (h). An alternative method (Becker & Bolton 2013), independent from the SFR densities consists of using the Lyα effective optical depth inferred from the comparison of the observed Lyα forest with synthetic data obtained from hydrodynamic simulations, modeling photon mean free paths for the assumed ionizing emissivities. These are the estimates plotted as filled circles in panel (k). Regarding to LyC emissivity from AGN, this is estimated from UV observations of a large sample of X-ray AGN. The results derived by Cowie et al. (2009) are plotted in the same panel as empty circles.
(l) The temperature of the singly ionized IGM is usually derived from the small-scale structure of the Lyα forest, using hydrodynamic simulations and galaxy and quasar emission models. This is the case for the data at z < 5 derived by Lidz et al. (2010) plotted in panel (l). We do not include the estimates inferred by Becker et al. (2011) as they greatly depend on the poorly known density-temperature relation. Unfortunately, all those measures are extremely challenging for z approaching 6. An alternative method that sidesteps such a difficulty consists of probing the IGM temperature within a proper distance of ∼ 5 Mpc of quasars by combining the cumulative probability distribution of Doppler broadening and synthetic Lyα spectra. The measurements plotted at z ∼ 6 were obtained in this way by Bolton et al. (2010) and Bolton et al. (2012), using one and 7 quasars, respectively.
(m) Galaxy stellar MFs are obtained from the galaxy LFs in rest-frame UV. The conversion involves the assumptions and procedures mentioned in item (b). For consistency, we adopt the MFs corrected for absorption inferred by González et al. (2011) at z = 3.8 and z = 6.8, which were also used to derive the stellar mass densities plotted in panel (b) at those redshifts that typically bracket the range of interest.
(n) Similar comments hold for the MBH MFs. Their derivation from optical AGN LFs involves the assumptions and procedures mentioned in item (c). For consistency, we adopt the MFs inferred by Vestergaard et al. (2008) and Willott et al. (2010b) at z = 3.3 and z = 6.1, respectively, which were also used to infer the MBH mass densities plotted in panel (c) at those redshifts that bracket the range of interest.

FITTING STRATEGY
The goal is to adjust the 11 free parameters mentioned in Section 2 through the fit to the 13 independent cosmic histories described in Section 3, enforcing the constraints on z ion,H and Q HII (z) mentioned in Section 1 as priors.
Scanning a whole 11-D parameter space is not an easy task. Fortunately, we can take advantage of some particularities of the problem that greatly simplify the fitting procedure.
First, the characteristic dissipation density, ρ dis , is fully constrained by the data, regardless of the value of any other parameter. Indeed, equation (8) recovers the observed effective radius of massive spheroids at z = 2−3 formed in major mergers of two galaxies with the typical gas content and metallicity (see panels (i), (a), and (e) of Fig. 1) provided only log(ρ dis M ⊙ −1 Kpc 3 ) = 6.0 ± 0.6. Second, the properties of normal galaxies and MBHs decouple from those of Pop III stars . Indeed, massive Pop III stars lose metals through PISN-driven winds into the ionized bubbles they form around. When this polluted gas falls into halos and cools, normal galaxies form. In principle, the properties of those galaxies should thus depend on the metallicity of the ionized IGM set by Pop III stars through the values of Z c , f 2 , and f 3 . But normal galaxies eject such large amounts of metals into the hot gas in halos through type II SN-and AGN-driven winds that the hot gas rapidly loses the memory of its initial metallicity, and the properties of galaxies rapidly become independent of Z c , f 2 , and f 3 . Similarly, MBHs are seeded by the remnants of massive Pop III stars. But the dramatic growth of MBHs within spheroids, regulated by the amount of gas reaching them, causes MBHs to rapidly lose the memory of their seeds, and the properties MBHs also become independent of Z c , f 2 , and f 3 . We can thus concentrate in adjusting, in a first step, the parameters referring to Pop III stars and, afterwards, those referring to normal galaxies and MBHs. Notice that, even though the properties of the two kinds of sources are independent of each other, we cannot exchange the order of those fits. Some of the observables involving the parameters of the latter set refer to mean mass densities per unit volume, not per unit ionized volume, or to metallicities averaged over all regions, not over ionized regions. Consequently, they involve not only the properties of normal galaxies and MBHs, but also the volume filling factor of ionized regions, which depends on the mass fraction, f 2 + f 3 , of massive Pop III stars. Nonetheless, once the parameters f 2 and f 3 are fixed, those observables will depend only on the parameters referring to normal galaxies and MBHs.
Third, most of the parameters in these two sets can be adjusted sequentially.
For a given value of Z c , the ionization and photoheating of the IGM at high-z depend only on the mass fraction f 2 + f 3 of massive Pop III stars, while the mass of metals ejected in ionized regions depends only on the mass fraction f 2 . Thus it should be possible to adjust f 2 + f 3 by fitting the IGMTH at high-z, and then adjust f 2 by fitting the IGMMH also at high-z.
On the other hand, for a given couple of ǫ B and ǫ D values, the amount of gas in the disk or spheroid of a galaxy at any moment is equal to the amount of gas fallen into it, which depends on cosmology (through the halo characteristics, the hot gas cooling rate, and the galaxy merger rate), 17 minus the amount of gas used to form stars and ejected from the component (the latter depending on the former). 18 The inductive reasoning thus implies that the stellar masses of normal galaxies depend only on α G . Similarly, the mass of gas fueling the central MBH of a galaxy depends on the gas mass collected in the spheroid, which together with the stellar mass sets the dynamics, and hence, the mass of newly formed stars and of gas heated by SNe, the mass of the MBH determining its accretion rate, and the mass of gas heated by the enlightened AGN. Therefore, the inductive reasoning implies that the MBHH depends on both α G and ǫ AGN . Lastly, the metallicity of the hot gas in the halo that determines, through the gas cooling and the metal-enriched gas ejected from galaxies, the cold gas and stellar metallicities of those objects depends on cosmology, on α G , ǫ AGN , and on h mix setting the mixing rate of the reheated gas in the halo. Finally, the emissivity of galaxies and AGN depends on all the previous parameters setting the whole evolution of galaxies and AGN and their respective intrinsic ionizing photon production rates and the values of f esc from the two kinds of sources. Consequently, it should be possible to adjust α G by fitting the SFH (or STH), then ǫ AGN by fitting the MBHH, then h mix by fitting the HGMH (or any of the CGMH and STMH), and lastly the f esc values of galaxies and AGN by fitting their corresponding IEHs with the suited priors on z ion,H and d z Q HII (hereafter d z means z-derivative).
Of course, this sequential fitting procedure should be carried out for every possible set of Z c , ǫ B , and ǫ D values in the corresponding 3-D space. But this is a minor difficulty compared to directly looking for acceptable solutions in the whole 11-D parameter space.
17 Cooling is much more sensitive to the density of the hot gas than to its metallicity, dependent on h mix . 18 The gas mass fueling the MBH is negligible compared to the mass converted into stars.
Any "acceptable solution" should thus fit with acceptable χ 2 values the IGMMH, the IGMTH, the SFH (and the STH), the MBHH, the HGMH (and the CGMH and STMH), and the IEHs for galaxies and AGN, with z ion,H = 6.0 +0.3 −0.5 and d z Q HII (z = 7) < 0. Moreover, since ǫ B and ǫ D govern the amount of cold gas remaining in disks and of stars forming in spheroids, the CGH and GAMH should hopefully also be fitted. And given the way ρ dis is determined, the GASH should too. In principle, the solution so obtained should also recover the observed galaxy stellar and MBH MFs and, provided Z c is not degenerate so that there is still one degree of freedom to play with, we may expect the solution to also fulfill the constraints on τ , and y. Thus the problem is not only well-posed, but in principle also slightly overdetermined.
Unfortunately, the IGMMH and IGMTH at high-z are not well-constrained. The only available data on the IGM temperature refer to z 6 where Pop III stars play no significant role, and there are just a few loose bounds for the IGM metallicity. Consequently, f 2 and f 3 cannot be determined in the previous simple manner. Moreover, Z c is degenerate with f 2 and f 3 . The good news is that this degeneracy reduces the degrees of freedom to 10, which facilitates the fit at the cost, of course, of making it even harder to find any acceptable solution.
Indeed, massive Pop III stars are the first source of ionizing photons, with emissivity equal to m p ξ III (f 2 + f 3 )Ṁ III , whereṀ III and ξ III are the formation rate and LyC photon production rate of Pop III stars, respectively. Those stars also pollute with metals the ionized bubbles at the rate 0.5f 2ṀIII . Therefore, the formation of normal galaxies is possible provided the metallicity of such bubbles, 0.5f 2 /[m p ξ III (f 2 + f 3 )], is greater than Z c . But this condition is necessary although not sufficient. When normal galaxies form, they keep on ionizing and photo-heating the surrounding IGM, while they lose very few metals into it (they rather enrich the hot intrahalo gas). Consequently, the metallicity of ionized bubbles decreases, and their temperature increases, so the formation of new galaxies may be quenched, and reionization may be aborted. Only if the initial IGM metallicity is sufficiently large will the formation of normal galaxies last long enough for reionization to proceed successfully.
According to the previous discussion, the evolution of the Universe for some values Z c , f 2 , and f 3 will thus be essentially the same (with Z IGM in units of Z c ) as for other values Z ′ c , f ′ 2 , and f ′ 3 , provided the IGM reionization and metal-enrichment histories driven by Pop III stars are the same, that is provided the two sets of values satisfy the relations This leads to the following simple fitting procedure. For some fiducial values of Z c , and f 2 /(f 2 + f 3 ) warranting that reionization will not be aborted, we scan all possible values of f 3 , ǫ B , and ǫ D , and for each set of these three values we adjust α G , ǫ AGN , h mix , and f esc in the sequential manner explained above, with the priors on z ion,H and d z Q HII . Then, given any acceptable solution, the relations (15) and (16) provide with any other combination of Z c , f 2 , and f 3 values leading to it. Lastly, we can check the values of τ and y as well as the right behavior of the galaxy stellar and MBH MFs.

Additional Constraints on the Pop III Star IMF
Were the Pop III star IMF known, we could readily calculate the mass fractions f 2 and f 3 , and, as a bonus, leave the above mentioned degeneracy. Unfortunately, such an IMF is poorly determined. But we can still enforce some conditions such an IMF must satisfy in order to further constrain the problem and reduce the degrees of freedom to 9 (8 if we discount ρ dis , directly determined by the data).
The IMF of Pop III stars must be top-heavier than the IMF of ordinary stars (Larson 1998 The rationale for a top-heavier Pop III IMF is that the fragmentation of protostellar clouds is favored by metals, absent in the pristine gas. 20 The upper mass of Pop I & II stars is poorly determined. An upper limit of 150 M ⊙ is sometimes considered, although no stars more massive than 130 M ⊙ have ever been observed (Figer 2005; but see Crowther et al. 2010).
21 This is the reason why in most studies dealing with the Pop III star IMF the value of 2.35 is adopted (e.g. Schaerer 2002).
M up III should be larger than 260 M ⊙ , otherwise Pop III stars would not seed MBHs, and smaller than 500 M ⊙ or at most than 1000 M ⊙ (Heger & Woosley 2002;Schaerer 2002).
To account for all these restrictions parameters f 2 and f 3 in the fitting procedure described above must be replaced by the equivalent ones (for a fixed α , and adjust in the sequential manner explained above the parameters α G , ǫ AGN , h mix , and the two f esc (see Fig. 2), with the priors on z ion,H and d z Q HII . Using the relations (15) and (16), we find all sets of Z c , α III , M up III , and M lo III values leading to every acceptable solution. Lastly, we check whether the constraints on τ and y are also fulfilled, and whether the observed galaxy stellar and MBH MFs are well-recovered.
The result will depend, of course, on the somewhat uncertain masses (130 M ⊙ and 260 M ⊙ ) delimiting the ranges where Pop III stars have distinct behaviors. But this is a small drawback compared to the gain of better constraining Z c and the Pop III star IMF.

RESULTS
Following this fitting procedure, we have looked for any possible solution, in the sense of being compatible with the currently available data on the high-z Universe and the usual reionization constraints, and checked whether it has single or multiple reionization episodes.
The idea of double reionization was advocated by several authors (Cen 2003;Chiu, Fan, & Ostriker 2003;Hui & Haiman 2003;Sokasian et al. 2004;Naselsky & Chiang 2004;Furlanetto & Loeb 2005) to reconcile the apparent inconsistency between the value of the redshift of complete hydrogen ionization obtained from the spectra of distant sources (z ion,H ∼ 6; Becker et al. 2001;Djorgovski et al. 2001) with that found from the CMB optical depth derived from the First-Year WMAP data (z ion,H > 15 as implied by τ = 0.17; Kogut et al. 2003). Subsequent CMB-based measurements decreased this latter estimate until the present value of z ion,H of about 10. In addition, it became clear that an instantaneous reionization event at z ∼ 10 was not necessarily inconsistent with reionization being completed by z ion,H ∼ 6. Consequently, the double reionization scenario lost its interest. However, the possibility that reionization is double still remains. As a matter of fact, it is quite natural to expect.
Indeed, Pop III star clusters have a typical mass of ∼ 10 3 M ⊙ , which is much lower than the mass ∼ 10 6 M ⊙ of the most abundant dwarf galaxies at high redshifts (see below). However, massive Pop III stars emit 40-80 times more ionizing photons than normal O, B stars, so both kinds of sources compete in the IGM reionization. What is the dominant contribution? Pop III stars are the first to form, so they necessarily drive the initial ionization phase. And, since neutral pristine regions progressively disappear, the final ionization phase will be driven by normal galaxies. The SFR of Pop III stars is fully determined by cosmology, and their more or less top-heavy IMF depends on the unknown value of M lo III (for the fiducial values of α III and M up III ). Whereas the IMF of Pop I & II stars, set by the physics of protostellar clouds, is well-determined, and their more or less intense SFR depends on the unknown value of α G . Consequently, the answer to the question above depends on the relative values of the two independent parameters M lo III and α G . For the α G value giving good fits to the SFH of galaxies, the larger the value of M lo III , the earlier Pop III stars will begin to ionize the IGM, and the more probable the double reionization scenario will be. Figure 3 shows the reionization histories of the only acceptable solutions we find. They are concentrated in two narrow sets, one with one single hydrogen ionization episode at z ∼ 6, and the other one with two ionization episodes at z ∼ 10 and z ∼ 6, separated by a short recombination period. The best solutions of the two kinds, i.e., those that provide the best fit to the data on the SFH, MBHH, HGMH, IEHs, GAMH, and GASH with the priors z ion,H = 6.0 and d z Q HII (z = 7) < 0, are hereafter referred to as "S" (for single reionization) and "D" (for double reionization).

Reionization History
In all the solutions, the reionization of He II also shows two distinct phases separated by a short recombination period. The first phase leads to a maximum in the volume fraction of He II-ionized regions, although without complete ionization, at z ∼ 10 in the solutions with double hydrogen reionization, and at z ∼ 6 − 7 in the solutions with single hydrogen reionization. The first phase is also driven by Pop III stars, but the second one, completed at z ion,He ∼ 2.8, is driven by AGN instead of normal galaxies. Thus, the He II reionization history confirms the idea that, when two independent ionizing sources with distinct timing are at play, it is natural to have a non-monotonous reionization. Table 1 lists the values of M lo III , ǫ B , ǫ D , α G , ǫ AGN , h mix , f esc,G , and f esc,AGN defining models S and D for the fiducial values of Z c , α III , and M up III given in Table 2. The errors in M lo III delimit the acceptable solutions around models S and D; those in ǫ B and ǫ D mark the deviations from the best values of these two parameters still leading to acceptable solutions; and the errors in the remaining parameters are 1σ deviations from their best values in models S and D, estimated from the χ 2 -test of their respective fits. In principle, the degeneracy between Z c and the Pop III star IMF could make the value of M lo III to vary from that quoted in Table 1. But, for the reasons explained next, there is actually no much room for the possible values of Z c , α III , and M up III , which should thus be close to those given in Table 2.
All the combinations of Z c , α III , M up III , and M lo III leading to any given solution can be obtained from one particular combination leading to that solution and equations (15) and (16) through the relations between the couples of M lo III and M up III , and f 2 and f 3 values for a given α III value. That degeneracy is shown in Figure 4. See also Figure 5 for the degeneracy in Z c , f 2 , and f 3 . Since the minimum possible value of f 2 /(f 2 + f 3 ) for reionization not to be In thick lines models S and D giving the best solution of the respective types. In thin lines the solutions bracketing the two sets. The error bars centered at z = 2.5 and z = 8.8 along Q SII = 1.0 give the estimated limits for the redshifts of complete helium and hydrogen ionizations, respectively (see Sec. 1); the vertical dotted black line marks the redshift z = 7 where Q HII is found to decrease with increasing z. (In the online journal, the solutions for single and double reionization in all the following figures in are in red and green, respectively.)

Table 1
Values of the adjusted free parameters (besides ρ dis = 10 6 M ⊙ Kpc −3 ) for the fiducial values of the degenerate ones. aborted depends on Z c , and, for a given α III , f 2 /(f 2 + f 3 ) is a function of M up III alone, there is also a one-toone correspondence between the maximum possible value of M up III and Z c . Notice that, for the reasons explained in Section 5.3, the maximum possible value of M up III for each Z c is smaller in single than in double reionization. Specifically, it is equal to the maximum M up III value found in double reionization times 0.85/0.21, where 0.21 and 0.85 are the minimum f 2 /(f 2 + f 3 ) values in single and double reionization, respectively.
As shown in Figure 4, the upper limit of Z c set by the lower limit of M up III equal to 260 M ⊙ is 10 −3.93 Z ⊙ in double reionization, and 10 −4.53 Z ⊙ in single reionization. Whereas the lower limit set by the upper limit of M up III equal to 500 M ⊙ (1000 M ⊙ ) is 10 −4.17 Z ⊙ (10 −4.30 Z ⊙ ) in double reionization, and 10 −4.77 Z ⊙ (10 −4.90 Z ⊙ ) in single reionization. Therefore, the central value and allowed range in double reionization are: log(Z c /Z ⊙ ) = −4.0 ± 0.1 Z ⊙ , whereas in single reion-ization they are: log(Z c /Z ⊙ ) = −4.6 ± 0.1. These are the fiducial values of Z c adopted for the fit. Note that, in the case of double reionization, the observed properties of the high-z Universe imply a value of Z c of 10 −4 Z ⊙ , which coincides with that preferred on pure theoretical grounds (Santoro & Shull 2005;Smith et al. 2009;Schneider & Omukai 2010). Instead, in the case of single reionization, they imply a value of Z c closer to the lower limit of 10 −5 Z ⊙ .
In Figure 4 we also see that the degeneracy track associated with any given solution depends slightly on α III , and so does the exact upper limit of Z c . But, for all reasonable values of α III , this dependence is only significant for unrealistic values of M up III close to 1000 M ⊙ . We can thus adopt α III = 2.35 without any loss of generality, even though the true value of this parameter could be somewhat smaller.
Of course, all M up III values smaller than the maximum possible value for a given Z c (i.e. on the left of the corresponding vertical line in Fig. 4) together with the M lo III value associated with it along the degeneracy track yield essentially the same solution. (In Fig. 5 all the couples f 2 and f 3 along the degeneracy track upwards that value of Z c also lead to the same solution.) But the maximum possible value of M up III in both single and double reionization, 300 M ⊙ , is already close to the lower limit of 260 M ⊙ , so the true value of M up III must be close, indeed, to the fiducial one of 300 M ⊙ .
Having fixed the values Z c , α III , and M up III , the central  and limiting values of M lo III given in Table 1 have been obtained as follows. In the single reionization case, the best solutions they define have the central and limiting values of z ion,H , equal to 6.0 +0.3 −0.5 , the condition d z Q HII < 0 at z = 7 being automatically fulfilled in this case. In the double reionization case, the best solutions they define have z ion,H = 6.0 +0.3 −0.5 , but, in addition, the minimum of Q HII (z) is at a redshift larger than, although as close as possible to z = 7. This warrants d z Q HII < 0 at z = 7 as wanted, together with the minimum possible redshift of first complete ionization, which is very convenient to have the smallest possible value of τ (see below).
The values of ǫ B and ǫ D are well constrained by the observed GAMH and GASH, as expected. In fact, those two properties appear to be extremely sensitive to small changes in the values of these (as well as other) parameters. Their values close to one and zero, respectively, in the two models agree with the results of analytic and hydrodynamic studies on the effects of SN heating (e.g. Dekel & Silk 1986;Mac Low & Ferrara 1999;Strickland & Stevens 2000), and they are also consistent with some properties of the nearby Universe such as the high abundance of α-elements in the hot gas of halos or the correlations between the metallicity of the hot gas and the abundances of early and late type galaxies in galaxy clusters indicative that most gas outflows take place from spheroids (e.g. Silk 1997).
The values of α G and ǫ AGN are also well constrained by the observed SFH and MBHH. The larger value of α G in model S compensates the lower volume filling factor of ionized regions at high redshifts in this model compared to model D. Indeed, the SFH is averaged over the whole space, while stars lie in ionized regions only. This then leads to a smaller value of ǫ AGN in model S so as to recover the right MBH-to-spheroid mass ratio despite the larger amount of stars forming in individual spheroids.
The value of h mix is the least constrained. The reason for this is that the data on HGMH show a large scatter, and, even though the data on the CGMH and STMH are better determined, these cosmic histories are quite insensitive to h mix . But this result, far from being negative for our purposes, is heartwarming as the mixing of metals in the hot gas is hard to model, and the hot gas metallicity estimates are not very reliable (see Sec. 3).
The values of f esc,G and f esc,AGN are also well constrained by the observed IEHs from galaxies and AGN and, more importantly, by the priors z ion,H = 6.0 +0.3 −0.5 and z ion,He ∼ 2.8. Notice that the uncertainty of about ±0.8 in z ion,He (Ricotti et al. 2000;Schaye et al. 2000;Dixon & Furlanetto 2009;McQuinn 2009;Becker et al. 2011) has not been taken into account since any change of that order in z ion,He has no significant effect in the remaining cosmic properties. The values of f esc,AGN are one order of magnitude smaller than those of f esc,G , which is consistent with the extra absorption of ionizing photons produced in the circumnuclear region of galaxies. Regarding the values of f esc,G , they agree with the observational estimate (spanning from 0.03 to 0.07) at redshifts out to z ∼ 5 (see e.g. Wyithe et al. 2010 and references therein). But these empirical estimates are quite uncertain, while all recent theoretical studies of reionization based on the observed SFH find instead f esc,G ≥ 0.20 (Robertson et al. 2015;Ishigaki et al. 2015;Figure 6. Effective LyC photon production rates predicted in models S and D for galaxies with any mass and metallicity (solid lines) and for galaxies with stellar masses larger than ∼ 5 × 10 9 M ⊙ and near solar metallicities only (dashed lines; the predictions from the two models superpose). The dashed line coincides, except for a factor 1.61 due to the different IMF used, with the LyC photon production rate predicted in the reionization models by Robertson et al. (2015), Ishigaki et al. (2015), and Bouwens et al. (2015). (A color version of this figure is available in the online journal.) Bouwens et al. 2015). About half this difference arises from the fact that those studies neglect the effect of Pop III stars (AGN have a negligible contribution to the H I-ionizing emissivity at z > 6), and the other half from the notably larger effective value of the LyC photon production rate from galaxies, ξ ion , defined through the re-lationṄ HII = f esc,G ξ ionρ * , found in our self-consistent modeling compared to that found in those studies (see Fig. 6). As we will show in Section 5.3, that increased ξ ion value is due to a notably larger (by a factor ∼ 2) contribution of low-mass galaxies with much lower metallicities (more than one order of magnitude lower than near-solar).

Fits to the Observed Global Properties
In Figure 7 we show the fits to the global data that result in the best models S and D. As it can be seen, both solutions recover satisfactorily the observed evolution of all cosmic properties. The only exception is the STH (panel (b)) and, to a much lesser extent, the CGH (panel (a)). Since the STH is the time-integral of the SFH, and this latter curve is very well fitted (see panel (j)), the poor fit to the STH simply reflects the well-known inconsistency between the datasets of the two cosmic histories. We will comeback to this point in Section 5.3. The possible origin of the slight excess in the predicted cold gas mass density at z < 3 will also be discussed in Section 5.4.
The fits to the data are particularly good in the case of model D. The predicted SFH almost fully matches the data, while in model S it only does marginally (see panel 0.26 ± 0.01 * This empirical value of τ combines those inferred from the WMAP 9-years and the Planck three-years data. (j)). This is reflected in the respective χ 2 values of 1.07 and 2.65. 22 Likewise, the predicted CGH in model D is almost as shallow as the observed curve in the redshift range covered by the data, while that in model S is slightly steeper (see panel (b)); the respective χ 2 values are 1.90 and 2.35. The predicted HGMH, CGMH, and STMH curves for the metallicities averaged over all halos and galaxies are also substantially different in the two models. However, when the averages are restricted to the mass ranges included in the observations, the two models agree with the data (see panels (e)-(f)). Similarly, the predicted GAMH and GASH are somewhat different at high-z in the two models, but they both match the data at z ∼ 2−3 (see panels (h) and (i)). Lastly, the predicted IGMMH is larger in model D than in model S (see panel (g)). However, the empirical bounds on this quantity are little reliable (see below). A full comprehensive picture of the interconnected behavior of all these cosmic histories will be given in Sections 5. 3 and 5.4. And what about the observational constraints on reionization? In Table 3 we give the values predicted in the best models S and D. The error bars mark those predicted in the models bracketing the two sets of acceptable solutions around them. The values of z ion,H and d z Q HII satisfy, of course, the observed constraints because they have been enforced as priors. But the values of τ and y are checked for the first time. The values of the Compton distortion y-parameter predicted in models S and D pass the test without trouble: they are one order of magnitude lower than the observed upper limit. The predicted values of τ are also consistent with the empirical estimates drawn from the WMAP9 and Planck3 data at the 1σ level. However, while the value predicted in model S still agrees at the 1σ level with the estimate of τ = 0.058 ± 0.012 recently reported by the Planck team, that in model D only does at the 3σ level. We remind that the way the condition d z Q HII (z = 7) < 0 was enforced in this model warrants that the redshift of first complete reionization is as small as possible, so the only way to further decrease the predicted value of τ is to relax that condition. This possibility cannot be discarded, indeed: the typical outflow velocity of gas clouds in LAEs could decrease with increasing z so that the absence of LAEs at z > 7 would not imply that Q HII (z) is decreasing at z = 7. 23 In any event, there is also the caveat that the value of τ inferred from the CMB anisotropies under the instantaneous reionization approximation is quite un- Figure 7. Evolution of galaxies and IGM predicted in models S and D, compared to the observational data shown in Figure 1. We plot the contribution from normal galaxies (solid lines) and normal galaxies plus Pop III star clusters (dotted lines); in panel (k) we also plot the contribution of AGN (dashed lines). In panels (d), (e), and (f), we plot the average metallicities in halos and galaxies within the mass ranges covered by the data (thin solid lines), as well as the average metallicities over all halos and galaxies (thick solid lines) with no corresponding data. In panels (h) and (i) predictions are only shown at the redshifts where galaxies with masses around 10 11 M ⊙ are abundant enough. In panel (l) we show the average IGM temperatures in singly ionized (solid lines), doubly ionized (dashed lines), and neutral (dotted lines) regions, the latter shifted a factor 500 upwards at z higher than the redshift of first (or unique) complete ionization. Note that the sum at any z of the densities in the left-hand panels gives the total baryon mass density at that z (scaled to the constant critical density at z = 0), and the sum at any z of the product of the metallicities in the middle panels times the corresponding densities in the left panels gives the cumulative mass density of metals produced until that z. (A color version of this figure is available in the online journal.) certain.

Checking the Differential Properties
But the fact that models S and D fit all the observed global properties, in particular the SFH and MBHH, does not necessarily mean that they also recover the observed galaxy stellar or MBH MFs. In the present Section we check that, describing in detail the formation and properties of galaxies and MBHs of different masses.
At high-z halos are little massive and their concentration is high, so the hot gas cools very efficiently. But, for this to happen, halos must first trap gas, that is their virial temperature must be higher than the IGM temperature, T IGM . The minimum mass of halos trapping gas at z is where we have taken into account that the inner mean density of halos at z is a factor ∆(z) (close to 178 at high-z) times the mean cosmic densityρ(z). Therefore, new luminous objects (Pop III star clusters or galaxies) will form as long as new halos will trap gas (with undercritical or overcritical metallicity, respectively), that is as long as M trap (z) will decrease or increase less rapidly than the always increasing typical halo mass of collapse at z, M c (z), equal to the mass for which the rms density perturbation equals the critical overdensity for ellipsoidal collapse at z (e.g. Juan et al. 2014b). M c (z) is fully determined by cosmology, so the formation of new luminous objects depends only on the evolution of T IGM (z). The temperature of the neutral cosmic gas decouples from the CMB temperature at z ∼ 150, when Compton heating of residual electrons from CMB photons becomes negligible. Then, the cosmic neutral gas cools adiabatically as (1 + z) 2 until the first Pop III stars form. The X-ray background they produce (mostly through PISNe explosions) heats the neutral regions above the CMB temperature by Compton scattering from the diffuse residual electrons. However, Compton cooling from CMB photons prevents such a heating from being too marked, and T IGM (z) diminishes essentially as 1 + z until Compton cooling progressively declines by the dilution of CMB photons. Then, the temperature of neutral regions keeps on diminishing almost as (1 + z) 2 as X-ray heating is insufficient to balance adiabatic cooling (see Fig. 7, panel (l)).
Such an evolution of T IGM (z) in neutral pristine regions implies that M trap (z) decreases roughly as (1 + z) 3/2 all the time except for a short intermediate period where it is kept essentially constant. Meanwhile, M c (z) starts growing very rapidly, and progressively slows down from (1 + z) −2 above z = 20 to (1 + z) −1 below z = 6, being proportional to (1 + z) −3/2 at an intermediate redshift of about 10. As a consequence, Pop III star clusters are permanently forming in neutral pristine regions. The situation is different, however, in the neutral polluted nodes that develop in model D during recombination. There, the temperature is kept much more constant (see Fig 7, panel (l)), 24 causing M trap (z) to increase roughly as (1 + z) −3/2 , that is more rapidly than M c (z) since z ∼ 10. Consequently, no Pop III stars will form in this model after that redshift.
More specifically, M c (z) starts from a very tiny value, and does not reach what is believed to be the neutralino free streaming mass, 10 −6 M ⊙ , until z ∼ 23. This means that the only halos present at any higher redshift (necessarily more massive than the free streaming mass of DM particles) are much more massive than the typical collapse mass at that z. Thus their abundance is insignificant. In fact, the first mini-halo (of ∼ 5 × 10 5 M ⊙ ) able to trap gas and form a Pop III star cluster (of ∼ 2000 M ⊙ ) within the current horizon is found at z ∼ 48.
As the temperature of the neutral IGM decreases, the mass of mini-halos able to form Pop III star clusters diminishes, and, at z ∼ 45, it is low enough for such halos to lose their metal-enriched gas into ionized bubbles when the most massive Pop III stars explode in PISNe. Then the IGM metallicity in ionized regions rapidly increases above Z c . From that moment, when the ionized metal-enriched gas is trapped by halos, it rapidly undergoes atomic cooling, and gives rise to the formation of normal galaxies. But the IGM in ionized bubbles is also photo-heated to a much higher temperature than in neutral regions, so halos able to trap ionized gas are again very massive and rare. In fact, at those redshifts most galaxies form instead in slightly less massive halos (of 10 6 − 10 7 M ⊙ ) that retained part of the hot gas when Pop III stars exploded. Since the amount of gas remaining in such halos can be very small, most of those galaxies are "failed" in the sense that they have very low stellar masses (from 10 2 M ⊙ to 10 5 M ⊙ ).
As halos able to trap fresh metal-enriched ionized gas become increasingly abundant, more and more galaxies form with stellar masses above 10 5 M ⊙ . This leads to the development of a bimodal MF. But the formation of normal galaxies causes ionized regions to be photo-heated to an approximately constant temperature, so M trap (z) finally increases at the same rate, (1 + z) −3/2 , as M c (z) at z ∼ 10. Consequently, the formation of new galaxies with stellar masses above ∼ 10 6 M ⊙ (the lower mass limit is slowly increasing with decreasing z) is halted roughly at that redshift.
From that moment, the evolution of luminous sources differs in the two models. In model S, Pop III star clusters keep on forming in neutral pristine regions, and galaxies keep on growing, although not forming, in ionized ones. While, in model D, Pop III stars stop forming at z ∼ 10 when the first complete reionization takes place. (Remember that no Pop III stars form in the neutral polluted regions that develop during the subsequent recombination period.) There is thus a sudden lack of ionizing photons that could definitely interrupt the reionization process. Fortunately, recombination yields a big dip in the IGM temperature of ionized regions (see Fig. 7, panel (l)) that allows halos with masses as small as 10 8 M ⊙ to trap fresh ionized gas. As explained next, the metallicity of this gas is still overcritical, so an abundant amount of normal galaxies then form which partially reble reionization by looking at the IGMTH. Although there is no temperature estimates at z > 6, the spin temperature of hydrogen atoms is already coupled the neutral gas temperature through the Wouthuysen-Field effect, so such a difference between models S and D should leave an unequivocal signal in 21 cm line observations. place Pop III stars in the production of ionizing photons. This will allow reionization to recover by z ∼ 7.2.
Indeed, the IGM metallicity in ionized regions is initially kept roughly constant due to the fixed proportion of metals (∝ f 2 ) and ionizing photons (∝ f 2 + f 3 ) provided by Pop III stars. But, once normal galaxies begin to form, there is a larger contribution of ionizing photons than of metals (remember that metals ejected from galaxies mainly stay within halos), so the metallicity of ionized regions begins to decrease. This does not prevent, of course, the total average IGM metallicity from growing (Fig. 7, panel (g)) through the increasing weight of ionized regions. But the metallicity of ionized regions themselves can become undercritical. This is actually what happens in model S where the minimum value of f 2 /(f 2 + f 3 ) just needs to warrant the overcritical metallicity when ionized bubbles form. (The undercritical metallicity of these regions at z < 10 in model S has no effect in the reionization process because no galaxies can form anyway in this model since z ∼ 10.) The situation is different, however, in model D where the IGM temperature decreases since z ∼ 10. Provided the minimum value of f 2 /(f 2 + f 3 ) is large enough (the maximum value of M up III low enough) to warrant that the IGM metallicity in ionized regions is still overcritical at z ∼ 10, new galaxies will form, indeed, in this model from that redshift to about z ∼ 7 when the temperature of the ionized IGM increases again at a high enough rate for M trap (z) to grow more rapidly than M c (z). 25 The minimum mass of trapping halos in model D during that period is the largest ever been, and the metallicity of the initial hot gas is the lowest, so the cooling rate reaches very low values. As a consequence, a large fraction of these second generation normal galaxies have very small stellar masses, similar to those of the above mentioned failed galaxies.
Since z ∼ 10 in model S or since z ∼ 7 in model D, no new galaxies form. Galaxies just grow at the center of halos through the accretion of new cooled gas and satellite captures. Of course, some of them freeze out as they become satellites of other central galaxies when their host halos are captured by other halos. Consequently, at any z lower than 7, there are in both models galaxies with the whole spectrum of stellar masses. The only difference is that, while in model S all the very low mass satellites are failed galaxies, in model D many of them are second generation galaxies formed at z < 10. Notice that, even though halos where those two kinds of galaxies form have similar masses, M c (z) is much greater at z < 10 than at the formation of failed galaxies, so those second generation low-mass galaxies are more abundant than failed ones. Figure 8 shows the galaxy stellar MFs that result in models S and D at redshifts z = 3.8 and z = 6.8. To compare them with the empirical MFs inferred by González et al. (2011) at the same redshifts we need to account for the above mentioned inconsistency between the SFH and STH datasets. Indeed, at high-and low-z, the predicted STH curve, consistent with the empirical 25 In fact, since the maximum possible value of M up III is close to the lower limit of 260 M ⊙ , the average IGM metallicity at z ∼ 10 cannot be much larger than Zc either (see Fig. 7, panel(g)). Therefore, another characteristic of double reionization is that the IGM metallicity at 7 z 10 is close to Zc. SFH one, is close to the upper envelope of the stellar mass density estimates, and between z 6 and z ∼ 4 it is even higher (see Fig. 9). Wilkins, Trentham, & Hopkins (2008) put forward that such an inconsistency could be due to a z-varying IMF. However, the largest discrepancy is found at an intermediate redshift range, which would require the IMF to have a strange non-monotonous trend. A simpler explanation is that the correction for dust attenuation achieved in the galaxy UV luminosities is typically underestimated. Galaxies in the redshift range 4 z 6 would be substantially brighter than believed, the effect being less marked at higher and lower redshifts because of the lower metallicity of galaxies and the more accurate correction achieved, respectively. But, whatever the real cause of that inconsistency is, provided all the observed galaxies are similarly affected, the shape of the MFs derived from the UV LFs would be right except for an horizontal shift. Specifically, according to the stellar mass density estimated inferred by González et al. (2011) at redshifts 3.8 and 6.8 and the theoretical STH curve consistent with the empirical SFH one, the MFs used by González et al. (2011) to infer those estimates should be shifted in model D (S) by a factor of about 6.0 (7.9) and 1.3 (0.8), respectively. Taking into account that the extrapolations of those empirical MFs were underestimated by factors 1.6 (1.3) and 2.1 (1.4), respectively, the actual shifts to be applied at those redshifts are 9.5 (10.3) and 2.7 (1.2), respectively. As it can be seen in Figure 9, these shifts make the empirical MFs fully match indeed the predicted ones except for the incomplete low- mass bins (see below). Thus the predicted galaxy stellar MFs are in very good agreement with the observations. The fact that the MF predicted at z = 3.8 also matches the Schechter fit to the original (with no shift) galaxy stellar MF recently derived by Caputi et al. (2015) at 3 ≤ z ≤ 4 assuming a constant IMF also gives strong support to the interpretation that the inconsistency between the SFH and STH datasets is due to a deficient correction for dust attenuation. According to our results in model D (S), the UV luminosities of galaxies inferred by González et al. (2011) at z = 3.8 and z = 6.8 would be underestimated by 2.5 (2.7) mag and 1.0 (1.1) mag, respectively.
Regarding the behavior at lower-masses of the predicted MFs, we find the turnover at ∼ 3.5 dex from the knee in model S, and at ∼ 4.5 dex in model D. These results are consistent with the latest UV LFs of galaxies observed at z ∼ 3 (Livermore et al. 2016;Alavi et al. 2016), which find that the turnover is not yet detected at ∼ 6.5 mag (or ∼ 2.6 dex in mass) from the knee. We remark that all previous models of galaxy formation predicted the turnover at just ∼ 3.5 mag (∼ 1.4 dex in mass) from the knee (Jaacks et al. 2013;Kuhlen et al. 2013), which is contradicted by these observations. In this sense, our model is the first to predict MFs consistent with the observations. According to it, such LFs should still go ∼ 2.25 mag or ∼ 4.75 mag deeper in models S or D, respectively, for the turnover to be finally detected.
Model D also predicts the existence of a second bump at very low masses in the galaxy stellar MFs at z = 3.8 caused by the second generation low-mass galaxies that become satellites soon after their formation. This bump is far from the current observational limits of the galaxy stellar MFs at those redshifts. However, given the small mass of those satellites, a large fraction of them should still survive at z = 0. In model S, there is no such a bump, although the abundance of very low mass galaxies is also quite high. Thus a prediction of our model, regardless of the version, is that it should be possible to detect at z = 0 very old satellites with masses between 10 2 M ⊙ and 10 6 M ⊙ . We describe next the expected metallicity, morphology, and size of those predicted relics.
The hot gas metallicity increases in parallel to the cold gas and stellar metallicities (see Fig. 7,panel (d), (e) and (f)). The only effect that partly balances that increase is the accretion by halos of low-metallicity IGM in ionized regions. That increase is more marked in model S than in model D due to the steeper SFH caused by the larger α G value, and hence, the larger amount of stars that are permanently forming in galaxies. Another more subtle difference between the two models is the small dip in the hot gas metallicity in model D after z ∼ 10 due to the large amount of fresh gas that is being incorporated by halos during recombination owing to the dip in the temperature of the ionized IGM. This translates, of course, in less marked dips in the CGMH and STMH of model D. The stellar and cold gas metallicities of the second generation galaxies formed in this model during recombination, similar to those of the primordial failed galaxies formed in both models, are of about 10Z c , i.e. significantly lower than average. 26 On the contrary, the metallicities of galaxies more massive than 5 × 10 9 M ⊙ are, of course, larger than average (see also Fig. 7, panels (e) and (f)) and show a well-understood behavior with z: at the highest redshift where those galactic masses are assembled, such galaxies are extraordinarily massive, having been formed in the merger of evolved high metallicity galaxies; then the metallicities slowly decrease as the masses become increasingly ordinary and such galaxies form through other channels (including accretion of hot gas with lower metallicity); and, when the masses of the galaxies approach the average galaxy mass at z ∼ 4, they begin to increase just as the average metallicities do since z ∼ 6.
The morphology of galaxies follows from the following general trends. When a galaxy forms, the gas it collects comes from the very central part of the halo, so it has little angular momentum, and it is pilled up in a spheroid. Consequently, the failed galaxies in models S and D as well as the second generation low-mass galaxies in model D are spheroids. But, soon after their formation, the morphology of galaxies evolves at the center of halos. At very high-z, when cooling is very efficient, galaxies usually collect large amounts of cold gas that render disks unstable. In addition, halos are small, and galaxy mergers frequent, so, whenever disks are stable, they are short-lived and never reach large radii. But, as halos grow, cooling becomes less efficient and galaxy mergers less frequent, so disks grow larger and are more long-lived. As a consequence, the fraction of spheroiddominated galaxies slowly decreases with decreasing z. On the other hand, at any given z, the most massive galaxies tend to be spheroids arising from galaxy mergers of massive objects, while disks preferentially lie in intermediate mass galaxies. These trends explain the observed fraction of spheroid-dominated objects with masses above ∼ 10 11 M ⊙ at z ∼ 2 − 3. We remark, however, that about 25% of galaxies at those redshifts are red dead ellipticals with no star formation. These galaxies have been excluded from the results shown in panel (h) of Figure 7, as they would go undetected. Therefore, the real fraction of spheroid-dominated massive galaxies at those redshifts is ∼ 0.65 rather than ∼ 0.40.
The size of galaxies increases in parallel to their mass since the beginning. In fact, when they form as spheroids from a gas cloud with small angular momentum, even though the metallicity of the gas is low (about Z c ), they undergo a large dissipative contraction until they reach the density of molecular clouds. This leads to stellar systems with the typical density of globular clusters. The size of more evolved spheroids depends on the amount and metallicity of the gas they collect (see eq. [8]). The result is an evolution of the average equivalent radius (measured as the square root of the mass-weighted squared equivalent radii of the bulge and the disk, if any) of spheroid-dominated galaxies with masses above 10 11 M ⊙ that recovers the observations at z ∼ 2 − 3 (Fig. 7, panel (i)). Specifically, when those very massive galaxies appear for the first time in dry mergers of massive spheroids, their equivalent radius is even slightly larger than that of present day slightly contracted spheroids of the same mass; this radius begins to decrease as the mass becomes increasingly ordinary and those galaxies form from unstable disks and wet mergers; it then increases again up to a secondary maximum caused by the increasing lack of gas followed by the increase of the gas metallicity at z ∼ 5 (see Fig. 7, panel (e)); and it finally increases until z = 0 due to the definite small amount of gas present in galaxies.
Thus the properties of very low mass galaxies are essentially the same in both models. They are among the first objects with Pop II stars to form; they are spheroids with no gas; they are very compact, with the typical density of globular clusters; and their Pop II stars have very low metallicities (of about 10-100 times Z c , with Z c equal to ∼ 0.25×10 −4 Z ⊙ in model S, and ∼ 10 −4 Z ⊙ in model D). On the other hand, those of them that froze out as satellites soon after formation have very large baryon-to-DM mass ratios as their original DM halos are truncated by the potential well of the also quite dense capturing halo near the region where all the stars are strongly concentrated. And, the fraction of those satellites that will be captured long time after (because of the slight dynamic friction they suffer) by evolved, much less dense galaxies will survive inside them as small independent dynamical entities.
Clearly, all these properties greatly resemble those of present day globular clusters. Since such objects are often seen to freely orbit within halos, they should have formed outside galaxies and have a cosmological origin similar to them. It is thus not surprising that a selfconsistent model such as AMIGA that monitors the formation and evolution of luminous objects since the dark ages does not miss them. Notice that, if this is indeed the origin of globular clusters, the larger abundance of such  Figure 7 but extended down to z = 1. Also shown are the spheroid mass density history (thin dashed lines), and disk mass density history (thin dot-long dashed lines). Empty circles give the MBH mass density estimates obtained from rough (incomplete and uncorrected for obscuration and inactivity) MBH MFs. Brown dots at z = 3.3 and z = 6.1 correspond to the estimates derived by Kelly et al. (2010) and Willott et al. (2010b), respectively, from the MFs used in Figure 11. (A color version of this figure is available in the online journal.) objects compared to massive galaxies is more consistent with model D than with model S.
Let us now turn to the MFs of MBHs. MBHs grow since the remnant mini-MBHs of Pop III star clusters begin to develop at the center of galaxies (Fig. 7, panel (c)). There, they merge with the MBHs brought by accreted or merged galaxies having reached the center of the galaxy by dynamical friction, and accrete new material. Indeed, when some amount of gas reaches the spheroid it fuels the MBH. The associated AGN is then lit, and begins to reheat the circumnuclear gas, quenching the star formation that is going on in the spheroid. The more massive the MBH, the larger both the AGN bolometric luminosity and the fraction of gas reheated and expelled from the spheroid back into the halo. This leads to a self-regulated growth of MBHs and stellar spheroids, which tends to a constant MBHto-spheroid mass ratio as observed (Fig. 10). This thus gives strong support to the mechanism proposed by Di Matteo, Springel, & Herquist (2005) included in AMIGA for the coupled growth of MBHs and spheroids.
The parallel evolution of spheroids and MBHs goes together with that of disks. Of course, when gas replenishment begins to decline by z ∼ 7 (see Section 5.4), the disk and spheroid mass density histories become slightly shallower. (The change of slope is a little less marked in spheroids than in disks due to their unchanged behavior in major mergers.) But this has a minor effect in the parallel evolution of the spheroid and MBH stellar mass densities because mergers become at the same time increasingly dry. Indeed, in dry mergers spheroids grow through the addition of stars belonging to the merging galaxies, while MBHs grow by coalescence of the respective MBHs and, since the MBH-to-spheroid mass ratio is the same in all the progenitors, it remains unchanged in the final object. As a consequence, the MBH and spheroid mass ratio will be kept essentially unaltered until z = 0. Figure 11 shows the MBH MFs predicted at redshifts z = 3.3 and 6.1. Like in the case of the galaxy stellar MFs, to compare these MFs to the empirical ones derived at the same redshifts by Willott et al. (2010b) and Vestergaard et al. (2008) we must take into account any effect that could affect the latter. As mentioned, the observed AGN LFs are greatly affected by incompleteness, obscuration, and inactivity, which cause the MBH mass density estimates derived from them to be much lower than expected from the observed constant MBH-to-spheroid mass ratio. Incompleteness affects in an uneven way the different mass bins of the MBH MFs, but their large mass end should be complete. Obscuration is also larger for low-mass objects (Ueda et al. 2003), but the most massive bins should be similarly affected. Thus a vertical shift by a factor s 1 of a few hundreds should be enough to correct the empirical MBH MFs for obscuration and inactivity. In addition, there is the uncertainty in the MBH mass estimates themselves coming from the poorly determined virial factor, f σ or f FWHM . For instance, the MBH masses derived by Willott et al. (2010b) and Vestergaard et al. (2008) at z = 6.1 and 3.3, respectively, from emission-line re-verberation assume f FWHM = 1.17 (Collin et al. 2006), while, according to Yong et al. (2016), this factor could be up to a factor 5 smaller. Consequently, the empirical MBH MFs derived by those authors might also need to be shifted horizontally towards large masses by a factor s 2 in the range 1 ≥ s 2 ≥ 0.2. The value of s 2 is unknown, but it should be independent of z. On the other hand, the shifts s 1 /s 2 these effects would cause to the MBH mass density estimates derived by Kelly et al. (2010) and Willott et al. (2010b) should make them match the predicted MBHH curve consistent with the observed MBHto-spheroid mass ratio, equal to ∼ 400 and ∼ 250 at z = 3.3 and z = 6.1, respectively (see Fig. 10). 27 Taking into account the different degrees of completion in the extrapolated MFs used in the calculation of those MBH mass density estimates, the vertical shifts to be applied to the corresponding MFs at z = 3.3 and z = 6.1 should be 1.58 and 1.70 times larger, respectively. Therefore, with the only freedom of the z-independent value of s 2 (with 1 ≥ s 2 ≥ 0.2), it should be possible to shift the empirical MFs at the two redshifts so as to match those predicted at the two redshifts. As shown in Figure 11, this is indeed what we find for s 2 = 0.33 (or f FWHM = 0.39). Notice that, as a consequence of the small correction in f FWHM , the upper mass limit of MBHs at z = 3.3 becomes ∼ 10 10 M ⊙ , which is more reasonable than the original value of ∼ 3 × 10 10 M ⊙ if we compare with such an upper limit at z = 0.

Reliability of the Model
It might be argued that the model could be wrong and the solutions found for some appropriate values of the free parameters still recover the observed global and differential cosmic properties at z > 2. But such good results have been achieved with less than one degree of freedom per curve, which is unbelievable for any model that would not include the right physics. To see this we must simply think, for instance, in the fitting by means of a well-suited function with no physics of one single monotonous curve such as the observed SFH at z > 2 that requires four parameters (e.g. Robertson et al. 2015). Furthermore, not only does our model with only 9 degrees of freedom provide very good fits to the dozen observed curves at z > 2 (plus the galaxy stellar and MBH MFs), but it also automatically recovers all their extreme and saddle points at z ∼ 2 with no need to change the value of any parameter.
As mentioned, disks become increasingly abundant as z decreases, so the cold gas mass density increases (Fig. 7,  panel (a)). But cooling becomes increasingly inefficient, and by z ∼ 7 the rate of gas storage in disks begins to decline and the CGH becomes increasingly shallower. This effect is particularly marked in model D, where the CGH reaches a plateau between z ∼ 5 to z ∼ 3. The observed CGH is also flat at those redshifts, but the plateau reaches z ∼ 2 (see Fig. 7, panel (b)). There are several possible causes for that small discrepancy, such as the fact that the empirical cold gas mass density estimates do not include molecular gas or that observations may be biased towards large (massive) DLAs which have the lowest gas content. But the most likely reason is that satellite interactions, neglected in the present work, begin to have a significant impact in the amount of gas present in disks. In any event, the predicted stagnation of cold gas around z ∼ 4 corresponds to a saddle point, not to a maximum, in agreement with observation (e.g. Péroux et al. 2003;Prochaska, Herbert-Fort, & Wolfe 2005). At lower redshifts the temperature of the singly and doubly ionized IGM begins to diminish (see Fig 7, panel (l)), causing a rather constant slope in the HGMH, and the continuation of cooling.
Of course, the SFR density of ordinary (Pop II ) stars increases in parallel to the cold gas mass density (see Fig. 9), so, when gas replenishment begins to decline by z ∼ 7, the SFR density also slows down. Until that redshift, spheroids were much more abundant than disks, which never collected large amounts of gas. Nonetheless, disks were small, so their dynamical time was short, and the SFR density in disks was notable (see Fig. 12). The small amount of gas in disks also prevents the SFR density in mergers from being too large. It is just comparable to that in unstable disks at high-z, and reaches a maximum at z ∼ 6. This does not mean that mergers are less frequent from that moment; they just become increasingly dry. As the cooling rate declines, the SFR in individual unstable disks also does, but the corresponding global SFR density keeps on increasing because of the increasing abundance of spheroids contributing to it. The same happens with disks: even though the size of disks increases and the amount of cold gas they harbor decreases with decreasing z so that the SFR in individual disks decreases, the global SFR density in disks keeps on increasing because of the increasing abundance of disks. As a result, the ratio between the SFR densities in disks and spheroids is kept roughly constant. However, the declining amount of cold gas in disks ultimately causes the SFR in both unstable and stable disks to reach a maximum at z ∼ 2, and then to rapidly diminish towards z = 0 (see Fig. 12). Therefore, the origin of this maximum in the SFH is the same as for the plateau in the CGH: the lack of gas replenishment. However, in the SFH case, such a stagnation leads to a real maximum, not to a saddle point like in the CGH case. The drop in the SFR in disks and spheroids at z ∼ 2 is so marked that, even though at lower redshifts the cold gas mass density begins to increase again, the total SFR density will keep on decreasing until z = 0.
The stagnation of cold gas mass at z 5 has a similar effect on the SFH in models S and D. But the rest of the curve is slightly different in the two models. In model D, the SFH grows more steeply at high-z despite the smaller value of α G (see Fig. 7, panel (j)) because of the more rapid growth of ionized regions where normal galaxies lie. But, near z = 10 the formation of Pop III stars declines and ionized regions stop growing, causing the SFH to become shallower. This yields a concave curvature in the SFH of model D that perfectly fits the data. Instead, the SFH in model S has a more constant slope. This difference in the SFH of both models combined with the large abundance of spheroids with a very efficient SN heating also translates into a larger amount of hot gas in model D between z ∼ 12 to z ∼ 5 (see Fig. 7, panel (x)). Unfortunately, there are no observational data on this quantity yet, so the behavior of the HGMH in the two models cannot be checked.
The behavior of the IEH from galaxies and AGN follows from the evolution of those two kinds of sources described above. At very high-z, Pop III stars have the dominant contribution. But this contribution diminishes as neutral pristine regions progressively disappear. In model D it completely vanishes at z ∼ 10 when Pop III stars stop forming, whereas in model S it continues until z ∼ 6. On the other hand, the maximum contribution from AGN is one order of magnitude smaller at z ∼ 2.
Consequently, the H I-ionizing emissivity at z 6 is always dominated by normal galaxies, and its evolution closely follows the SFH associated with ordinary stars. In particular, the IEH also reaches a maximum at z ∼ 2 (see Fig. 13). What is less obvious is the maximum in both models S and D of the H I-and He II -ionizing (as well as the X-ray) emissivities from AGN at that redshift (Fig. 13): the MBH mass density is ever increasing at z < 2 (see Fig. 10) so we could naively expect their emissivities to also do. What causes these maxima is the fact that the most massive MBHs stop being fueled by that redshift (their galaxies undergo dry mergers). Consequently, the bright AGN stop radiating as quasars, and remain in the more quiescent radio phase. This has a major impact in the global AGN emissivities because bright AGN have the largest bolometric luminosity (e.g. Hatziminaoglou et al. 2003). However, massive MBHs keep on growing through dry mergers and less massive ones through wet mergers of low-mass galaxies that still harbor gas.

SUMMARY AND DISCUSSION
Observational data on the high-z universe are nowadays numerous and accurate enough to severely constrain the interconnected evolution of galaxies and IGM. In fact, they would completely determine it provided the IGM metallicity and temperature were known at z 6. These pieces of information are needed, indeed, to determine the threshold metallicity Z c for atomic cooling and the Pop III star IMF that play a crucial role in the reionization process. Assuming that IMF a power-law, the conditions that the metallicity of ionized regions must be high enough for triggering sustained galaxy formation and Pop III stars massive enough for seeding MBHs, greatly constrain those two properties.
The only acceptable solutions we find are in two narrow sets, one for small M lo III values (moderately top-heavy Pop III star IMFs) leading to single reionization, and another one for large M lo III values (top-heavier Pop III star IMFs) leading to double reionization. The latter solutions have a first ionization phase driven by Population III stars completed at z ∼ 10, and after a short recombination period a second ionization phase driven by normal galaxies completed at z ∼ 6. In the former solutions with single reionization, both kinds of luminous sources contribute in parallel to the ionization of the IGM, completed at z ∼ 6.
The best solution with double reionization gives excellent fits to all the observables. Instead, the best solution with single reionization predicts a too monotonous SFH, and a slightly too steep CGH at z ∼ 4. However, the CMB Thomson optical depth found in double reionization (τ = 0.102 +0.001 −0.001 ), though consistent with the observational estimates from the WMAP 9-year and Planck three-year data, is 3σ larger than the most recent estimate (τ = 0.058 ± 0.014) from the Planck data taking into account the large-scale polarization anisotropies. Instead, the CMB optical depth found in single reionization (τ = 0.072 +0.004 −0.005 ) is consistent with that estimate. There is thus some tension between the new estimate of τ and the observed properties of the Universe at high-z. However, the empirical values of the CMB optical depth derived under the instantaneous reionization approxima-tion must be taken with caution.
The fact that with only 9 degrees of freedom it is possible to fit a dozen of independent cosmic histories with non-trivial features (extreme and saddle points) confers strong reliability to our self-consistent model of the EoR.
As a byproduct, we have clarified the origin of some interesting features shown by the data on the high-z Universe: 1) the very small sizes of spheroids at z ∼ 2 − 3; 2) the sudden drop of the SFH at z < 2; 3) the similar behavior of the UV (and X-ray) emissivity from AGN at the same redshift; 4) the plateau in the CGH between at 1 < z < 5; and 5) the rough constancy of the MBH-tospheroid mass ratio. Our results also suggest that globular clusters are but the relics of the lowest mass galaxies with primordial origin or, more probably, formed during recombination in model D, which soon become frozen as satellite galaxies.
Finally, we have provided several predictions that should allow one to confirm the validity of the present model, and tell between single and double reionization. The most interesting one is that, in the case of double reionization, LAEs should be detectable near z = 10. This could explain the recent detection (to be confirmed) of a galaxy with Lyα emission at z = 8.68 (Zitrin et al. 2015).