© ESO 2018

1 Introduction

Around redshift z ≃ 1090, during the recombination era, protons paired with free electrons to form neutral atoms: the ionisation level of the intergalactic medium (IGM) fell to 0.0001% and remained at this level for several billions of years (Peebles 1968; Zel’dovich et al. 1969; Seager et al. 2000). Nevertheless, observations of the Gunn–Peterson effect (Gunn & Peterson 1965) in quasar spectra inform us that at z ~ 6, 99.96 ± 0.03% of the IGM hydrogen atoms are ionised (Fan et al. 2006). What happened in the meantime, during the Epoch of reionisation (EoR), is an essential source of information about the evolution of the Universe, the formation of large cosmic structures and the properties of early galaxies, to cite only a few. Thanks to improved observations of the cosmic microwave background (CMB), luminosity functions of galaxies, damping wings of quasars and Ly-α emissions (e.g. Schenker et al. 2013; Schroeder et al. 2013; Madau & Dickinson 2014; Planck Collaboration Int. XLVII 2016), more and high-quality data are available. Now, the generally accepted scenario is that, the first star-forming galaxies reionised neutral regions around them between z ≃ 12 and z ≃ 6 and then the ionised regions progressively overlapped (e.g. Aghanim et al. 1996; Becker et al. 2015) so that IGM neutral hydrogen fraction rapidly decreased until quasars took over to reionise helium from z ≃ 3 to 4 (Mesinger 2016).

Yet, some doubts remain about the sources of reionisation: some support the hypothesis that quasars could have led the process(Madau & Haardt 2015; Khaire et al. 2016; Grazian et al. 2018) but star-forming galaxies are often preferred. For instance, Robertson et al. (2015) argue that they were sufficient to maintain the IGM ionised at z ~ 7. The most recent value of the integrated Thomson optical depth, deduced from observations of the CMB, equals τ_Planck = 0.058 ± 0.012 and is obtained considering an instantaneous reionisation of duration δz = 0.5 (Planck Collaboration Int. XLVII 2016). It is much lower than previous observations by the Wilkinson Microwave Anisotropy Probe (WMAP) τ_WMAP = 0.088 ± 0.014 for z_reio = 10.5 ± 1.1 (Hinshaw et al. 2013). This decrease, according to Robertson et al. (2015), reduces the need for a significant contributionof high-redshift galaxies and allows them to extrapolate galaxies luminosity functions for 10 ≲ z ≤ 30.

Like Robertson et al. (2015), a number of recently published papers assume redshift-independent values of the escape fraction of ionising photons f_esc and of the clumping factor C_H II (Bouwens et al. 2015a; Ishigaki et al. 2015; Greig & Mesinger 2017), which is a questionable hypothesis. The escape fraction depends on numerous astrophysical parameters and, for this reason, it is often a generalised, global and redshift-independent value that is used, for an order of magnitude of 0.1. Some simulations give expressions of f_esc as a function of redshift (Haardt & Madau 2012; Kuhlen & Faucher-Giguère 2012) or of various parameters such as halomass or star formation rate (Wise et al. 2014; Paardekooper et al. 2015), but these models are rarely combined with observational constraints, aiming to deduce a certain history of reionisation. The situation is similar for the clumping factor: its evolution with redshift can be considered in simulations through various models (e.g. Mellema et al. 2006; Pawlik et al. 2009; Sobacchi & Mesinger 2014), but these are rarely compared with observations. We must, however, refer to Price et al. (2016) who constrain parametrised models of the escape fraction f_esc (z) with Thomson optical depth and low multipole E-mode polarisation measurements from Planck Collaboration Int. XLVII (2016), SDSS BAO data and galaxy observationsfor 3 ≲ z ≲ 10.

We first describe in Sect. 2 the observables of the reionisation process we will use throughout the analysis: the cosmic star formation rate density, the ionised fraction of the IGM and the Thomson optical depth, for which observational data is available – described in Sect. 3; as well as the two key-parameters of this study, the escape fraction of ionising photons and the clumping factor of IGM ionised hydrogen. Then we look in Sect. 4for the redshift-evolution we will further consider for the star formation rate (SFR) density, extrapolating luminosity functions at z ≳ 10. Doing this, we study the impact of our observational constraints on ρ_SFR. Investigations are then made on the escape fraction value and on how observations can constrain it: we try several parametrisations out – a redshift-independent one, where f_esc is free to vary in [0.1, 0.4], and a power-law function of z. We proceed the same for C_H II, but this time considering several possible parametrisations of its evolution with redshift, mainly from Iliev et al. (2007) and Pawlik et al. (2009). We conclude with a discussion of our results in Sect. 5, including a test of different values for the magnitude limit, and a summary in Sect. 6.

Throughout this paper, all cosmological calculations assume the flatness of the Universe and use the Planck cosmological parameters (Planck Collaboration XIII 2016): h = 0.6774, Ω _m = 0.309, Ω _b h² = 0.02230, and Y _p = 0.2453. Unless otherwise stated, all distances are comoving.

2 Observables of reionisation

2.1 Drawing the history of reionisation

Clues about the reionisation process can be derived from various observables. Under the assumption that star-forminggalaxies provided the majority of the photons which ionised the IGM, the star formation rate density, ρ_SFR, can logically give precious information about the EoR. Values of SFR density with redshift are deduced from luminosity functions (LF) of star-forming galaxies. LF can be observed down to a certain magnitude, but needs to be extrapolated to consider the contributionof unobserved fainter galaxies. Equation (1) shows how the comoving ionisation rate ṅ _ion is computed from the LF. $$\begin{array}{ll}{\dot{n}}_{\text{ion}}\hfill & ={\displaystyle {\int }_{M_{\text{lim}}}^{\infty }\varphi }(M_{\text{UV}})f_{\text{esc}}(M_{\text{UV}}){\xi }_{\text{ion}}(M_{\text{UV}})\text{d}{M}_{\text{UV}}\hfill \\ \hfill & \simeq \langle f_{\text{esc}}{\xi }_{\text{ion}}\rangle {\displaystyle {\int }_{M_{\text{lim}}}^{\infty }\varphi }(M_{\text{UV}})\text{d}{M}_{\text{UV}}\hfill \\ \hfill & \simeq f_{\text{esc}}{\xi }_{\text{ion}}{\rho }_{\text{SFR}}.\hfill \end{array}$$(1)

The final expression directly relates ρ_SFR to the cosmic reionisation rate ṅ_ion, in units of photons per unit time per unit volume, and is the version we will use in our models. We see that the choice of M_lim is fundamental as it directly impacts the value of ρ_SFR. Bouwens et al. (2015a) state that faint galaxies must contribute to the total UV radiation from galaxies but, assuming they do not form efficiently for lower luminosities (see Rees & Ostriker 1977; Mac Low & Ferrara 1999; Dijkstra et al. 2004), Robertson et al. (2015) choose to use M_lim = −13 rather than M_lim = −17, a choice we will discuss in this paper.

Two important parameters are used in Eq. (1): f_esc and ξ_ion. They describe the fact that only a limited amount of the photons produced by star-forming galaxies eventually end up ionising the IGM: first, they need to have sufficient energy – above the Ly-α limit, and second, they must escape their host galaxy and reach the IGM. The first condition is conveyed by ξ_ion, the quantity of Lyman continuum photons produced per second and per unit SFR for a typical stellar population. According to Robertson et al. (2015), we take ξ_ion = 10^53.14 Lyc photons s${}^{-1} M_{\odot }^{-1}$ yr. The second condition is conveyed by f_esc, the fraction of ionising radiation coming from stellar populations which is not absorbed by dust and neutral hydrogen within the host galaxy and so does contribute to the process. We note that in Eq. (1) we chose to consider values of f_esc and ξ_ion averaged over magnitude, i.e. the effective values.

Aiming to reproduce observations on the star formation history from z ~ 30 to z ~ 1, we choose the four-parameter model suggested by Robertson et al. (2015), updated from Madau & Dickinson (2014, Sect. 5, Eq. (15)) and described in Eq. (2) below. According to data, ρ_SFR (z) follows a first rising phase, over 3 ≲ z ≲ 15, which is expressed in our parametrisation by an evolution ${\rho }_{\text{SFR}}(z)\propto {(1+z)}^{b-d}$, up to a peaking point around z ~ 1.8, that is, when the Universe was around 3.6 Gyr old. It then declines as ${\rho }_{\text{SFR}}(z)\propto {(1+z)}^b$ until z =0. To stay consistent with observations, we set b > 0 and b− d < 0. $${\rho }_{\text{SFR}}(z)=a\frac{{(1+z)}^b}{1+{\left(\frac{1+z}{c}\right)}^d}.$$(2)

In orderto put our results in perspective, we consider different values of the magnitude limit for our study and therefore use another parametrisation of the star formation history, suggested by Ishigaki et al. (2015) and designed to reproduce the rapid decrease of ρ_UV(z) from z ~ 8 towards higher redshifts and but not the bump on luminosity density observed around z ~ 2 $${\rho }_{\text{UV}}(z)=\frac{2{\rho }_{\text{UV}}(z=8)}{{10}^{a(z-8)}+{10}^{b(z-8)}}.$$(3)

Here, ρ_UV(z = 8) is a normalisation factor, and a and b characterise the slope of ρ_UV(z). This model is more adapted to the study of reionisation in itself, as the process is known to end before z = 4 and so before the star formation bump. However we cannot limit our analysis to this late-redshift model since the former carries more information about the star formation history and is therefore more interesting when considering a large amount of free parameters. We note that for M_lim = −10 and M_lim = −17, we use ξ_ion = 10^25.2 erg⁻¹ Hz, following Ishigaki et al. (2015).

Other observations can lead to estimations of the fraction of ionised IGM Q_H II, also called filling factor, which relates to the SFR density via Eq. (4). In this equation, the time-related evolution of Q_H II depends on two contributions: an ionisation source term, proportional to ṅ _ion, and a sink term due to the competition of recombination. t_rec is the IGM recombination time defined in Eq. (5) and ⟨n_H ⟩ is the mean hydrogen number density, defined by $\langle n_{\text{H}}\rangle =\frac{X_{\text{p}}{\Omega }_{\text{b}}{\rho }_{\text{c}}}{m_{\text{H}}}$, with ρ_c the critical density of the Universe. $${\dot{Q}}_{H_{II}}=\frac{{\dot{n}}_{\text{ion}}}{\langle n_{\text{H}}\rangle }-\frac{Q_{H_{II}}}{t_{\text{rec}}},$$(4) $$\frac{1}{t_{\text{rec}}}=C_{H_{II}}{\alpha }_{\text{B}}(T)\left(1+\frac{Y_{\text{p}}}{4X_{\text{p}}}\right)\langle n_{\text{H}}\rangle {(1+z)}^3.$$(5)

In Eq. (5), X_p and Y _p are the primordial mass fraction of hydrogen and helium respectively. α_B(T) is the case B recombination coefficient at a fiducial IGM temperature of T = 20 000 K, often considered as the mean temperature around a newly ionised atom. This value is consistent with measurements at z ~ 2−4 (Lidz et al. 2010) but has been estimated to T ≲ 10⁴ K at z ~ 5−6 (Becker et al. 2011; Bolton et al. 2012). It fluctuates by a factor of between one and two, depending on the spectrum of the sources and on the time passed since reionisation (Hui & Haiman 2003). Yet, α_B is expressed as ${\alpha }_{\text{B}}(T)\approx 2.6\times {10}^{-13}{T}_4^{-0.76}{\text{cm}}^3{\text{s}}^{-1}$ with T₄ = T∕10⁴ K (Osterbrock 1989), in other words, it is a weak function of T so that its variations do not affect our results significantly. We note that, rather than case A, we considered case B recombinations in order to exclude recombinations to the ground state and because we consider that ionisations and recombinations are distributed uniformly throughout the IGM, so that each regenerated photon soon encounters another atom to ionise (Loeb & Furlanetto 2013, Sect. 9.2.1). The clumping factor C_H II expresses how ionised hydrogen nuclei are distributed throughout the IGM. C_H II and t_rec are inversly proportional: the more the matter is aggregated in clumps, the easier for ionised atoms to recombine in these very same clumps. To compare with the evolution derived from Eq. (4), we considered two parametrisations of the time evolution of the filling factor Q_H II, that we willthen use to calculate the integrated Thomson optical depth from data. The first depicts the reionisation process as a step-like and instantaneous transition with a hyperbolic tangent shape (Eq. (6)). The second is a redshift-asymmetric parametrisation, described in Eq. (7), inspired by Douspis et al. (2015). It uses a power-law defined by two parameters i.e. the redshift at which reionisation ends z_end and the exponent α: $$Q_{H_{II}}(z)=\frac{f_{\text{e}}}{2}\left[1+\text{tanh}\left(\frac{y-y_{\text{re}}}{\delta y}\right)\right],$$(6) $$Q_{H_{II}}(z)=\{\begin{array}{ll}{f}_{\text{e}}\hfill & \text{for}z<z_{\text{end}},\hfill \\ f_{\text{e}}{\left(\frac{z_{\text{early}}-z}{z_{\text{early}}-z_{\text{end}}}\right)}^{\alpha }\hfill & \text{for}z>z_{\text{end}}.\hfill \end{array}$$(7)

where $y(z)={(1+z)}^{\frac{3}{2}}$, y_re = y (z = z_re) for z_re the redshift of instantaneous reionisation and $\delta y=\frac{3}{2}{(1+z)}^{\frac{1}{2}}\delta z$. z_early corresponds to the redshift around which the first emitting sources form, and at which Q_H II (z) is matched to the residual ionised fraction ($\overline{\text{x}}={10}^{-4}$). To be consistent with observations, which give Q_H II (z ≤ 6.1) ≃ 1 with very low uncertainty (McGreer et al. 2015), we choose z_end = 6.1. Furthermore, when comparing our findings with the Planck results we set z_re at equal to 8.5, z_early = 20, and also α =6.6 (Planck Collaboration Int. XLVII 2016).

Observations of CMB satellites allow us to estimate the Thomson optical depth τ, integrated over the electron column density to the last scattering surface. It expresses the fraction of photons scattered along the line of sight by free electrons and thus is a direct indicator of the global ionisation rate of the IGM. It is related to the two previously described observables Q_H II and ρ_SFR via Eq. (8), where c is the speed of light in vacuum, σ_T the Thomson scattering cross-section, H(z) the Hubble constant and f_e the number of free electrons per hydrogen nucleus. We have assumed that helium is doubly ionised at z ≤ 4 (Kuhlen & Faucher-Giguère 2012) and thus have f_e = 1 + ηY _p∕4X_p with η = 2 for z ≤ 4 and η = 1 for z > 4. $$\tau (z)=c\langle n_{\text{H}}\rangle {\sigma }_{\text{T}}{\displaystyle {\int }_0^z{f}_{\text{e}}}\frac{Q_{H_{II}}(z^{\prime })}{H(z^{\prime })}{(1+z^{\prime })}^2\text{d}{z}^{\prime }$$(8)

2.2 Configuring the key-parameters of reionisation

Among the various parameters cited in Sect. 2.1, two key-parameters of the reionisation history are still under a lot of investigations: the escape fraction and the clumping factor. As mentioned before, f_esc expresses the fraction of the ionising radiation produced by stellar populations which is not absorbed by dust and neutral hydrogen within its host galaxy, and thus contributes to the ionisation of the IGM. In our approach, it is an effective value, averaged over stochasticity, halo mass dependencies in the source populations and, most importantly, over all sources considered in the Universe. This averaged value is hard to compare with observations of lone galaxies or haloes, which usually give much lower values. For instance, Steidel et al. (2001) and Iwata et al. (2009) estimate the escape fraction of some z ~ 3 galaxies to be ≳ 1%. On the contrary, overall values of f_esc can be derived from simulations but are still highly uncertain. According to Finkelstein et al. (2015) and to agree with Ly-α forests measurements (Bolton & Haehnelt 2007), it should not be higher than 0.13; Fernandez et al. (2013) use a value of 0.1 from a simulation; Robertson et al. (2015) deduce from their analysis that, in order to have star-forming galaxies driving the reionisation process at high redshift, f_esc must equal at least 0.2; Inoue et al. (2006) find that, if recent values of the escape fraction can be as low as f_esc= 0.01 at z~ 1, f_esc increases quickly with redshift to reach 10% at z≳ 4. Finally, Dunlop et al. (2013) assure that, considering the spectral energy distributions observed from high-redshift galaxies, it should be ≈ 0.1−0.2. Yoshiura et al. (2017) summarise results on f_esc by saying that if it is generally acknowledged that, among all dependencies, the escape fraction decreases with the mass of the galaxy, there is a variance within one or two orders of magnitude among simulations results. For instance, a simulation from Yajima et al. (2014), on which assumptions of Robertson et al. (2015) are based, shows that, amidst all types of photons produced in star-forming galaxies (Ly-α, UV-continuum and ionising photons), the escape fraction of ionising photons is the only one which seems to depend neither on the redshift nor on the galaxy properties: it keeps a constant value of 0.2 with time, that we use for our first analysis.

However, photons from different ranges of energy are subject to different physical phenomena and thus escape more or less easily from their host galaxy. For instance, dust extinguishes ionising, Ly-α and UV continuum photons similarly, but only ionising photons are also absorbed by neutral hydrogen clumps. Thus, at high redshifts, when there is little dust around the galaxy, photons of all energy ranges escape as easily; on the contrary, at low redshift, ionising photons experience more difficulties to escape than others (Yajima et al. 2014). We can then infer an increase of f_esc with redshift that we parametrise in Eq. (9), defined for z ≥ 4 and inspired by Kuhlen & Faucher-Giguère (2012); Chisholm et al. (2018). This evolution corresponds to either an evolution of the SFR of galaxies themselves and its associated feedback, or by a redshift evolution in the make up of the galaxy population. Here, owing to the UV spectral slope constraints, we set a maximum of 1 for the f_esc value, corresponding to a situation where all ionising photons escape. $$f_{\text{esc}}(z)=\alpha {\left(\frac{1+z}{5}\right)}^{\beta }.$$(9)

In this parametrisation, also close to the one used in Price et al. (2016), α is the value of f_esc at z = 4 and αβ∕5 of its derivative at z = 4, redshift at which we expect the hydrogen ionising background to be dominated by star-forming galaxies (Kuhlen & Faucher-Giguère 2012). We take β positive in order to have an increasing escape with redshift, as anticipated earlier.

The second key-parameter of the reionisation process which we are going to investigate is the clumping factor of ionised hydrogen in the IGM C_H II, used in Eq. (5). It expresses how ionised hydrogen nuclei are gathered in heaps throughout the IGM. This parameter is essential because it is the growth of these clumps that allows the reionisation front to progress in the IGM and because competing recombinations will predominantly take place there. A precise estimate of C_H II can be difficult to obtain. Simulations do indeed have several obstacles to overcome: getting a sufficient precision for the gas distribution, a correct topology of ionised and neutral matter, and an accurate model of the evolution of gas clumps themselves during the reionisation process. Besides, C_H II is often firstdefined on a single ionisation bubble and then summed on all bubbles to get the global volume-averaged value used here: the simulation must consider an extremely wide range of scales (Loeb & Furlanetto 2013, Sect. 9.2).

Most recent studies use values ranging from one to six at the redshifts of interest, i.e. for 6 ≲ z ≲ 30 (Sokasian et al. 2003; Iliev et al. 2006; Raičević & Theuns 2011; Shull et al. 2012; Robertson et al. 2015; Finkelstein et al. 2015; Bouwens et al. 2015a). Other studies predict a redshift-dependent evolution (Iliev et al. 2007; Pawlik et al. 2009; Haardt & Madau 2012; Finlator et al. 2012; Sobacchi & Mesinger 2014), justified by the fact that during the late stages of EoR, ionisation fronts penetrate into increasingly overdense regions of the IGM, which have higher recombination rates and so drive a rapid increase of C_H II (Furlanetto & Oh 2005; Sobacchi & Mesinger 2014). In our study, besides constant values of C_H II, we consider two parametrisations on the redshift range 3 ≤ z ≤ 30¹: $$\begin{array}{lll}\hfill & \hfill & C_{H_{II}}(z)=\alpha +a{\left(\frac{z}{8}\right)}^b,\hfill \\ \hfill & \hfill & C_{H_{II}}(z)=a{\text{e}}^{b(z-8)+c{(z-8)}^2}.\hfill \end{array}$$

The first expression comes from Haardt & Madau (2012). We update it in order to have a = C_H II (z = 8) − α because Q_H II is close to 0.5 at z = 8. The second one comes from Mellema et al. (2006) and Iliev et al. (2007) and shows a different behaviour: it is convex and has a minimum at z_min = −b∕2c. As explained earlier, it is generally admitted that the clumping factor only decreases with z, and therefore we set z_min ≳ 30 so that C_H II does not reach its minimum on our analysis range. For the same reason, a and b from Eq. (10) have to be of opposite signs and more precisely we take a > 0 and b < 0 in order to have $C_{H_{II}}(z)\underset{z\to 0}{\to }+\infty$.

The formal definition of the clumping factor is (Bouwens et al. 2015a; Robertson et al. 2015): $C_{H_{II}}=\langle n_{H_{II}}^2\rangle /{\langle n_{H_{II}}\rangle }^2=1+{\delta }_{H_{II}}$, if we define the overdensity of ionised hydrogen as ${\delta }_{H_{II}}=\left(n_{H_{II}}-\langle n_{H_{II}}\rangle \right)/\langle n_{H_{II}}\rangle$. Long beforethe EoR, most of the hydrogen was neutral so that fluctuations in the ionised hydrogen overdensity were very weak. In this perspective, we consider in our models that δ_H II(z →∞) = 0 and so take C_H II(z = 100) = 1.

3 Data

The SFR density can be estimated via the observed infrared and rest-frame UV LFs. We use the luminosity densities and SFR densities compiled by Robertson et al. (2015), computed from Madau & Dickinson (2014), Schenker et al. (2013), McLure et al. (2013), Oesch et al. (2015), and Bouwens et al. (2015a). Robertson et al. (2015) also use HST Frontier Fields LF constraints at z ~ 7 by Atek et al. (2015) and at z ~ 9 by McLeod et al. (2015). Estimates of Madau & Dickinson (2014) derived from Bouwens et al. (2012) are updated with newer measurements by Bouwens et al. (2015a). For the calculation of ρ_SFR, as a start, luminosity functions of star-forming galaxies are extended to UV absolute magnitudes of M_lim = −13. Then we compared this with results for minimal and maximal magnitude limits M_lim = −17 and M_lim = −10. We note that if Robertson et al. (2015) express ρ_SFR in M_⊙ yr⁻¹ Mpc⁻³, Ishigaki et al. (2015) use UV luminosity units, i.e. ergs s⁻¹ Hz⁻¹ Mpc⁻³. In order to compare results, we used the conversion factor used in Bouwens et al. (2015a) and first derived by Madau et al. (1998): $$L_{\text{UV}}=\frac{\text{SFR}}{M_{\odot }{\text{yr}}^{-1}}\times 8.0\times {10}^{27}\text{ergs}{\text{s}}^{-1}{\text{Hz}}^{-1}.$$

UV luminosity densities used in this work are the ones detailed in Ishigaki et al. (2015), namely they come from Schenker et al. (2013); McLure et al. (2013); Bouwens et al. (2007, 2014, 2015b); Oesch et al. (2015).

Observations related to the ionised fraction of the IGM Q_H II used as constraints to our fits include the Gunn-Peterson optical depths and the dark-gap statistics measured in z ~ 6 quasars (McGreer et al. 2015), damping wings measured in z ~ 6−6.5 quasars (Schroeder et al. 2013) and the prevalence of Ly-α emission in z ~ 7−8 galaxies (Schenker et al. 2013; Tilvi et al. 2014; Faisst et al. 2014). We note that in the figures, further data points, not used as constraints in the fit, are displayed to use as comparison. These include observations of Lyman-α emitters (Konno et al. 2018; Ouchi et al. 2010; Ota et al. 2008; Caruana et al. 2014), of near-zone quasars (Mortlock et al. 2011; Bolton et al. 2011), and of a gamma-ray burst (Chornock et al. 2014).

Last, we consider estimations of the Thomson optical depth derived from Planck Satellite observations: τ_Planck = 0.058 ± 0.012 for a redshift of instantaneous reionisation z_reio = 8.8 ± 0.9 (Planck Collaboration Int. XLVII 2016). We compare it to the asymptotic value τ obtained from our model calculations at high redshift.

4 Results

4.1 Cosmic star formation history

Since we are interested in the reionisation history both up to and beyond the limit of the current observational data, we adopt the four-parameter model from Eq. (2) into a Monte Carlo Markov chain (MCMC) approach. We perform a maximum likelihood (ML) determination of the parameter values assuming Gaussian errors on a redshift range of [0, 30], extrapolating current observations on star formation history from z = 10.4 to z = 30. We fit to the star formation data described in Sect. 3 and then compute the range of credible reionisation histories for every value of the ρ_SFR model parameters by solving the differential equation of Eq. (4). Filling factor data is used as an additional observational prior for the fit. Finally, we evaluated the Thomson optical depth as a function of z via Eq. (8) and compare its “asymptotic” value, at z = 30, to τ_Planck = 0.058 ± 0.012 (Planck Collaboration Int. XLVII 2016) as a last constraint on the fit. Because we want to know what observable constrains reionisation history the most, all constraints are not always used: the run ALL uses all three sets of data as constraints; NOQ skips Q_H II data; NORHO skips star formation data, and ORHO uses only star formation history in the fit.

In this first step, we adopt the fiducial, constant with redshift values f_esc= 0.2, log₁₀ ξ_ion = 53.14 [Lyc photons ${\text{s}}^{-1}{M}_{\odot }^{-1}$ yr] and C_H II = 3 (e.g. Pawlik et al. 2009; Shull et al. 2012; Robertson et al. 2013, 2015). Results are summarised in Fig. 1 and in Table 1. Figure 2a shows resulting star formation history and Fig. 2b resulting reionisation history. We find that star formation history constrains reionisation the most: both figures show that ALL and ORHO runs give similar evolutions and close ML values for a, b, c, and d (see Table 4). We note that our constraints with ORHO and ALL are dominated by the ρ_SFR data points at a redshift of approximately five and the fixed functional form assumed for ρ_SFR(z); they are fully consistent with Robertson et al. (2015). On the contrary, for NORHO, the shape of ρ_SFR(z) is changed and reionisation begins muchlater, around z ~ 12 rather than z ~ 15 for other runs. NORHO results must be handled carefully as its parameters probability density functions (PDFs) are extremely spread-out; the NORHO line drawn on figures corresponds to the median values of parameters. All we can conclude is that, when star formation history constraints are skipped, there is a much wider range of possible scenarios.

Interestingly, Fig. 2b shows that for each run considering star formation history constraints, the process begins as early as z = 15. This is hardly compatible with WMAP results which stated that, if we consider reionisation as instantaneous, it should occur at z_reio ≃ 10.5 ± 1.1 (Hinshaw et al. 2013) and so cannot begin before z = 12. Observations also have an influence on the Thomson optical depth values, as NORHO gives a slightly lower value of τ (0.053 ± 0.003 compared to 0.061 ± 0.001 for ALL). Yet, all results remain in the 1-σ confidence interval of τ_Planck.

In the rest of the study we used the ALL run as our definitive parametrisation for ρ_SFR evolution with redshift: definitive parameters for Eq. (2) are (a = 0.146, b = 3.17, c = 2.65, d = 5.64) from Table 1. ML parameters for other runs can be found in Table 4.

Fig. 1 Results of the MCMC analysis for the ALL case. The contours show the 1, 2, 3σ confidence levels for a, b, c, d and the derived parameter τ.

Fig. 2 Panel a: star formation rate density ρSFR with redshift. Data points are determined from infrared (plotted in red) or ultraviolet (in grey) luminosity densities (Sect. 3). Maximum likelihood parametrisations (continuous lines) are shown for various set of constraints: blue when all constraints are used; coral when only data on star formation are used; green when τ and reionisation history data are used. The 68% confidence interval on ρSFR (light blue region) is drawn for the blue model. We note that the interval, corresponding to statistical uncertainties, is very narrow. These inferences are compared with a model forced to reproduce results from Robertson et al. (2015), cited as R15 in the legend, drawn as the purple dotted line. The horizontal dashed-dotted line corresponds to the upper limit on a hypothetical constant value of ρSFR for z > 10.4 (Sect. 5.2). Panel b: ionised fraction of the IGM QH II with redshiftfor same models as (a). Details on the origin of data points are given in the legend and Sect. 3. Inferences are also compared with the two evolutions used in Planck Collaboration Int. XLVII (2016, cited as P16) to model the reionisation process: a redshift-symmetric hyperbolic tangent as the brown dashed-dotted line and a redshift-asymmetric power-law in black.

Fig. 3 Ionised fraction of the IGM QH II with redshiftwhen fesc is introduced as a parameter. Details on the origin of data points are given in the legend. ML models (continuous lines) are shown for various set of constraints: blue when all constraints are used, coral when QH II constraints are skipped, for which the 68% and 95% confidence intervals are drawn in salmon. These inferences are compared with a model forced to reproduce results from Robertson et al. (2015, R15, purple dotted line) and with the two evolutions used in Planck Collaboration Int. XLVII (2016, P16): redshift-symmetric as the dashed-dotted brown line and redshift-asymmetric in black.

4.2 Escape fraction of ionising photons f_esc

In order to study the role of the escape fraction in this analysis we chose, as detailed in Sect. 2.2, to first consider it as a fifth parameter of the fit – on top of (a, b, c, and d) from Eq. (2), free to vary between 0 and 0.4. We name ALL the run which uses ρ_SFR, Q_H II and τ constraints, and NOQ the one that skips ionisation level constraints. f_esc is involved only in the ṅ _ion calculation of Eq. (1) and not in the one of ρ_SFR so that star formation history takes no part in the computation of f_esc. This explains why for all runs, results on the SFR density are close to the ones of Sect. 4.1 (see Tables 4 and 5 for details). For ALL, we get ML parameters (a = 0.0147, b = 3.14, c = 2.69, d = 5.74). Figure 3 shows that Q_H II constraints have a strong influence on f_esc: confidence intervals are much wider for NOQ than for ALL (see Table 5). Besides, the NOQ PDF of f_esc is almost flat: standard deviation is equal to 0.079, that is, around 30% of the mean value and two times more than for ALL. For now, we chose to use f_esc = 0.19 ± 0.04, in other words, the median value of the escape fraction for the ALL run, when a redshift-independent value is needed for f_esc. The full triangle plot for the ALL case is shown in Fig. A.1.

We now turn to the possibility of a redshift evolution in f_esc for z ∈ [4, 30]. We performan MCMC maximum likelihood sampling of the two-parameter parametrisation described inEq. (9). For the reasons explained above on the lack of relation between ρ_SFR and f_esc, we do not use star formation data as a constraint any more and assume that the time evolution of the SFR density follows Eq. (2) using parameters (a, b, c, d) resulting from Sect. 4.1. We used parameters corresponding to the set of constraints that is used on f_esc: if only τ priors are considered here, we use (a, b, c, and d) resulting from a NOQ run (see Table 4 for values).

We findthat priors on the IGM ionisation level have a much stronger influence on results than the Thomson optical depth. Indeed, Fig. 4 shows that ML evolutions using both Q_H II and τ constraints or only Q_H II are very similar: mean values for z ≥ 4 are similar by ~3% and in bothcases, the evolution with redshift is rather weak, as values range from 0.15 around z ~ 4 to 0.24 around z ~ 30. We note that if Mitra et al. (2015) draw a similar conclusion of an almost constant f_esc value with redshift from their modelling, they obtain lower values of the escape fraction with an average of about 10% in the redshift range six to nine. For NOQ, the optical depth remains surprisingly close to other models and to τ_Planck = 0.058 ± 0.012, around 0.061. The difference is apparent in the evolution of the ionised fraction, as reionisation begins and ends later, around z = 6 rather than z = 6.4 in this case; on the contrary, when Q_H II data is used, the history tends to be the same as in previous analysis. Our results when only τ_Planck constraints are considered are quite similar to those of Price et al. (2016) in which authors study the evolution of f_esc with redshift. They mainly use constraints from τ_Planck, concluding to a strong increase of f_esc from about 0.15 to about 0.55, depending on the observational constraints used. ML parameters for Eq. (9) when all constraints are considered are (α = 0.14 ± 0.02, β = 0 ± 0.3) and give a mean value for f_esc of about 0.20, which is extremely close to the 0.19 ± 0.04 found when considering the escape fraction constant with redshift (see Table 5 for details).

Fig. 4 Possible evolutions of fesc with redshift. ML models are shown for various set of constraints: blue when all constraints are used; coral when τ constraints are skipped; green when QH II constraints are skipped. Horizontal dotted lines represent the mean value of fesc over 4 ≤ z ≤ 30 for the model of the corresponding colour.

4.3 Clumping factor of ionised hydrogen in the IGM C_H II

Following the definition of Sect. 2.2, we now investigate the constraints on C_H II set by observations. As we did in Sect. 4.2 for f_esc, we added C_H II as a fifth parameter of the fit on ρ_SFR using Eq. (2), apart from (a, b, c, d). It is free to vary between zero and ten, the order of magnitude of fiducial values most commonly used in publications (e.g. Shull et al. 2012; Robertson et al. 2013, 2015). Here again, we call ALL the run using all constraints in the fit, and NOQ the one that skips Q_H II constraints.

After performing the MCMC ML sampling of the five parameters (see Table 4 for details), we get a quite spread PDF for C_H II with ALL: the standard deviation is equal to 1.85 for a medianvalue of 4.56. Even with such a wide range of possible values, the range of possible reionisation histories remains very narrowand the Thomson optical depth PDF is almost exactly the same as when we take C_H II = 3: τ_ALL = 0.0570 ± 0.0019 to be compared with ${\tau }_{C_{H_{II}}=3}=0.0612\pm 0.0013$ (see Table 2 and Fig. A.2). Besides, for NOQ, the range of possible reionisation histories is wider than for ALL: the width of the 95% confidence area is about 0.6 when ML reionisation model is halfway through (Q_H II = 0.5) for NOQ but 0.16 for ALL. We also note that for NOQ τ takes lower values (τ_NOQ = 0.0561 ± 0.0064) but remains, as others, in the 1-σ confidence interval of τ_Planck. This confirms that IGM ionisation level data are compatible with Planck observations and that the value of C_H II constrains only slightly the reionisation history.

We now successively test the two redshift-dependent models of the clumping factor given in Eqs. (10) and (11). C_H II is not involved in the calculation of ρ_SFR but only of the recombination time. Thus, as for f_esc, star formation history data have no influence over it: the ALL run is now constrained by Q_H II and τ_Planck only. We also note that, for low values of z (precisely for z ≤ 6.8), Q_H II becomes higher that 1 in our calculations, which is physically irrelevant so we ignore results in this range.

Once again, IGM reionisation level data constrain results more than τ_Planck. The redshift-evolution of C_H II and Q_H II for the two parametrisations presented in Sect. 2.2 and for ALL runs are shown in Fig. 5a and b. We see on the left panel that there are a lot of possible output evolutions for both models but this does not translate in significant variations of Q_H II(z) whose 68% confidence intervals are found to be very narrow. All scenarios remain quite close, with reionisation beginning around z = 16 and ended by z = 6. This means that, as in previous paragraph where C_H II was assumed constant with redshift, its exact value has no significant impact on the reionisation history. In fact, variations in C_H II have some impact on the computed Thomson optical depth: as seen in Table 2, higher values of C_H II allow for a lower value of τ – consistentwith Eqs. (5) and (8). All values remains in the 1-σ confidence interval of τ_Planck.

Finally, it seems that the fiducial constant value often used in papers, C_H II = 3, and which lies between the mean values of our models (~3 for HM12, 1.8 for I07, and 4.5 for Free), is a reasonable choice. More generally, and in accordance with Bouwens et al. (2015a), as long it remains in a range of [1.4, 8.6], which is the 95% confidence interval of C_H II from first paragraph (Free fit), results are consistent with the three sets of constraints available. This result corroborates the work of Price et al. (2016), who also note that their analysis is almost completely independent of the clumping factor over the prior range 1 < C_H II < 5.

Fig. 5 Panel a: possible evolutions of CH II with redshift. ML models are shown for the two models of Sect. 2.2: blue for the first, coral for the second. Dotted horizontal lines correspond to the mean value of CH II (z) for z > 6.8, where outputs of the model are used in calculations, for the model of the corresponding colour. The vertical line is located at z = 6.8. Lines of lighter colours represent various outputs of the sampling of the corresponding model. Panel b: redshift evolution of QH II for the same models of CH II(z). Inferences are compared to a result with CH II(z) = 3 in purple dashed line, and to the theoretical models of Planck Collaboration Int. XLVII (2016): a redshift-symmetric model in black and a redshift-asymmetric model in brown. References. HM12: Haardt & Madau (2012) or Eq. (10). I07: Iliev et al. (2007) or Eq. (11).

4.4 Varying both f_esc and C_H II

Now we have studied the impact of f_esc and C_H II separately, we set the evolution of ρ_SFR(z) according to Eq. (2), using parameters a, b, c, and d resulting from the analysis of Sect. 4.1. We performed an MCMC maximum likelihood sampling of the two parameters f_esc and C_H II, considered constant with redshift. The first is allowed to vary between 0.001 and 1, the other between one and seven. We show parameter distributions for f_esc and C_H II in Fig. 6. We constrain the fit with all three data sets.

If we consider the median value of each parameter distribution as its maximum likelihood value, we find f_esc= 0.193 ± 0.026 and C_H II = 4.43 ± 1.11. We see results are pretty similar to the previous analysis: if the escape fraction is well constrained, with a standard deviation of about 13%, the clumping factor can take a much wider range of values, between 2 and 6. We note that there seems to be a strong upper bound for the escape fraction around 0.26, which we can compare to the asymptotic value of f_esc when it is allowed to change with redshift (see Fig. 4). Because parameters take values close to previous results, the resulting ionisation histories are also close to the ones observed in Fig. 2b and are hence in good agreement with observations.

Finally, we considered the case when the four parameters describing the evolution of ρ_SFR(z) are set free in the same time as f_esc and C_H II, using all datasets. We assumed the same prior as Price et al. (2016) on C_H II considering values between one and five. The full triangle plot is shown in Fig. A.3 and best fit parameters are reported in Table 4. The values found are in agreement with previous runs, with an undetermined value of C_H II at the 2σ level. As in Price et al. (2016) the degeneracy between f_esc and C_H II and the current data do not allow to constrain strongly all free parameters. However the evolution of the filling factor (Fig. 7) and thus the derived value of τ remain quite well constrained (τ = 0.058 ± 0.002) and in agreement with Planck (τ_Planck = 0.058 ± 0.012).

Fig. 6 MCMC distribution for fesc and CHII when both are taken as fit parameters (other parameters fixed). The escape fraction is allowed to vary between 0.1% and 100%, the clumping factor between one and seven.

5 Discussion

5.1 Influence of the magnitude limit

In order to study the influence of the choice of magnitude limit on our results, we adopt the model of Eq. (3) into an MCMC approach similar to Sect. 4. We fit the model to our three data sets adapted to the corresponding magnitude limit as described in Sect. 3. M_lim = −17 and M_lim = −10 correspond to the analysis performed in Ishigaki et al. (2015), and M_lim = −13 corresponds to Robertson et al. (2015).

We compute the star formation and reionisation histories compatible with the three sets of observational data, for the maximum likelihood parameters (here, median values) of the parametrisation in Eq. (3) and for the three M_lim values. Results can be found in Fig. 8a and b where two cases have been considered: f_esc fixed, takento have the value used in corresponding references (left panels) and f_esc allowed to vary between 0 and 1 (right panels). In both cases, the effect of the two additional sets of data used as constraints here, Q_H II and τ_Planck, which were not used in Ishigaki et al. (2015), is to lower the quantity of ionising sources needed at high redshift to reach a fully ionised IGM by z ~ 6. We note, however, that here the values of some parameters were taken from Ishigaki et al. (2015) and hence quite different from the ones used in Sect. 4.1. For instance, Ishigaki et al. (2015) found C_H II values of 1.9 and 1.0 for respectively M_lim = −17 and M_lim = −10 whereas we used C_H II = 3 before and consequently in the analysis for M_lim = −13.

However, this comparison illustrates the systematic uncertainties on reionisation history due to the choice in the magnitude limit, but also in f_esc and C_H II values. We see these are much wider than the statistical uncertainties observed in Fig. 2b while still being reasonable. In particular, they mainly concern high redshifts. Indeed, we see in Fig. 8b that the 68% confidence interval on star formation histories widens with redshift. However, few observations are available on this redshift range so we may expect that once data on earlier times is available, we will be able to improve constraints on the magnitude limit. In this perspective we can mention the work of Mason et al. (2018), who derived a new constraint on reionisation history from simulations and models of the effects of IGM radiative transfer on Lyman-α emissions. They find an IGM ionised fraction at z ~ 7 of $x_{H_{II}}={0.41}_{-0.11}^{+0.15}$ in better agreement with our model for M_lim = −17 (see Fig. 8a).

From a different point of view, Price et al. (2016) consider a varying value of M_lim with redshift, and find that M_lim varies in order to match the value of τ_Planck and to balance the increasing value of f_esc with redshift allowed by their model. Here we find that, overall, the model combines star formation history and ionised fraction with difficulties when M_lim = −17. Indeed, Fig. 9 shows the probability distribution functions of the parameters log(ρ_SFR) and f_esc and the corresponding distribution of derived optical depths for the three choices of magnitude limit. We see that for M_lim = −17 the value of f_esc is not well constrained and tends to be high. For lower values of the escape fraction, the reionisation process needs to start way earlier than in most of our results in order to have enough radiation to fully ionise the IGM and to reach a sufficient value of τ. In fact, leaving the escape fraction as a free parameter balances the uncertainty in the choice of M_lim: Fig. 8a shows a narrower range of uncertainties when we do not fix f_esc, confirming the correlation mentioned in Price et al. (2016).

Fig. 7 Redshift evolution of QH II when all parameters (a, b, c, d, fesc, CH II) are free andall datasets used. Figure A.3 shows the corresponding constraints on assumed parameters.

5.2 Reionisation sources at z > 10

Some doubts remain about the sources of reionisation: if Robertson et al. (2015) found that star-forming galaxies are sufficient to lead the process and to maintain the IGM ionised at z ~ 7 – assuming C_H II = 3 and f_esc = 0.2, their analysis extrapolates luminosity functions between z ≃ 10 and z ≃ 30, overlooking the possibility that other sources may have taken part in the early stages of reionisation process. Besides, they argue that low values of the Thomson optical depth reduce the need for a significant contribution of high-redshift galaxies and Planck Collaboration Int. XLVII (2016) give much lower values than WMAP did (Hinshaw et al. 2013): τ_Planck = 0.058 ± 0.012 vs. τ_WMAP = 0.088 ± 0.014. Thus, now that we have investigated the possibility of this extrapolation, we chose to try the one of a constant SFR at z ≳ 10.

We performed an MCMC maximum likelihood sampling of the 4-parameter model of ρ_SFR (z) in Eq. (2) and add as a fifth parameter the value of SFR density at z > 10.4, our last data point corresponding to a redshift of 10.4. We refer to it as ρ_asympt and chose to use all observations cited in Sect. 2.1 as constraints. Final values of parameters a, b, c, and d are close to the ones from Sect. 4.1. We find that there is a strong correlation between ρ_asympt and τ, because of the direct integration in Eq. (8) and so expect higher values of the optical depth for high values of ρ_asympt. Yet, τ values are limited by Q_H II data points and they have more impact on the global scenario. Indeed, models where Q_H II equals 30% as soon as z = 10 are allowed, whereas it is closer to 20% at the same redshift when ρ_SFR is extrapolated. The correlation observed in our model parameters likelihood functions between ρ_asympt and τ had alreadybeen noticed by Robertson et al. (2015), as a correlation between τ and the averaged value of ρ_SFR for z > 10. A linear regression gives $${\langle {\rho }_{\text{SFR}}\rangle }_{z>10.4}=0.51\tau -0.026[M_{\odot }{\text{yr}}^{-1}{\text{Mpc}}^{-3}],$$(12)

with a correlation coefficient r = 0.98.

In this parametrisation, ρ_asympt can take very low values (down to 10⁻⁴ [M_⊙ yr⁻¹ Mpc⁻³]) meaning that reionisation sources are almost completely absent at z > 10. It also has an upper limit of 0.016 [M_⊙ yr⁻¹ Mpc⁻³]. This is close to the redshift-independent evolution of ρ_SFR (≃10^−1.5 [M_⊙ yr⁻¹ Mpc⁻³]) considered by Ishigaki et al. (2015) for z > 3 in order to reproduce ${\tau }_{\text{2014}}=0.091{}_{-0.014}^{+0.013}$ (Planck Collaboration XVI 2014), when usual decreasing models only gave them τ ≃ 0.05. We can compare Sect. 4.1 results with this upper limit in Fig. 2a. Despite the wide range of possible values for ρ_asympt, all results are consistent with our data and in particular, optical depths always remain in the 68% confidence interval of τ_Planck.

Fig. 8 Panels a: redshift evolution of QH II for various choices of the magnitude limit in luminosity data: brown for Mlim = −10, orange for Mlim = −13 and beige forMlim = −17. The light orange region represents the 68% confidence level for the worst case scenario, i.e. Mlim = −17. Left panel: escape fraction fixed to the values used by corresponding references. Right panel: escape fraction varying between zero and one. Panels b: UV luminosity density ρUV with redshift in logarithmic scale for three values of the magnitude limit: Mlim =−10 in the upper panel, Mlim = −13 in the middle panel and Mlim = −17 in the lower panel. Data points are from Ishigaki et al. (2015) or adapted from Robertson et al. (2015). Maximum likelihood parametrisations corresponding to Eq. (3) (continuous lines) are shown for fits using all observational constraints. The 68% confidence interval is represented as the light blue region. These results are compared with a model forced to reproduce results fromcorresponding references, drawn as the purple lines. Left panel: escape fraction fixed. Right panel: escape fraction allowed to vary between zero and one.

Fig. 9 Result of the MCMC analysis for the parametrisation described in Eq. (3) with fesc added as a fitting parameter. Here, all three sets of observational data were used as constraints. The contours show the 1-, 2-, and 3-σ confidence levels for log(ρz=8), fesc and the derived parameter τ. Colours correspond to the different values of the magnitude limit used in the analysis: purple for Mlim =−10, blue for Mlim = −13 and orange for Mlim = −17. These results can be compared to reionisation histories displayed in Fig. 8a and b.

5.3 How are f_esc, ṅ _ion, and ρ_SFR correlated?

We expect a correlation between the amplitude a of the star formation rate density parametrisation Eq. (2) and the escape fraction. Indeed, f_esc takes no part in the estimation of ρ_SFR but they both take part in the calculation of ṅ_ion in Eq. (1) and then in the integration of Q_H II in Eq. (4). Thus, they must be constrained by the same data, so that the parameter a can be a proxy for variations in the escape fraction value. To investigate this possible correlation, we plotted the distributions of a × f_esc for various sets of constraints and in different models: with (PAR) and without (CST) the escape fraction as a fifth fit parameter and with all constraints.We find that CST gives a lower value than PAR with a relative difference of 2.8%. This hints at a correlation between a and f_esc but more tests are needed to confirm or infirm this result.

To further investigate the link between f_esc, ṅ _ion and ρ_SFR, we considered values of the reionisation rate at various redshifts, used in Kuhlen & Faucher-Giguère (2012) and Robertson et al. (2013), and inferred from measurements and calculations of Faucher-Giguère et al. (2008); Prochaska et al. (2009); Songaila & Cowie (2010). We call NION the run using these new constraints – in addition to the others – and FREE the one skipping them, corresponding to ALL from Sect. 4.2.

We compare in Table 3 values of the reionisation rate at various redshifts for NION, FREE and Kuhlen & Faucher-Giguère (2012). NION gives results close to data points, increasing with z, whereas FREE values are significantly higher and decrease with redshift. This difference in the evolutions of ṅ _ion (z) is directly linked to f_esc (z). We see in Fig. 10 that f_esc,FREE is almost constant with redshift and therefore ṅ_ion (z) decreases on this redshift range: because f_esc values remains quite high, there is no need for many ionising sources at high redshift. On the contrary, when the constraints on ṅ_ion are included in the fit, the reionisation rate takes overall lower values (see Table 3) so that f_esc,NION has to take higher values at high redshift (saturating to 1 for z ≥ 15) to compensate for the lack of ionising sources. However, this is still hardly sufficient and we find that for NION, the reionisation process needs to start as early as at z = 18 to fully ionise the IGM, with Q_H II = 1 being reached later than others models, around z = 5.5. This behaviour leads to a high value of τ = 0.082, at the edge of the 3-σ confidence interval of τ_Planck and therefore, hardly compatible with observational results (Planck Collaboration Int. XLVII 2016). Removing the constraints on the filling factor, f_esc remains low on the whole redshift range (< 20). We then get values of the optical depth in agreement with Planck (0.058 ± 0.011) but reionisation does not end before z ~ 4. Thus, the estimations on the reionisation rate from Faucher-Giguère et al. (2008); Prochaska et al. (2009); Songaila & Cowie (2010) are compatible with one observable at a time: either the ionisation level – leading to a higher value of τ –, or the Thomson optical depth – so that reionisation ends around z ~ 4 –, but cannot match all observations in a coherent way.

6 Conclusions

We used the latest observational data available on reionisation history, i.e. cosmic star formation density, ionised fraction of the IGM and Thomson optical depth derived from Planck observations to find that they are all compatible with a simple and credible scenario where reionisation begins around z = 15 and ends by z = 6. Among all data, star formation history seems to be the most constraining for the EoR.

An investigation of various parametrisations of the escape fraction of ionising photons has lead us to conclude that it is very well constrained by observations: when considered constant with redshift, values allowed by the fit range from 20% to 28%; when considered redshift-dependent, from f_esc (z = 4) ≃ 17% to f_esc (z = 30) ≃ 26% following a low increase with z. The fiducial constant value of 20% often used in papers seems then to be perfectly consistent with our data. However, one must keep in mind that these results strongly depend on the hypothesis we make about the magnitude limit as a lower value of M_lim will require higher values of f_esc and vice versa. While the constraints on τ are unaffected by the assumption on M_lim, the confidence range on f_esc is enlarged for M_lim = 10. Furthermore, our different sets of observations seem to be in tension with each other for M_lim = −17 or for values of f_esc ≲ 10%.

On the contrary, the clumping factor of ionised hydrogen in the IGM can take a wide range of different values without impacting the reionisation observables significantly. For instance, when take C_H II as a redshift-independent parameter, its relative standard deviation is 41% whereas it isat most 7.6% for Q_H II (z)². The result is the same when we consider that C_H II depends on redshift: a great variety of possible evolutions gives the same scenario in terms of ionisation level. There is no greater impact on Thomson optical depth values, which vary of a maximum of a few percent compared to ${\langle \tau \rangle }_{C_{H_{II}}=3}$ and always remains in the 1-σ confidence interval of τ_Planck. Observational constraints are thus extremely robust to variations of the clumping factor. We nevertheless find a correlation between the averaged value of C_H II for z ∈ [6.8, 30] and τ: the linear fit $${\langle C_{H_{II}}\rangle }_{z>6.8}=-350\tau +24.4$$(13)

provides a good description of their connection³. This supports the use of a redshift-independent clumping factor to study the EoR. A possible choice, consistent with observations, would then be C_H II = 3, the fiducial value often used in papers, because it lies in the range of the ML C_H II values found in Sect. 4.3.

Last, a quick study on the possible reionisation sources at z ≳ 10 showed that there is no need for exotic sources such as early quasars (Madau & Haardt 2015) or for an artificial increase in star formation density at high redshift (Ishigaki et al. 2015). When their luminosity functions are extrapolated, a hypothesis still recently strongly supported by Livermore et al. (2017), star-forming galaxies provide enough photons to have a fully ionised IGM at z = 6.

Fig. 10 Evolution of fesc with redshiftwhen ṅion data points are used. ML models are shown for various set of constraints: blue when all constraints are used (NION); coral when ṅion constraints are skipped (FREE). Horizontal dashed lines represent the mean value of fesc over 4 ≤ z ≤ 30 for the model of the corresponding colour. Inferences are compared to results of Kuhlen & Faucher-Giguère (2012; KFG12, purple dashed line).

Acknowledgements

The authors thank B.E. Robertson for kindly providing us with his compilation of star formation rate densities. They thank the referee for useful comments. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013; 2018); matplotlib, a Python library for publication quality graphics (Hunter 2007) and emcee, an implementation of the affine invariant MCMC ensemble sampler (Foreman-Mackey et al. 2013). This work was partly supported by Programme National de Cosmologie et Galaxies (PNCG). AG acknowledges financial support from the European Research Council under ERC grant number 638743-FIRSTDAWN as well as from an STFC PhD studentship.

Appendix A MCMC multidimensional plots

We show in this appendix the additional triangle plots of the runs ALL corresponding to the studies with f_esc as additional free parameter (see Sect. 4.1), with C_HII as additional free parameter (see Sect. 4.3), and finally with both free (see Sect. 4.4).

Fig. A.1 Results of the MCMC analysis for the ALL case when fesc is added as a free parameter. The contours show the 1, 2, and 3σ confidence levels for a, b, c, d, fesc, and the derived parameter τ.

Fig. A.2 Results of the MCMC analysis for the ALL case when CH II is added as a free parameter. The contours show the 1, 2, and 3σ confidence levels for a, b, c, d, CH II, and the derived parameter τ.

Fig. A.3 Results of the MCMC analysis for the ALL case when both fesc and CHII are added as a free parameter. The contours show the 1, 2, and 3σ confidence levels for a, b, c, d, fesc, CH II, and the derived parameter τ.

References

All Tables

All Figures

	Fig. 1 Results of the MCMC analysis for the ALL case. The contours show the 1, 2, 3σ confidence levels for a, b, c, d and the derived parameter τ.
In the text

Fig. 2 Panel a: star formation rate density ρSFR with redshift. Data points are determined from infrared (plotted in red) or ultraviolet (in grey) luminosity densities (Sect. 3). Maximum likelihood parametrisations (continuous lines) are shown for various set of constraints: blue when all constraints are used; coral when only data on star formation are used; green when τ and reionisation history data are used. The 68% confidence interval on ρSFR (light blue region) is drawn for the blue model. We note that the interval, corresponding to statistical uncertainties, is very narrow. These inferences are compared with a model forced to reproduce results from Robertson et al. (2015), cited as R15 in the legend, drawn as the purple dotted line. The horizontal dashed-dotted line corresponds to the upper limit on a hypothetical constant value of ρSFR for z > 10.4 (Sect. 5.2). Panel b: ionised fraction of the IGM QH II with redshiftfor same models as (a). Details on the origin of data points are given in the legend and Sect. 3. Inferences are also compared with the two evolutions used in Planck Collaboration Int. XLVII (2016, cited as P16) to model the reionisation process: a redshift-symmetric hyperbolic tangent as the brown dashed-dotted line and a redshift-asymmetric power-law in black.

In the text

Fig. 3 Ionised fraction of the IGM QH II with redshiftwhen fesc is introduced as a parameter. Details on the origin of data points are given in the legend. ML models (continuous lines) are shown for various set of constraints: blue when all constraints are used, coral when QH II constraints are skipped, for which the 68% and 95% confidence intervals are drawn in salmon. These inferences are compared with a model forced to reproduce results from Robertson et al. (2015, R15, purple dotted line) and with the two evolutions used in Planck Collaboration Int. XLVII (2016, P16): redshift-symmetric as the dashed-dotted brown line and redshift-asymmetric in black.

In the text

	Fig. 4 Possible evolutions of fesc with redshift. ML models are shown for various set of constraints: blue when all constraints are used; coral when τ constraints are skipped; green when QH II constraints are skipped. Horizontal dotted lines represent the mean value of fesc over 4 ≤ z ≤ 30 for the model of the corresponding colour.
In the text

Fig. 5 Panel a: possible evolutions of CH II with redshift. ML models are shown for the two models of Sect. 2.2: blue for the first, coral for the second. Dotted horizontal lines correspond to the mean value of CH II (z) for z > 6.8, where outputs of the model are used in calculations, for the model of the corresponding colour. The vertical line is located at z = 6.8. Lines of lighter colours represent various outputs of the sampling of the corresponding model. Panel b: redshift evolution of QH II for the same models of CH II(z). Inferences are compared to a result with CH II(z) = 3 in purple dashed line, and to the theoretical models of Planck Collaboration Int. XLVII (2016): a redshift-symmetric model in black and a redshift-asymmetric model in brown. References. HM12: Haardt & Madau (2012) or Eq. (10). I07: Iliev et al. (2007) or Eq. (11).

In the text

	Fig. 6 MCMC distribution for fesc and CHII when both are taken as fit parameters (other parameters fixed). The escape fraction is allowed to vary between 0.1% and 100%, the clumping factor between one and seven.
In the text

	Fig. 7 Redshift evolution of QH II when all parameters (a, b, c, d, fesc, CH II) are free andall datasets used. Figure A.3 shows the corresponding constraints on assumed parameters.
In the text

Fig. 8 Panels a: redshift evolution of QH II for various choices of the magnitude limit in luminosity data: brown for Mlim = −10, orange for Mlim = −13 and beige forMlim = −17. The light orange region represents the 68% confidence level for the worst case scenario, i.e. Mlim = −17. Left panel: escape fraction fixed to the values used by corresponding references. Right panel: escape fraction varying between zero and one. Panels b: UV luminosity density ρUV with redshift in logarithmic scale for three values of the magnitude limit: Mlim =−10 in the upper panel, Mlim = −13 in the middle panel and Mlim = −17 in the lower panel. Data points are from Ishigaki et al. (2015) or adapted from Robertson et al. (2015). Maximum likelihood parametrisations corresponding to Eq. (3) (continuous lines) are shown for fits using all observational constraints. The 68% confidence interval is represented as the light blue region. These results are compared with a model forced to reproduce results fromcorresponding references, drawn as the purple lines. Left panel: escape fraction fixed. Right panel: escape fraction allowed to vary between zero and one.

In the text

Fig. 9 Result of the MCMC analysis for the parametrisation described in Eq. (3) with fesc added as a fitting parameter. Here, all three sets of observational data were used as constraints. The contours show the 1-, 2-, and 3-σ confidence levels for log(ρz=8), fesc and the derived parameter τ. Colours correspond to the different values of the magnitude limit used in the analysis: purple for Mlim =−10, blue for Mlim = −13 and orange for Mlim = −17. These results can be compared to reionisation histories displayed in Fig. 8a and b.

In the text

Fig. 10 Evolution of fesc with redshiftwhen ṅion data points are used. ML models are shown for various set of constraints: blue when all constraints are used (NION); coral when ṅion constraints are skipped (FREE). Horizontal dashed lines represent the mean value of fesc over 4 ≤ z ≤ 30 for the model of the corresponding colour. Inferences are compared to results of Kuhlen & Faucher-Giguère (2012; KFG12, purple dashed line).

In the text

	Fig. A.1 Results of the MCMC analysis for the ALL case when fesc is added as a free parameter. The contours show the 1, 2, and 3σ confidence levels for a, b, c, d, fesc, and the derived parameter τ.
In the text

	Fig. A.2 Results of the MCMC analysis for the ALL case when CH II is added as a free parameter. The contours show the 1, 2, and 3σ confidence levels for a, b, c, d, CH II, and the derived parameter τ.
In the text

	Fig. A.3 Results of the MCMC analysis for the ALL case when both fesc and CHII are added as a free parameter. The contours show the 1, 2, and 3σ confidence levels for a, b, c, d, fesc, CH II, and the derived parameter τ.
In the text