\section{Introduction}

Photodetectors are one of the main devices (together with logical gates and waveguides) that are essential for applications in integrated optics. Over the last 20 years, improvements in nanofabrication techniques have enabled the creation of new photodetector devices based on nanobelts, nanowires and nanotubes \cite{intro_1D_ZnO}. These structures have different characteristics to those obtained from thin films or bulk materials. One dimensional (1D) metal oxide structures have been of great interest to the scientific community over recent years. Due to the high surface to volume ratio, these components benefit from properties specific to nanoscaling. ZnO nanowires have already proven their usefulness in the realisation of multiple components, such as field effect transistors \cite{FET_ZnO}, gas detectors \cite{biblio_gas_sensor}, photodetectors \cite{fast_response_ZnO_water}, LEDs \cite{biblio_LED} or even solar cells \cite{biblio_solar_cells}. \\
Zinc oxide is a wide band-gap semiconductor (3.37 eV) and as such, ZnO nanowires (nw) have been shown to be very promising UV detectors thanks to their significant photoconductive gain, as high as $10^{8}$ \cite{gain_ZnO}. This feature gives them sufficient sensitivity to detect small amounts of light (potentially at the single photon level).
Photoconduction is the physical effect that will interest us during this study. Under optical illumination, electrical charges are liberated in the semiconductor, increasing the density of free carriers and therefore decreasing the resistivity of the device. \\
So far, almost all of these detectors consist of an $n$ doped ZnO nw. ZnO nanowire photodetectors are known to rely on surface states, usually with the absorption of $\text{O}_2$ molecules at the surface, which contribute to the existence of gain \cite{prades_PCC_explaination}.
Unfortunately, this is accompanied by a persistence of the photocurrent (PPC) which can last for several hours, preventing the device from going back to its initial off-state after being illuminated. As a result, the device can not be used at high repetition rates. In this work, we present a complete study of ZnO based nw photodetectors. It combines electrical measurements at different pressures, temperatures, and gaseous compositions. A study of the nw photoluminescence in air and in vacuum was also performed. We addressed the photocurrent persistence and demonstrate a way to circumvent it. This step is based on a change in the polarity of the gate voltage, which facilitates the transit of charges to the electrodes and stems the persistence of the photocurrent.

\section{Experimental details}

Our ZnO nanowires were grown on a $c$-plane sapphire substrate by MOCVD (metalorganic chemical vapour deposition) in a horizontal reactor operating at a reduced pressure of $65\,  \rm mbar$.
Oxygen and zinc were supplied by nitrous oxide ($8.6\, \text{mbar}$) and diethylzinc ($9.1\, \mu \text{bar}$) respectively, using heloium as a carrier gas. The hexagonal nanowires were grown at a temperature of $875$°C, and ended up forming a dense vertical forest of nanowires.
In order to create these photodetectors, ZnO nanowires were detached from the substrate by sonification in an ultrasonic bath with ethanol for $5\, \rm s$. Then, the nanowires in suspension were drop-casted on a silicon substrate ($p$ doped with an $\text{SiO}_2$ insulating layer of $300 \, \rm nm$ on top) \cite{carac_zno_iop}. \\
The electrical contacts ($\text{Ti/Au}$, of respectively $30\, \rm nm$ and $300\, \rm nm$) were fabricated using e-beam lithography on PMMA resist. For the alignement of the contacts to the nw, mechanical indents were drawn on the surface and then used as landmarks to localize the nanowires' positions \cite{Wei_Geng_fab}. Argon ion beam etching (IBE) was performed for $80\, \rm s$ at $5\, \rm mTorr$ and $50\, \rm W$ before the evaporation. The main purpose of this IBE step is to create oxygen vacancies in order to reduce an effective barrier at the metal/semiconductor interface. This step helps in the formation of good ohmic contacts, which are necessary for photoconduction with a gain higher than 1 \cite{ohmic_necessity}. Figure \ref{photo_PD19}a shows a SEM image of a ZnO nw contacted with Ti/Au electrodes. The average diameter of the nanowires is $250\, \rm{nm}$, their average length is $8\, \mu$m and the distance between the two electrodes is $2.5\, \mu$m. Figure \ref{photo_PD19}b is a schematic of the electrical contact setup. The gate voltage $V_{\rm g}$ is applied through the silicon substrate and the bias voltage $V_{\rm ds}$ is applied via the metallic electrodes.

\begin{figure}
    \centering
    \includegraphics[height=6cm]{PD19_photo.png}
    \includegraphics[height=6cm]{nanowire_PD_edit.png}
    \caption{\textbf{a.} SEM image of a single nanowire photodetector. Scale bar is 5 $\mu$m. \textbf{b.} Schematic of the device, the drain-source bias $V_{\rm ds}$ is applied using the Ti/Au electrodes. A gate voltage $V_{\rm g}$ is applied below the sample, by polarising the silicon substrate.}
    \label{photo_PD19}
\end{figure}

\section{Device Characterisation}

After nanofabrication, electrical properties were studied. Dark current measurements were performed at room temperature (approximately $300\,\rm K$ and $1\, \rm bar$) but in a pure O$_2$ atmosphere (see Figure \ref{IV_PD19}). In the transfer curve, we observe a drop of the dark current $I_{\rm d}$ over 7 decades below a gate voltage threshold around -$10\, \rm  V$, demonstrating a very low dark current capability and a very high ON/OFF ratio. Hysteresis is observed, the origin of which is attributed to carriers trapped at the nw surface, as discussed later. It is noteworthy that the hysteresis is highly dependent upon the exact experimental conditions (operating gas, pressure, speed of the voltage sweep). The nanowire behaves like a transistor with a threshold voltage around $V_{\rm g}= -10\,V$, below which the current no longer circulates in the nanowire.
The nanowire presents an almost linear behaviour for the $I_{\rm d}-V_{\rm ds}$ characteristic over a range of gate voltages $V_{\rm g} = [-20,20]\, \rm V$ (see inset Figure \ref{IV_PD19}). Using Ohm's law and a fit of these curves we can estimate the resistance of the device.
The resistances are respectively $36.5\, \rm{T}\Omega$, $1.25\, \rm{M}\Omega$ and $296\, \rm{k}\Omega$ for $V_g = -20\,\rm{V}$, 0 and $20\, \rm V$. These values are consistent with our previous measurements \cite{Wei_Geng_fab}.
The ZnO nanowire resistivity can also be estimated from $\rho = R \frac{\phi}{l_{\rm ds}}$, with $R$ the electrical resistance ($\Omega$), $\phi$ the nanowire cross section area ($\rm{m}^2$) and $l_{\rm ds}$ the distance between the two electrodes ($\rm m$).
The average nanowire diameter $D$ was determined from SEM images to be around $250\, \rm nm$. The nanowires are homogeneous in diameter and the length of the device $l_{\rm ds}$ is $2.5\, \mu \rm m$ (see Figure \ref{photo_PD19}).
The resistivities found are $3\,10^{8}\, \Omega . \rm cm$, $9.8\, \Omega .\rm  cm$, $2.32\, \Omega . \rm cm$  for respectively $V_{\rm g}=-20\,\rm V$, $V_{\rm g}=0\, \rm V$ and $V_{\rm g}=20\, \rm V$.
In the literature a wide range of values can be found (from $10^{-3}$ to $10^5 \, \Omega . \rm cm$
\cite{difficult_nw_characterisation}). It is difficult to make a comparison though as in most reported results, the value of the gate voltage is not specified or even controlled during the experiment.
Also, it is important to highlight that the contact resistances are not taken into account in these resistivity calculations. Consequently the nanowire resistivities discussed here are an overestimate. Nevertheless, $I_{\rm d}-V_{ \rm ds}$ measurements have been performed on a previous sample using a 4-point contact design (fabrication details can be found in \cite{these_wei}), where the contact resistance contribution can be separated from the intrinsic nanowire contribution. The intrinsic nanowire resistivity was then estimated to be around $\rho_{\rm nw}=0.32\, \Omega . \rm cm $ for $V_{\rm g}=0\, \rm V$. These values are consistent with the literature \cite{difficult_nw_characterisation}, where resistivities ranging from 2 to 12 $\Omega . \rm cm$ are found for similar nanowire diameters. To conclude, the contact resistances are the main contribution to our device total resistance. Therefore there is an opportunity to further improve the performance of our PD by reducing the contact resistances.

\begin{figure}
    \centering
    \includegraphics[height=6cm]{hysteresis_Vg_vs_I_d.png}
    \caption{Room temperature (dark) $I_{\rm d}-V_{\rm g}$ transfer characteristic  obtained at $V_{\rm ds} =1 \rm V$. The measurements are performed in 1 bar oxygen atmosphere. The flow of electrons through the nanowire can be modulated by the gate voltage. Inset: $I_{\rm d}-V_{\rm ds}$ curves for -20, 0 and 20 V gate voltage showing an almost ohmic behaviour.}
    \label{IV_PD19}
\end{figure}

Photoluminescence (PL) measurements were performed both at room temperature under vacuum and in ambient air. Measurements have been done using a frequency-tripled, mode-locked Ti:Sapphire laser as an excitation source ($\lambda = 266\, \rm nm$, pulse duration = $150\, \rm fs$, repetition rate = $76 \, \rm MHz$). The beam was focused on the sample with a $36\times$ reflecting objective (numerical aperture of 0.5), the photoluminescence was then collected through the same objective and directed to an Andor $30\, \rm cm$ Shamrock spectrometer (with a spectral resolution of $0.07\, \rm nm$) For an excitation at $266\, \rm nm$ (above the bandgap energy of ZnO), one can observe a luminescence maximum at $378\, \rm nm$. The position of this maximum is in agreement with the results for PL from ZnO \cite{PL_T_ambiante} ($383\, \rm nm$).
Figure \ref{comparaison_air_vide}a presents the PL intensity of a collection of ZnO nanowires for wavelengths between $\lambda = 350\, \rm nm$ and $\lambda = 420\, \rm nm$, in air and under vacuum. The PL shape is similar in air and under vacuum, however the intensity in vacuum is higher than in air, which indicates that non-radiative recombination is enhanced by the presence of $\text{O}_2$ (as discussed later). The PL intensity was also recorded at various excitation power densities (see Figure \ref{comparaison_air_vide}b). The fits were carried using the equation $I_{\rm PL}=a P_{\rm exc}^b+c$. The PL intensities do not increase strictly linearly with the excitation power density. This tends to indicate the presence of nonlinear effects, especially for the sample in air, since its curve follows $I_{\rm PL}^{\rm air} \propto P_{\rm exc}^{1.25}$ (compared to its counterpart in vacuum $I_{\rm PL}^{\rm vac} \propto P_{\rm exc}^{1.15}$). These nonlinear effects bypass the classical band to band recombinations and can imply bound-excitons, biexcitons, etc. \cite{example_PL_ZnO}.\\

\begin{figure}
    \centering
    \includegraphics[height=6cm]{PL_293K_in_Air.png}
    \includegraphics[height=6cm]{comparaison_air_vide_PL.png}
    \caption{\textbf{a.} Photoluminescence spectra of ZnO nanowires under vacuum and ambient air, at room temperature for an incident optical power of $1.5 \times 10^{6} \, \rm{W.m}^{-2}$ \textbf{b.} Nanowire photoluminescence intensity variation versus excitation optical power density, in air and in vacuum. The fits are carried out with equation $I_{\rm PL}=a P_{\rm exc}^b+c$. }
    \label{comparaison_air_vide}
\end{figure}

The nanowire photocurrent time response has been investigated in various atmospheres in order to investigate the persistence of the photocurrent. The study of photoconduction in O$_2$ is relevant because it has been shown \cite{prades_PCC_explaination} that oxygen leads to a reduction in the PPC. Illumination was provided by various fiber pigtailed LEDs, synchronised with a Keysight E5270 precision I(V) Analyser. The nanowires were probed in a cryostat under various temperatures and atmospheres. When exposed to photons of energy greater than the gap of ZnO, the photoresponse evolution drastically changed depending on the surrounding gas and pressure. \\

\begin{figure}
    \centering
    \includegraphics[height=6cm]{effet_pression.png}
    \includegraphics[height=6cm]{loglog_data_article.png}\\
    \includegraphics[height=5.2cm]{data_rise_fall_article.png}
    \caption{\textbf{a.} Nanowire photocurrent response under a pure oxygen atmosphere at $5\, \rm mbar$, $50\, \rm mbar$ and $750\, \rm mbar$ (respectively blue, red and orange lines). Illumination at $\lambda = 365\, \rm nm$ starts at t=0 s and stops at t=300 s.\textbf{b.} Same curve on a log scale. \textbf{c.} and \textbf{e.} Raw data for rise and fall transients respectively. \textbf{d.} and \textbf{f.} Normalized rise and fall transients respectively with their fitted curves.}
    \label{effet_oxygen}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[height=6.7cm]{ZnO_nw_in_air.png}
    \includegraphics[height=6.8cm]{schema_depletion_nw_O2_UV.png}
    \caption{Schematic of the influence of oxygen on the electrical properties of the nanowire. \textbf{a.} In the dark, oxygen molecules are absorbed on the surface and create potential barrier which reduce the density of free electrons near the surface. \textbf{b.} UV light contributes to the desorption of molecules previously trapped, electrons are released and so electrical conductivity is enhanced.}
    \label{schema_photoconduction}
\end{figure}

The nanowire photoresponse was first studied in a pure oxygen atmosphere (see Figure \ref{effet_oxygen}a and b). These graphs represent the photocurrent intensity versus time. The measurement begins in the dark with a dark current below $10^{-13}\, \rm{A}$. In the interval $t \in [0, 300\, \rm{s}]$ the sample is illuminated with UV light at $ \lambda = 365 \, \rm nm $. At t=$300\, \rm s$ the illumination is stopped and the photocurrent decay is measured up to t=$1800\, \rm s$.
The photocurrent transients have been recorded under 3 oxygen pressures ($5\, \rm mbar$, $50\, \rm mbar$ and $750\, \rm mbar$). The transients (rise and decay of the photocurrent) are slow (much slower than a second) and vary with the oxygen pressure. The lower the pressure, the higher the photocurrent, and the longer its decay time. \\
The relatively slow response of the nanowire may seem to contradict with expectations. Indeed in photoconductors, the photocurrent gain is due to the relative transit times of electrons and holes \cite{fundamentals_photonics}. In order to maintain neutrality, the fastest carrier (here an electron) is injected into the photoconductor from one electrode and transits to the other electrode, again and again, as long as the slower photogenerated carrier (here a hole) does not reach its collecting electrode. Therefore several carriers transit through the biasing circuit for each single absorbed photon, generating a photoconduction gain. In our case where holes are much slower than electrons, the gain is given by $ G \approx \frac{t_{\rm{p-transit}}}{t_{\rm{n-transit}}} $. The transit time of the holes from one electrode to the other is $t_{\rm{p-transit}}$. The transit time of electrons from one electrode to the other can be approximated by $ t_{\rm {n-transit}} = \frac{L^2}{\mu_n V_{\rm ds}} \approx 0.3 \, \rm ns $ and the gain of the detector was measured as $ G \sim 10^6 $ (see later in the article). As a result the photocurrent should reach its maximum after about $ t_{\rm{p-transit}} = G \times t_{\rm{n-transit}} \approx 0.3 \, \rm ms $. Our measurements show a time response 3 orders of magnitude slower. Another process is obviously responsible for the observed photocurrent transients.
A more detailed analysis of rise and fall transients appeared necessary, as illustrated Figure \ref{effet_oxygen}a-f where the raw and normalized rise and fall transients are presented. \\

 Considering the rise transients, we observed two distinct contributions to the photocurrent, each with its own characteristic time constant and saturation value. To the best of our knowledge, this is the first time observation of two distinct photocurrent responses has been reported and discussed. For the sake off analysis, the second rise transient and the photocurrent decay have been plotted on Figure \ref{effet_oxygen} c-f on a log scale and from their own time origin. \\

 A first fast transient occurs as soon as illumination starts, up to a first rise plateau. The fast rise has a response time of the order of the calculated one $ (0.3 \, \rm ms) $ and is too fast to be resolved using our equipment (limited to a few ms). It contributes almost completely to the value of the photocurrent gain. It is likely that this first and fast photo-response does not involve  reactions occurring at the surface of the nanowire. It is mainly due to the time required to reach an equilibrium of the free carrier concentrations, electrons and holes, under illumination. The carrier concentrations are determined by optical generation of photo-carrier pairs (electron-hole pairs), photo-carrier recombination, and drift-diffusion of electrons and holes towards the contacts. This first rise time depends on factors such as the optical power, the recombination rates, and the carrier mobilities, and is too short for our measurement setup. This is why only a first plateau is observed, with only a slight dependence on the oxygen pressure. \\

The second rise transient, which we will call secondary photocurrent evolution, is shown in Figures \ref{effet_oxygen}b and c. The secondary photocurrent contributes only weakly to the gain (only by a factor of about 5), and can take a very long time to reach a plateau. At low oxygen pressure $(5\, \rm mbar)$ the second plateau is not reached after $300\, \rm s$.
The evolution of the secondary photocurrent rise transient can be approximated by:

    \begin{equation}
 I = I_0 \left( 1+\frac{1}{\frac{1}{\left(t/{\tau}_{\rm {rise}}\right)^{\alpha}}+1/A}   \right)
 \label{eq_1}
    \end{equation}

where $ I_0 $ is the first (fast) photocurrent plateau, $ A = \frac{I_{\infty}}{I_0} -1 $ accounts for the maximum increase of the secondary photocurrent, and ${\tau}_{\rm rise}$ is the characteristic time of the second rise transient. Fitting the secondary rise transient using this expression gives the various constants given in Table \ref{rise_fall_time_table}.
For all transients $ \alpha $ is close to unity. Therefore, at the very beginning of this process and before saturation occurs, the secondary rise transients increase almost linearly with time. The time constant ${\tau}_{\rm rise}$ of the process responsible for this transient is of the order of 3 seconds, and is almost independent of the oxygen pressure. The saturation value of the photocurrent depends slightly on the oxygen pressure. The secondary gain $I_{\infty}/I_0$ is about 5. In comparison the gain of the fast first photocurrent transient ranges from 6 to $ 12\times 10^4 $ for this experiment (see Table \ref{rise_fall_time_table}a), contributing to most of the total gain. \\

The second transient is attributed to the modification of the nanowire surface. In the presence of oxygen molecules and under dark conditions, oxygen can be chemisorbed at the nanowire surface in the form ${\rm O}_{2}^{-}$ \cite{comini}.
As illustrated in Figure \ref{schema_photoconduction}, these molecules are responsible for a surface depletion, the free electron density $n$ is reduced to ensure neutrality of the nanowire, and consequently the nanowire conductivity decreases. Indeed, in dark conditions, neglecting surface level pinning, neutrality requires the free electron density $n$ to compensate exactly for the residual donor concentration: ${N}_{D}^{+}=n$. In the presence of chemisorbed ${\rm O}_{2}^{-}$ molecules, the averaged free electron concentration is reduced to maintain global neutrality following: $n = {N}_{D}^{+} - S_{O_{2}}\frac{4}{D}$ where $S_{O_{2}}$ is the surface density of chemisorbed ${\rm O}_{2}^{-}$  molecules and $\frac{4}{D}$ a coefficient to account for the surface to volume ratio of a nanowire slice. The free electrons are pushed away from the surface due to the negative charge of the ${\rm O}_{2}^{-}$ molecules, which creates a depletion region all around the nanowire surface. The nanowire conductivity is therefore reduced in the presence of chemisorbed ${\rm O}_{2}^{-}$ molecules. The application of a negative gate voltage further reduces the free electron concentration to a point where the dark current becomes lower than the measurement noise, well below $0.1 pA$. When illuminated, photogenerated holes are attracted to the surface, where they recombine with the excess electrons of the chemisorbed ${\rm O}_{2}^{-}$ molecules, contributing to the oxygen molecule desorption (in a neutral form). The surface depletion disappears and the nanowire conductivity increases, rising the photocurrent. This loss of chemisorbed oxygen can be partly counterbalanced by the chemisorption of new oxygen molecules, if available in the ambient atmosphere. The surface density of adsorbed oxygen molecules will reach an equilibrium due to these two competing processes. The lower the oxygen pressure, the deeper the surface oxygen desorption, the higher the electron density, the higher the conductivity, and the higher the secondary photocurrent plateau. The speed of the ${\rm O}_{2}^{-}$ desorption process does not depend on the surrounding oxygen pressure within the experimental range. A high surrounding oxygen pressure helps reaching low dark current conditions, however a low oxygen pressure ensures deeper ${\rm O}_{2}^{-}$ desorption and higher photocurrent gain. \\

 An examination of the photocurrent decay after illumination has been turned off reveals a clear dependence of the fall transients on the oxygen pressure: the higher the pressure the faster the photocurrent decay. The transient can be approximated by the following expression:
\begin{equation}
I = I_{\infty} \left( \frac{1}{\left( t/{\tau}_{\rm fall}\right) ^{\alpha}+1}   \right)
\label{eq_2}
\end{equation}

The parameters obtained for the three oxygen pressures are reported in Table \ref{rise_fall_time_table}b. Again $ \alpha $ is close to unity for all pressure conditions, indicating a photocurrent decaying closely to $ 1/t $. The decay process time constant ${\tau}_{\rm fall}$ is however very sensitive to the oxygen pressure, as already observed in the literature \cite{absorption_O2_pression}. The time constant ${\tau}_{\rm fall}$  increases from $ 35\, \rm s $ to $ 360\, \rm s $ when reducing the oxygen pressure from $(750\, \rm mbar)$ to $(5\, \rm mbar)$. The decay process involves a regeneration of the ZnO surface with fresh oxygen. When fresh oxygen is lacking due to external conditions the regeneration process becomes very slow or can even stop. \\

\begin{table}[]
\centering
\begin{tabular}{c|c|c|c|c|cccccl}
\cline{2-5} \cline{7-10}
$\textbf{a.}$ & $\textbf{Rise}$                                  & 5 mbar          & 50 mbar         & 750 mbar        & \multicolumn{1}{c|}{$\textbf{b.}$} & \multicolumn{1}{c|}{$\textbf{Fall}$}   & \multicolumn{1}{c|}{5 mbar} & \multicolumn{1}{c|}{50 mbar} & \multicolumn{1}{c|}{750 mbar} &  \\ \cline{2-5} \cline{7-10}
              & $\alpha$                                         & 0.8             & 1.05            & 0.9             & \multicolumn{1}{c|}{}              & \multicolumn{1}{c|}{$\alpha$}          & \multicolumn{1}{c|}{0.8}    & \multicolumn{1}{c|}{0.85}    & \multicolumn{1}{c|}{0.9}      &  \\ \cline{2-5} \cline{7-10}
              & $\tau_{\rm rise}$ (s)                                & 3.6             & 2.8             & 3.6             & \multicolumn{1}{c|}{}              & \multicolumn{1}{c|}{$\tau_{\rm fall}$ (s)} & \multicolumn{1}{c|}{360}    & \multicolumn{1}{c|}{100}     & \multicolumn{1}{c|}{35}       &  \\ \cline{2-5} \cline{7-10}
              & $I_{\infty}/I_0$, $2^{\text{nd}}$ transient gain & 6               & 4.5             & 4.6             &                                    &                                        &                             &                              &                               &  \\ \cline{2-5}
              & $1^{\text{st}}$ transient gain                   & $1.2\times 10^{5}$ & $1.1\times 10^{5}$ & $6.2\times 10^{4}$ &                                    &                                        &                             &                              &                               &  \\ \cline{2-5}
                            & Total gain $G$                   & $7.2\times 10^{5}$ & $5.0\times 10^{5}$ & $2.9\times 10^{5}$ &                                    &                                        &                             &                              &                               &  \\ \cline{2-5}
\end{tabular}
\caption{\textbf{a.} Extracted constants of rise \textbf{a.} and fall \textbf{b.} secondary photocurrent transients. All transients $\alpha$ values are close to 1, indicating an almost linear time dependence of all processes. The various rise constants are almost independent of the oxygen pressure. Conversely, the fall transient is faster by a factor of 10 when increasing of the oxygen pressure from $ 5 \, \rm mbar $ to $ 750 \, \rm mbar $. The photocurrent gains have been estimated using Equation \eqref{eq_Gain} as described in the text.}
\label{rise_fall_time_table}
\end{table}

Figure \ref{I_ph_effet_atm} presents similar experiments but for different gases at room pressure (argon, air, $\rm{O}_2$) and under vacuum (approximately $ 10^{-4}\, \rm mbar $). The initial dark currents depend on the atmosphere. The dark currents in air and oxygen are below $10^{-13} \, \rm{A}$, and are higher in vacuum and argon (respectively $14 \, \mu \rm{A}$ and $2 \, \mu \rm{A}$).
Data have been analysed using the methodology previously described. The transient phases (increase and decrease of the photocurrent $I_{\rm ph}$) appear rapid in the presence of oxygen and air, and the decay is even faster in air, suggesting that humidity plays some role in the oxygen chemisorption process, increasing the chemisorption rate. An increase of a factor of 100 is observed between the decay times in oxygen and vacuum. This confirms that the persistence of the photocurrent can be reduced in the presence of oxygen. For the rise processes, analysis show that similar time constants are obtained at ambient pressure regardless of the atmosphere, but the process seems slower in vacuum, unexpectedly. We believe that the extracted time constant in vacuum is not accurate since the photocurrent is already very high when the illumination starts. In vacuum the surface hardly regenerates its oxygen density and the photocurrent almost never returns to a comparable initial state than in other experiments. Actually in the present experiment the dark current under vacuum was about $ 14\, \mu \rm A $, about half of the saturation photocurrent.\\

From these various experiments the rise and fall transients of the photocurrents can be well understood. Ambient $\text{O}_2$ molecules are physisorbed on the surface of the nanowire, usually because of surface defects such as oxygen vacancies \cite{origin_doping_ZnO}, as observed in previous works \cite{prades_PCC_explaination}. Once physisorbed, oxygen molecules can trap electrons, becoming chemisorbed in the form $\rm{O}_{2}^{-}$ \cite{comini}, which increases the band bending at the surface to ensure neutrality. The band bending involves ionised donors, the presence of which is confirmed by conductivity measurements under various gate voltages, as shown above. The surface depletion shrinks the conductive channel (see Figure \ref{schema_photoconduction}) and reduces the number of free electrons. The nanowire conductivity tends to be reduced by the presence of chemisorbed oxygen at the nanowire surface. This is an efficient effect in nanowires because of their high surface to volume ratio. A similar process can arise from dipolar $\text{H}_2\text{O}$ molecules, which in the presence of $\text{O}_2$ can react to form two $\text{OH}^{-}$ groups \cite{TiO2_absorption_eau}. This might be one explanation for the very fast recovery process of the dark photocurrent in ambient air.\\

\begin{figure}
    \centering
    \includegraphics[height=6cm]{effet_atm.png}
    \includegraphics[height=6cm]{valeur_fit_all_atm.png}
    \caption{\textbf{a.} Photocurrent versus time under different atmospheres at room temperature. Illumination at $365\, \rm nm$ starts at $t=0\, \rm s$ and stops at $t=300\, \rm s$. The photocurrent persistence is much longer under vacuum and argon than in presence of oxygen. \textbf{b.} Values obtained from the fits of the rise and fall transients of the photocurrent using the formula \eqref{eq_1} and \eqref{eq_2} described in the text. }
    \label{I_ph_effet_atm}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[height=5.4cm]{I_ph_vs_Temperature_article_final.png}
    \caption{\textbf{a.} Photocurrent measurements at different temperatures ($100\,\rm K$-$293\, \rm K$), in oxygen atmosphere ($0.5\, \rm mbar$), illuminated at $365\, \rm nm$. First photocurrent rise plateau is indicated. \textbf{b.} Data for rise transient on a log scale. \textbf{c.} First photocurrent plateau versus $1000/T$: activation energy is determined from the formula $I_{\rm ph} = I_0 \exp{ \left( -E_{\rm a}/k_BT \right)}$.}
    \label{richardson}
\end{figure}

 Temperature measurements have been performed to better characterize the photoconduction mechanism.  The photocurrent response has been measured at $V_{\rm g} =-20\, \rm V$ (in order to reduce dark current) for temperatures ranging from $100\,\rm K$ to $293\, \rm K$, under low oxygen pressure of $0.5\, \rm mbar$. As previously we observed two successive increases of the photocurrent (see Figure \ref{richardson}a and b). A first plateau is reached almost instantly, then the current slowly increases towards a second plateau. Unfortunately the second plateau is not reached after $(300\, \rm{s})$ under a low oxygen pressure $(0.5\, \rm mbar)$. This can be understood since the second plateau involves oxygen desorption at the nanowire surface. The activation energy of this process is high, about $(0.7\, \rm eV)$ as determined by \cite{chemisorption_O2_ZnO}, therefore the process becomes very slow at low temperatures, as confirmed by our experiments. However, as shown Figure \ref{richardson}b, the first plateau is clearly reached and could be analysed. From a plot of the photocurrent plateau $I_{\rm ph}$ versus $1000/T$ an activation energy of $E_a = 115\, \rm meV$ is determined (see Figure \ref{richardson}c).

 This activation energy could be a signature of the surface equilibrium reached under dark conditions. The surface depletion at the nanowire surface, due to chemisorbed oxygen molecules, is responsible for the separation of photo-generated electrons and holes. The depletion width reduces the electron-hole recombination rate, and consequently influences the equilibrium concentrations of free electrons and holes. This equilibrium is responsible for the high photocurrent gain of nanowires compared to bulk photoconductors. The detailed analysis of the equilibrium reached in dark conditions, such as the surface depletion width and the chemisorbed oxygen surface density, is complex and beyond the scope of this paper. It would probably require complementary experiments. However the present activation energy determination could help in understanding the surface dynamics of ZnO nanowires.\\

In order to gather more information about the photoconduction in our component, the photoresponse has been studied as a function of the incident optical power density in an oxygen atmosphere. In Figure \ref{I_ph_vs_P}a, the maximum photocurrent is presented as a function of the incident optical power density for three nanowires. The maximum photocurrents are within the same order of magnitude for the three nanowires studied, and increase with incident optical power as expected. The photocurrent gain can be estimated from the following expression:
\begin{equation}
 G = \frac{I_{\rm ph}/q }{\eta _{\rm ext} \, P_{\rm exc}/(h \nu)}
 \label{eq_Gain}
 \end{equation}
where $q$ is the elementary charge, $\eta _{\rm ext}$ is the external quantum efficiency, $I_{\rm ph}$ is the photocurrent, $P_{\rm exc}$ is the incident optical power, and $h \nu$ is the photon energy. In a first approximation $\eta _{\rm ext}$ can be estimated considering a rectangular nanowire with a width and a thickness of $d_{\rm nw}=250\, \rm nm$ and a length $L=2.5\, \mu \rm m$. Photon absorption in the ZnO nanowire is estimated using the Beer-Lambert equation, with an absorption coefficient $ \alpha = 3 \times 10^{3}\,  \rm{cm}^{- 1} $ as measured previously \cite{Wei_absorption} for a nanowire of similar diameter. The calculated absorption is almost complete ($>92\%$). The reflection at the air/nanowire interface is approximated to $18\%$, as it would be for a plane wave under normal incidence. Finally, taking into account the illumination setup ($50\, \mu \rm m$ diameter cleaved multimode fiber) and the nanowire surface, the external quantum efficiency is evaluated to be  $\eta_{\rm ext} = 2.4\times 10^{-4}$. \\

From the measured photocurrent $I_{\rm ph}$, the gain can be calculated. Results are shown in Figure \ref{I_ph_vs_P}b. One notices for all the nanowires a similar decrease of the gain when increasing the incident optical power. The gain variation can be fitted by an inverse power law  $G = P^{-k}$, with $k=0.94$. The photocurrent gain is therefore almost inversely proportional to the incident optical power. This is in agreement with Soci \textit{et al.} \cite{gain_ZnO} where $k$ was found to be $0.7$. According to this work, the photocurrent can be expressed from the photon absorption rate F ($\rm{s}^{-1}$) as:
\begin{equation}
I_{\rm ph} = q \frac{\tau_l^{0}}{\tau_{t}} \frac{F}{1+(F/F_0)^k}
\label{eq_4}
\end{equation}
where $q$ is the elementary charge, $\tau_t$ the electron transit time, $\tau_l^0$ is the hole transit time (or the electron lifetime), and $F_0$ is a constant. The photon absorption rate is directly proportional to the incident optical power: $F \propto P_{\rm exc} $. Consequently the photocurrent varies like $I_{\rm ph} \propto P_{\rm exc}^{1-k} $ with the incident optical power. This is a counter-intuitive result when $k$ is close to 1, as $I_{\rm ph}$ becomes almost independent of the incident optical beam power. Our device is in this situation: as can be seen Figure \ref{I_ph_vs_P}a, the photocurrent is only multiplied by a factor 2 when the incident optical power $P_{\rm exc}$ is increased by more than 450. It means the sensitivity is increased at low optical power, and therefore, the device is suited for low intensity light detection.\\

\begin{figure}
    \centering
    \includegraphics[height=6.5cm]{I_ph_vs_power_at_365nm.png}
    \includegraphics[height=6.5cm]{gain.png}
    \caption{\textbf{a.} Maximum photocurrent at $293\, \rm K$ in 1bar $\rm{O}_2$ atmosphere, as a function of light excitation power density for three different nanowires. \textbf{b.} Nanowire photodetector gain as a function of the excitation power density for the three nanowires.}
    \label{I_ph_vs_P}
\end{figure}

The persistence of the photocurrent is detrimental to many applications of the device since it prevents high frequency operation. We discuss now a method to quickly restore the nanowire photodetector to its initial state once illumination has been turned off. As shown Figure \ref{I_ph_vs_t} (note the log scale), the photocurrent decay can be very long after illumination is turned off in normal conditions: it may take several hours. This prevents using the photodetector at useful rates. Luckily, an appropriate gate voltage pulse happens to stop the persistent photoconductivity and restore very low dark currents. The persistent photocurrent is due to a very slow return of the nanowire to its neutral state: after illumination, the surface is free of chemisorbed oxygen and the surface depletion is low. The conductivity of the nanowire is high is this state. A positive gate bias helps oxygen becoming trapped at the nanowire surface and restoring the surface depletion faster. The time required to get back to the initial state depends on the dark current requirements for a given application. For a dark current lower than $10^{-13}\, \rm A$, the time is shortened by the gate pulse from hours to a few minutes (see Figure \ref{I_ph_vs_t}), which considerably improves the device operability.  \\

\begin{figure}
    \centering
    \includegraphics[height=6cm]{retour_etat_initial_I_ph_vs_t.png}
    \caption{Nanowire response to illumination at $365\,$nm in an oxygen atmosphere and at room temperature (note the log scale). Return to the initial off-state is accelerated using an appropriate gate pulse.}
    \label{I_ph_vs_t}
\end{figure}

\section{Conclusion}

Single ZnO nanowire photodetectors have been fabricated, and their electrical properties have been investigated in the dark and under illumination. Experiments were performed under various atmospheres and in a temperature range of $100\, \rm K$--$293\, \rm K$. For the first time the existence of a two step photocurrent response was measured and discussed. A high photoconductive gain of $7.8\times 10^{7}$ was determined, associated with very low dark currents, well below $0.1\,\rm pA$, when the device is appropriately biased. These properties demonstrate the potential of nanowires for very high sensitivity photodetection.
The oxygen chemisorption process shows an activation energy of $E_{\rm a} = 115\, \rm meV$.
Finally, an appropriate gate pulse was demonstrated to shorten the persistence of the photocurrent, and restore the photodetector operating condition. This is an important step towards future application of these high gain photodetectors.
Further work is now underway in order to investigate higher rate$/$high gain photodetection operation.

\section{Acknowledgements}
The authors wish to thank the Dr V. Sallet and Dr A. Lusson from the GeMac (University of Versailles) for providing us the nanowires. We acknowledge the financial support of the STSM programme from the COST Action MP1403 'Nanoscale Quantum Optics'. e acknowledge the financial support of the  Labex ACTION (ANR-11-LABX-0001-01).
We acknowledge the financial support of the Grand Est region as all the fabrication was realised using the regional platform “Nanomat” under project “Nanogain”.

\bibliographystyle{iop-num}


\section{Introduction}
\label{sec:intro}
The ultra-violet (UV) radiation of the first populations of stars and accreting black holes reionized hydrogen in the intergalactic medium (IGM) in the Universe's first billion years. Understanding the reionization process enables us to learn about the properties of these first light sources \citep[e.g.,][]{Loeb2001,Robertson2010}.

However, to model reionization and interpret observations, significant assumptions are often made about the properties of these sources. Were they galaxies or quasars? What was their number density? How clustered were they? How hard were their ionizing spectra? A common approach, assuming reionization is driven by galaxies, is to parameterise the comoving density of ionizing photons, $\dot{N}_\textrm{ion}$, as a product of: the galaxy UV luminosity density, $\rho_L$; a conversion from observed UV photons (rest-frame $\sim 1500$\,\AA) to unobserved hydrogen-ionizing photons ($< 912$\,\AA), $\xi_\mathrm{ion}$; and the fraction of photons which escape galaxies to ionize the IGM, $f_\textrm{esc}$ \citep[e.g.,][]{Robertson2010}. However, in the reionization epoch, $z\simgt6$, $f_\textrm{esc}$ cannot be directly measured due to the high opacity of the IGM to ionizing photons \citep[e.g.,][]{Madau1995}, and the most straightforward way to estimate $\xi_\mathrm{ion}$ requires dust-corrected H$\alpha$ emission \citep{Leitherer1995} which calls for $\simgt4.6\,\mu$m spectroscopy e.g. with \textit{JWST}\ NIRSpec. Thus these parameters are often treated as constants in models \citep[e.g.,][]{Robertson2015,Bouwens2015}. One may ask how valid such assumptions are in the context of current constraints on reionization.

In the past year, several independent measurements provided evidence that reionization was a late ($z<9$), moderately extended ($\Delta z \simlt 4$), process. The optical depth to cosmic microwave background (CMB) photons provides an integral constraint: the \citet[][hereafter P18]{PlanckCollaboration2018} results suggest reionization's midpoint is at $z=7.7\pm0.7$. Observations of reduced Lyman-alpha (Ly$\alpha$) emission from high-redshift sources can be used to measure the neutral content of the IGM \citep[e.g.,][]{McQuinn2007a,Dijkstra2011,Mason2018}, as Ly$\alpha$\ photons are absorbed by neutral hydrogen. Damped Ly$\alpha$\ emission from galaxies and quasars currently provide the most precise constraints on the timeline of reionization, with most observations indicating a substantially neutral $z>7$ IGM \citep{Banados2017a,Davies2018,Mason2018,Mason2019c,Hoag2019a,Greig2019}.

The low number density of quasars observed at $z>6$ \citep[e.g.,][]{Parsa2018} has motivated recent work on galaxy-driven reionization. Fixing $f_\textrm{esc}=0.2$ and $\xi_\mathrm{ion} = 10^{25.2}$\,Hz/erg, \citet{Bouwens2015} and \citet{Robertson2015} found $z<10$ galaxies could reionize the IGM by $z\sim6$, providing undetected faint galaxies contribute significantly (extrapolating the galaxy luminosity function (LF) to $M_\textsc{uv} < -13$ with steep LF faint-end slopes, $\alpha \simlt -2$). However, these studies used parametric models, which could rule out possible evolutionary pathways \textit{a priori}: \citet{Bouwens2015} parameterised $\log_{10}{\dot{N}_\textrm{ion}} \propto z$; and \citet{Robertson2015} modelled the cosmic star formation rate (SFR) density in a parametric form based on $z\simlt8$ measurements.

Other works have explored models for $f_\textrm{esc}$ and $\xi_\mathrm{ion}$, allowing them to vary with galaxy properties. \citet{Finkelstein2019} used $f_\textrm{esc}$ as a function of halo mass, obtained in hydrodynamical simulations, where the lowest mass galaxies have the highest escape fractions, and modelled the galaxy UV luminosity density based on extrapolating $z\leq10$ UV LF fits. By shifting the load of reionization to the lowest mass galaxies \citet{Finkelstein2019} find a more extended reionization process and relatively flat $\dot{N}_\textrm{ion}(z)$, and also require a contribution from quasars at the end of reionization. \citet{Duncan2015} modelled the evolution of $f_\textrm{esc}$ and $\xi_\mathrm{ion}$ as functions of UV luminosity and spectral slope $\beta$ based on stellar population modelling, aiming to reproduce the $\dot{N}_\textrm{ion}(z)$ inferred by \citet{Bouwens2015}. \citet{Duncan2015} predict that faint blue galaxies produce and emit the most ionizing photons, which softens the requirement to extrapolate the LF in their work. \citet{Naidu2019}, based on the empirical galaxy evolution model of \citet{Tacchella2018a}, model $f_\textrm{esc}$ as a function of SFR surface density \citep[see also,][]{Sharma2017b,Seiler2019}, motivated by observations of Lyman Continuum leakers, and $\xi_\textrm{ion}$ using stellar population synthesis. With this model, \citet{Naidu2019} find $M_{\mathrm{UV}}<-18$ galaxies contribute $>80\%$ of the reionization budget \citep[see also,][]{Sharma2018}. However, \citet{Greig2017c} demonstrated that significant degeneracies exist between reionization parameters, which could not be broken with current reionization constraints, motivating an integrated approach to modelling the timeline and sources of reionization.

In this work, we seek to update our knowledge of $\dot{N}_\textrm{ion}(z\simgt5)$, by constraining it from recent estimates of the reionization timeline, and ask how consistent simple parametric models are with the allowed form of $\dot{N}_\textrm{ion}(z)$. This method enables us to estimate $\dot{N}_\textrm{ion}(z)$ at redshifts higher than is possible with the Ly$\alpha$\ forest, which becomes too heavily absorbed at $z>6$ \citep[e.g.,][]{Fan2006,Becker2013}. We aim to be maximally agnostic to the evolution of $\dot{N}_\textrm{ion}$ and so in a novel step we fit for a \textit{non-parametric} form of $\dot{N}_\textrm{ion}(z)$ allowing us to flexibly explore the allowed space. This offers an advantage over parametric models which rule out physically possible evolution \textit{a priori}.

The paper is structured as follows: Section~\ref{sec:reion} describes the reionization history model; Section~\ref{sec:nion} describes the method for inferring the redshift evolution of $\dot{N}_\textrm{ion}$; Section~\ref{sec:results} presents our inferred $\dot{N}_\textrm{ion}(z)$ and reionization history. We discuss our results in Section~\ref{sec:disc}, and summarise in Section~\ref{sec:conc}.

We use the \citet{PlanckCollaboration2015} cosmology: $(\Omega_\Lambda, \Omega_\textrm{m}, \Omega_\textrm{b}, n,  \sigma_8, H_0) =$ (0.69, 0.31, 0.048, 0.97, 0.81, 68 km s$^{-1}$\ Mpc$^{-1}$), as implemented in \verb|astropy|. Magnitudes are in the AB system. Distances, volumes, and densities are comoving. Confidence levels are 68\% unless otherwise stated.

\section{Modelling the reionization history}
\label{sec:reion}

Reionization progresses as ionizations overtake recombinations in the IGM. The reionization history: the volume-averaged ionized hydrogen fraction as a function of redshift, ${\overline{x}_\textrm{{\textsc{hii}}}}(z)$, can be calculated by solving the following differential equation ~\citep[e.g.,][]{Madau1999}:
\begin{equation}  \label{eqn:reion_Q}
	\frac{\mathrm{d}{\overline{x}_\textrm{{\textsc{hii}}}}}{\mathrm{d}t} = \frac{\dot{N}_\textrm{ion}}{\langle n_\textrm{H}\rangle} - \frac{{\overline{x}_\textrm{{\textsc{hii}}}}}{t_\textrm{rec}}
\end{equation}
where $\dot{N}_\textrm{ion}$ is the number density of ionizing photons in Mpc$^{-3}$ s$^{-1}$ and $\langle n_\textrm{H}\rangle$ is the mean number density of hydrogen atoms. The recombination time of the IGM is
$t_\textrm{rec}(z) = \left[C\alpha_B(T)n_e(1+z)^3 \right]^{-1}$,
where $\alpha_B(T)$ is the case B hydrogen recombination (i.e. opaque IGM) coefficient, $n_e$ is the number density of electrons (assuming singly ionized He), and $C=\langle n_H^2 \rangle/\langle n_\textrm{H}\rangle^2$ is the `clumping factor' which accounts for inhomogeneity in the IGM \citep{Madau1999}. For computational efficiency we fix $C=3$, motivated by IGM modelling and simulations \citep[e.g.,][]{Shull2012,Finlator2012,Kaurov2015}.  Appendix~\ref{app:clumping} discusses how our results are insensitive to this assumption. Throughout most of this work we use the \textit{neutral} fraction ${\overline{x}_\textrm{{\textsc{hi}}}} = 1 - {\overline{x}_\textrm{{\textsc{hii}}}}$.

Assuming galaxies produce the bulk of ionizing photons during reionization $\dot{N}_\textrm{ion}$ can be expressed as:
\begin{equation} \label{eqn:reion_nion_G}
    \dot{N}_\textrm{ion,G} = f_\textrm{esc} \xi_\textrm{ion} \rho_L
\end{equation}
where $f_\textrm{esc}$ is the fraction of ionizing photons which escape galaxies to the IGM; $\xi_\textrm{ion}$ is the production rate of ionizing photons per UV luminosity, in Hz/ergs, which depends on the stellar populations' initial mass function, metallicity, age and dust content; and $\rho_L$ is the dust-corrected UV luminosity density \citep[e.g.,][]{Robertson2010}. These parameters are likely functions of at least mass and redshift, but $f_\textrm{esc}$ and $\xi_\textrm{ion}$ are commonly treated as constant for simplicity.

The ionizing photon production rate from quasars can be derived from their spectral energy distribution (SED):
\begin{equation} \label{eqn:reion_nion_Q}
    \dot{N}_\textrm{ion,Q} = f_\textrm{esc,Q} \int d\nu \, \frac{\epsilon_{\nu}}{h\nu}
\end{equation}
where the integral limits are from $1-4$\,Ry (the energy of hydrogen-ionizing photons) and $f_\textrm{esc,Q}$ is the ionizing photon escape fraction from quasars (usually assumed $\sim1$). The UV SED of quasars is assumed to follow a double-power law, $\epsilon_\nu \propto \nu^\alpha$, with a pivot at 912\,\AA. For $<912$\,\AA\ we use $\alpha=-1.70$ \citep{Madau2015,Kulkarni2018}. This value was derived from $M_\textrm{UV} \sim -27$ $z\sim2$ quasars \citep{Lusso2015}. $\alpha$ may be steeper for lower luminosity quasars \citep[e.g., -0.56 -- -1.4,][]{Stevans2014,Scott2004}, which could increase the normalisation of $\dot{N}_\textrm{ion}$ by up to $\sim0.5$ dex.

\section{Non-parametric inference of $\dot{N}_\textrm{ion}$}
\label{sec:nion}

To be maximally agnostic about the form of $\dot{N}_\textrm{ion}$ we fit for it as a non-parametric function of redshift, i.e. $\dot{N}_\textrm{ion}(z)$ can take \textit{any} value at any redshift. We use redshift bins of width $\Delta z =1$, assuming smooth evolution on that scale, and fit for $\dot{\mathbf{N}}_\textbf{ion}(\mathbf{z})$ at $\mathbf{z}= \begin{bmatrix} 4 \; \dots \; 14 \end{bmatrix}$, i.e. 11 parameters. In the following we use $\dot{N}_\textrm{ion}$ to refer to the ionizing emissivity in general and the vector notation $\dot{\mathbf{N}}_\textbf{ion}$ to refer specifically to our inferred values of $\dot{N}_\textrm{ion}$ in each redshift bin. Using Bayes' theorem, the posterior probability for $\dot{\mathbf{N}}_\textbf{ion}(\mathbf{z})$ is:
\begin{equation} \label{eqn:nion_post}
p(\dot{\mathbf{N}}_\textbf{ion} \,|\, \textrm{data}) = \mathcal{L}(\textrm{data} \,|\, \dot{\mathbf{N}}_\textbf{ion})\, p(\dot{\mathbf{N}}_\textbf{ion})
\end{equation}
where $\mathcal{L}(\textrm{data} \,|\, \dot{\mathbf{N}}_\textbf{ion})$ is the likelihood of obtaining observed data given our model $\dot{\mathbf{N}}_\textbf{ion}$ (described in Section~\ref{sec:nion_like}). We set the prior, $p(\dot{\mathbf{N}}_\textbf{ion})$, such that:
\begin{itemize}
    \item $\log_{10}\dot{\mathbf{N}}_\textbf{ion}(z=5)$ is uniformly distributed between 49 and 53 ($\dot{N}_\textrm{ion}$ in units of Mpc$^{-3}$ s$^{-1}$), motivated by measurements of $\dot{N}_\textrm{ion}$ at $z\simlt5$ \citep{Becker2013}.
    %
    \item The gradient between redshift steps, $d\log_{10}\dot{\mathbf{N}}_\textbf{ion}/dz$, is uniformly distributed between -1 and 1 dex. This is motivated by the SFR/luminosity density of galaxies and quasars which fall by $<1$ dex per redshift at these redshifts \citep[e.g.,][]{Bouwens2015b,Finkelstein2015a,Oesch2018,Kulkarni2018}, even if bright galaxies dominate reionization \citep{Sharma2018}. The luminosity density of observable galaxies falls only $\sim0.2$ dex per redshift at these redshift, thus we do not expect enormous jumps in $\dot{N}_\textrm{ion}(z)$.
\end{itemize}
Figure~\ref{fig:nion} shows the posterior is not restricted by this prior.

In each likelihood call we sample $\dot{\mathbf{N}}_\textbf{ion}$, linearly interpolate it over a redshift grid $\Delta z=0.2$ and solve the reionization history Equation~\ref{eqn:reion_Q}.
While e.g. $f_\textrm{esc}$ in individual galaxies could fluctuate on these timescales due to supernova feedback \citep{Trebitsch2017}, the average over the ensemble galaxy population should be smooth. The likelihood for the data, given $\dot{\mathbf{N}}_\textbf{ion}$, is described in the next section.

To obtain the posterior distribution for $\dot{N}_\textrm{ion}(z)$ we use dynamic nested sampling via \verb|dynesty|\footnote{\url{https://dynesty.readthedocs.io/}} \citep{Speagle2019}, with the sampling settings optimised for posterior estimation.

\subsection{Measurements used in $\dot{N}_\textrm{ion}$ likelihood}
\label{sec:nion_like}

We use two sets of reionization constraints to create two likelihood functions. The total likelihood is the product of the individual likelihoods $\mathcal{L}_i$ for the below data given our model $\dot{\mathbf{N}}_\textbf{ion}$. Following \citet{Greig2017c}, we define a `Gold' sample of model-independent measurements, and an additional sample of more model-dependent measurements (described below). The `Gold' sample constraints are:

\begin{enumerate}
    \item \textbf{CMB electron scattering optical depth.} CMB photons scatter off a fraction $\tau_\mathrm{es}$ of free electrons created during reionization, suppressing CMB anisotropies by $\exp(-\tau_\mathrm{es})$, below angular scales corresponding to the size of the cosmological horizon at reionization.
    %
    \begin{equation}  \label{eqn:reion_tau}
    	\tau_\mathrm{es}(z) = \int_0^z dz' \; \sigma_\textsc{t} n_e (1+z')^2 \, {\overline{x}_\textrm{{\textsc{hii}}}}(z') \frac{c}{H(z')}
    \end{equation}
    %
    where $c$ is the speed of light, $\sigma_\textsc{t}$ is the Thomson scattering cross section, $n_e$ is the comoving number of free electrons \citep[assuming doubly-ionized Helium at $z<4$, following][the exact timing of Helium reionization has a negligible impact on our results]{Kuhlen2012,Bouwens2015}, and $H(z)$ is the Hubble parameter. As $\tau_\mathrm{es}$ is measured at $z\sim1100$ for the CMB it provides only an integral constraint on reionization and cannot give precise information about its timing.

    In each call of our likelihood we solve for $\tau_\mathrm{es}(\dot{\mathbf{N}}_\textbf{ion}, z=15)$ (assuming fully neutral IGM) and calculate the likelihood for obtaining the \citetalias{PlanckCollaboration2018} optical depth: $\tau_\mathrm{es} = 0.054\pm 0.007$, from our model, assuming Gaussian uncertainties.

    \item \textbf{Ly$\alpha$\ forest dark fraction.} The fraction of dark pixels in the Ly$\alpha$\ and Ly$\beta$ forest of $z\sim6$ quasars (hereafter, `dark fraction') provides a model-independent constraint on ${\overline{x}_\textrm{{\textsc{hi}}}}$, as neutral hydrogen in the IGM produces fully saturated absorption in the Ly$\alpha$\ and Ly$\beta$ forest \citep{Mesinger2010}. This is degenerate with absorption by self-shielded neutral gas post-reionization, so measurements provide upper limits. We use dark fraction measurements by \citet{McGreer2015}: ${\overline{x}_\textrm{{\textsc{hi}}}}(z=5.6) < 0.04+0.05 \; (1\sigma)$, ${\overline{x}_\textrm{{\textsc{hi}}}}(z=5.9) < 0.06+0.05$ and ${\overline{x}_\textrm{{\textsc{hi}}}}(z=6.1) < 0.38+0.20$. Following \citet{Greig2017c} we implement these via a likelihood with a uniform probability for ${\overline{x}_\textrm{{\textsc{hi}}}} < {\overline{x}_\textrm{{\textsc{hi}}}}_\mathrm{lim}$ and a half-Gaussian distribution with $\mu={\overline{x}_\textrm{{\textsc{hi}}}}_\mathrm{,lim}$ and $\sigma=\sigma_\mathrm{lim}$ for higher values of ${\overline{x}_\textrm{{\textsc{hi}}}}$.
\end{enumerate}

The second likelihood adds the following measurements of ${\overline{x}_\textrm{{\textsc{hi}}}}(z)$ from Ly$\alpha$\ damping of high-redshift sources, which account for patchy reionization in statistically robust methods. ${\overline{x}_\textrm{{\textsc{hi}}}}$ can be measured from an observed reduction in Ly$\alpha$\ as neutral hydrogen in the IGM absorbs Ly$\alpha$\ photons. However, such measurements rely on modelling the  sources observed as backlights to the neutral IGM:
\begin{enumerate}
    \item \textbf{Ly$\alpha$\ equivalent width (EW) distribution.} We use ${\overline{x}_\textrm{{\textsc{hi}}}}$ constraints at $z\sim7,7.5$ and 8 from \citet{Mason2018,Hoag2019a,Mason2019c}. These measurements infer ${\overline{x}_\textrm{{\textsc{hi}}}}$ from the redshift evolution of the Ly$\alpha$\ EW distribution in Lyman-break galaxies via comparison with realistic IGM simulations. As these constraints were relative to $z=6$ (i.e. assuming the IGM was fully ionized at $z=6$) we use the relative measurement $\Delta_{6-z}{\overline{x}_\textrm{{\textsc{hi}}}} = {\overline{x}_\textrm{{\textsc{hi}}}}(z) - {\overline{x}_\textrm{{\textsc{hi}}}}(6)$ in our likelihood. We use the posteriors $p({\overline{x}_\textrm{{\textsc{hi}}}}(z_\mathrm{meas}))$ from those works to directly calculate the likelihood $p(\Delta_{6-z}{\overline{x}_\textrm{{\textsc{hi}}}} \,|\, \dot{\mathbf{N}}_\textbf{ion})$ from our model. To account for the redshift uncertainty in these measurements due to using photometric samples ($\Delta z \sim 1$), we draw 1000 redshifts from the photometric redshift distribution for sources used in each measurement, maintaining detailed balance by using same drawn redshifts in every likelihood calculation, and use the median of the likelihoods calculated for the drawn redshifts.
    %
    \item \textbf{Ly$\alpha$\ emitter clustering.} Reionization increases the clustering of Ly$\alpha$\ emitting galaxies. We use the constraint of ${\overline{x}_\textrm{{\textsc{hi}}}}(z=6.6) \leq 0.5\,(1\sigma)$ by \citet{Sobacchi2015} obtained from comparing the angular correlation function of Ly$\alpha$\ emitters by \citet{Ouchi2010} to reionization simulations. We implement the likelihood as half-Gaussian distribution with $\mu=0$ and $\sigma=0.5$.
    %
    \item \textbf{Ly$\alpha$\ damping wings in quasar spectra.} The neutral IGM creates smooth Ly$\alpha$\ damping absorption in quasar spectra. We use recent constraints on ${\overline{x}_\textrm{{\textsc{hi}}}}$ from two $z>7$ quasars by \citet{Davies2018} and \citet{Greig2019}. We use the ${\overline{x}_\textrm{{\textsc{hi}}}}$ posteriors from those works to calculate likelihoods for our model \citep[using the `Intermediate HII regions' IGM simulation from][]{Greig2019}. We note these works obtained slightly different constraints (though consistent within $1-1.5\sigma$) primarily due to differences in modelling the intrinsic quasar spectra. We include both constraints to `marginalise' over differences in modelling.
\end{enumerate}

These Ly$\alpha$\ damping ${\overline{x}_\textrm{{\textsc{hi}}}}$ measurements are all obtained by comparison to large-scale reionizing IGM simulations. By construction these rely on parametric `sub-grid' models for the properties of reionizing sources, which set the timing and morphology of reionization. However, currently these ${\overline{x}_\textrm{{\textsc{hi}}}}$ measurements are relatively insensitive to the properties of reionizing sources: there is little difference between ${\overline{x}_\textrm{{\textsc{hi}}}}$ results obtained using different simulation setups \citep{Sobacchi2015,Mason2018,Greig2019}, meaning our results should be robust to the simulation assumptions. This is due to the structure of the reionizing IGM depending primarily on ${\overline{x}_\textrm{{\textsc{hi}}}}$, with a lesser some dependence on the clustering scale of the primary reionizing sources, and only weak dependence on the redshift of reionization \citep[as the luminosity-weighted power spectra of galaxies does not change much $z\sim6-10$,][]{McQuinn2007}. The impact of reionization morphology is mostly diluted in the measurements from galaxies, which span a redshift range greater than typical size of ionizing bubbles ($<\Delta z = 0.1$), and in the case of quasars due to the large uncertainties in ${\overline{x}_\textrm{{\textsc{hi}}}}$ measurements from single objects.

Additional constraints on reionization come from the patchy kinetic Sunyaev-Zel'dovich (kSZ) effect, which requires large 3D simulations to model, and is beyond the scope of this work. \citet{Greig2017c} demonstrate current patchy kSZ results \citep{George2015} favour late reionization, but do not provide more constraining power on ${\overline{x}_\textrm{{\textsc{hi}}}}(z)$ than $\tau_\mathrm{es}$ and the dark fraction.

\section{Results}
\label{sec:results}

\begin{figure*}
    \includegraphics[width=\textwidth]{figure1.pdf}
    \vspace{-16pt}
    \caption{Redshift evolution of $\dot{N}_\textrm{ion}$.
    %
    \textbf{(a)}
    1000 draws from our posteriors for $\dot{\mathbf{N}}_\textbf{ion}$. Dark grey lines show $\dot{\mathbf{N}}_\textbf{ion}$ inferred using only the \citetalias{PlanckCollaboration2018} optical depth, $\tau_\mathrm{es}$, and dark fraction constraints. dark purple lines show $\dot{\mathbf{N}}_\textbf{ion}$ inferred using the additional ${\overline{x}_\textrm{{\textsc{hi}}}}$ constraints from Ly$\alpha$\ damping. Light grey lines show 500 draws from the prior. The posteriors in both cases are substantially narrower than the prior, especially $5 < z < 10$, demonstrating that we have updated our knowledge on $\dot{\mathbf{N}}_\textbf{ion}$.
    %
    \textbf{(b)}
    The distribution of gradients in $\log_{10}\dot{N}_\textrm{ion}$ between $z=6$ and $z=8$. Colours same as above, shading show 68\% confidence region. Prior distribution shown as light grey cross-hatch. The grey diagonal hatch region shows range of the galaxies model described in panel (c). Orange points show models by \citet{Bouwens2015}, \citet{Finkelstein2019} and \citet{Naidu2019}.
    %
    \textbf{(c)}
    Same as (a) but now shaded regions shows the 68\% confidence region for our inferred $\dot{\mathbf{N}}_\textbf{ion}$. Circular points show the $\dot{N}_\textrm{ion}$ measurements by \citet{Becker2013}. We also plot the range of $\dot{N}_\textrm{ion}(z)$ allowed by the evolving number density of galaxies \citep[with constant $f_\textrm{esc}$ and $\xi_\mathrm{ion}$, integrating the UV LF down to $M_\textsc{uv} < -18, -15, -13$,][dashed lines]{Mason2015} and quasars \citep[constant $f_\textrm{esc}$ and EUV slope,][circular hatch]{Kulkarni2018} -- described in Section~\ref{sec:results_nion}.
    %
    \textbf{(d)}
    Same as (c), comparing $\dot{N}_\textrm{ion}$ models from \citet[][cross hatch]{Bouwens2015}, \citet[][light green region]{Finkelstein2019} and \citet[][orange dashed line]{Naidu2019}.}

    \label{fig:nion}
\end{figure*}

\subsection{Evolution of $\dot{N}_\textrm{ion}(z)$}
\label{sec:results_nion}

Figure~\ref{fig:nion} shows $\dot{\mathbf{N}}_\textbf{ion}(\mathbf{z})$ inferred from the measurements described in Section~\ref{sec:nion_like}. Panel (a) shows draws from our posteriors and the prior, to demonstrate how $\dot{N}_\textrm{ion}$ evolution is constrained when the likelihood is included (Section~\ref{sec:nion_like}). In both cases the inferred $\dot{\mathbf{N}}_\textbf{ion}$ is a smooth, declining function of redshift, with the tightest posteriors obtained at $z\sim6-8$ where the majority of the reionization history constraints are found (Figure~\ref{fig:history}). When the Ly$\alpha$\ damping ${\overline{x}_\textrm{{\textsc{hi}}}}$ constraints are included, the posterior for $\dot{\mathbf{N}}_\textbf{ion}$ is noticeably tighter and prefers a more rapid decline with redshift.

This is clearer in panel (b) which shows the distribution of derivatives of $\log_{10}\dot{\mathbf{N}}_\textbf{ion}$ between $z=6\rightarrow8$ from our $\dot{\mathbf{N}}_\textbf{ion}$ posteriors. $d\log_{10}\dot{\mathbf{N}}_\textbf{ion}/dz|_{z=6\rightarrow8} = -0.31\pm0.35$ inferred from the $\tau_\mathrm{es}$ and dark fraction constraints, and $-0.44\pm0.22$ when the Ly$\alpha$\ damping constraints are added. We also compare to models by \citet{Bouwens2015}, \citet{Finkelstein2019} and \citet{Naidu2019}, also shown in panel (d) and described in more detail below.

Panel (c) compares our inferred $\dot{\mathbf{N}}_\textbf{ion}$ to simple `Galaxies' and `Quasars' models. For the `Galaxies' model we take the galaxy UV luminosity density from the \citet{Mason2015} UV luminosity function (LF) model \citep[consistent with observations from $0 < z < 10$, see also,][]{Tacchella2013}, and calculate $\dot{N}_\textrm{ion,G}$ (Equation~\ref{eqn:reion_nion_G}), fixing $f_\textrm{esc} = 0.2$ and $\log_{10}(\xi_\mathrm{ion}) = 25.2$ for comparison with previous works \citep{Bouwens2015,Robertson2015}. To encompass uncertainty in the faint-end cut-off of the galaxy LF we show the range of $\dot{N}_\textrm{ion,G}$ allowed by integrating the LFs to $M_\textsc{uv} < -13$ \citep[comparable to the atomic cooling limit at $z\sim8-10$, $M_h\simgt10^8\,M_\odot$,][resulting in a flatter luminosity density and $\dot{N}_\textrm{ion}$]{Bromm2011}, $M_\textsc{uv} < -15$ \citep[the depth reached in deep HST imaging of galaxy cluster lensed fields, e.g.,][]{Livermore2017,Bouwens2017a} and $M_\textsc{uv} < -18$. For the quasar model we take the quasar ionizing emissivity at 912\,\AA\ from a homogenised compilation by \citet{Kulkarni2018} and assume a power-law quasar UV SED to calculate the hydrogen-ionizing emissivity (Equation~\ref{eqn:reion_nion_Q}). We assume the ionizing photon escape fraction from quasars, $f_\textrm{esc,Q}=1$, though this may be optimistic \citep[e.g.,][]{Micheva2017}, thus our model shows the upper limit of the quasar contribution. To encompass the uncertainty in the quasar LF faint-end cut-off we show the range of $\dot{N}_\textrm{ion,Q}$ obtained by integrating to $M_\textsc{uv} < -18$ or $M_\textsc{uv} < -21$. The galaxy model is consistent with $\dot{\mathbf{N}}_\textbf{ion}$ inferred from $\tau_\mathrm{es}$ and the dark fraction constraints. However, including the Ly$\alpha$\ damping constraints requires a steeper evolution in $\dot{N}_\textrm{ion}$ between $z\sim6-8$, as described above and shown in panel (b). The majority of the posterior inferred for the slope $d\log_{10}\dot{\mathbf{N}}_\textbf{ion}/dz|_{z=6\rightarrow8}$ favours steeper values than the galaxy luminosity density: 72\% of the posterior has $d\log_{10}\dot{\mathbf{N}}_\textbf{ion}/dz|_{z=6\rightarrow8} < -0.31$\,dex (steeper than the slope of the galaxy luminosity density if $M_\textsc{uv} < -18$),  80\% of the posterior has $d\log_{10}\dot{\mathbf{N}}_\textbf{ion}/dz|_{z=6\rightarrow8} < -0.25$\,dex, steeper than the galaxy luminosity density for $M_\textsc{uv} < -15$, and 83\% has $d\log_{10}\dot{\mathbf{N}}_\textbf{ion}/dz|_{z=6\rightarrow8} < -0.23$\,dex (steeper than the galaxy luminosity density for $M_\textsc{uv} < -12$). The inferred steep $\dot{N}_\textrm{ion}(z)$ gradient thus makes a model with constant $f_\textrm{esc}$ and $\xi_\mathrm{ion}$ and/or abundant faint galaxies less likely.

Panel (d) compares our inferred $\dot{\mathbf{N}}_\textbf{ion}$ to models by \citet{Bouwens2015}, \citet{Finkelstein2019} and \citet{Naidu2019}. The \citet{Bouwens2015} model was derived using older constraints on the reionization history, and used tighter priors on $\dot{N}_\textrm{ion}$, but is mostly consistent with our $\dot{\mathbf{N}}_\textbf{ion}$ inferred from the \citetalias{PlanckCollaboration2018} $\tau_\mathrm{es}$ and dark fraction. This model is similar to the `Galaxies' model in panel (c), with its evolution consistent with the galaxy luminosity density, dropping $\sim0.2$ dex per redshift. It is difficult for such a model, where $f_\textrm{esc}$ and $\xi_\mathrm{ion}$ are constant, to match the $\dot{N}_\textrm{ion}(z)$ evolution inferred from the full set of ${\overline{x}_\textrm{{\textsc{hi}}}}$ constraints which show later, more rapid, reionization -- $d\log_{10}\dot{\mathbf{N}}_\textbf{ion}/dz|_{z=6\rightarrow8} = -0.44\pm0.22$ dex.

The models by \citet{Finkelstein2019} and \citet{Naidu2019} model $f_\textrm{esc}$ as functions: of halo mass (with the highest $f_\textrm{esc}$ in low mass halos); and SFR surface density \citep[highest $f_\textrm{esc}$ in high SFR surface density objects, see also][]{Sharma2017b,Seiler2019} respectively. The \citet{Finkelstein2019} $\dot{N}_\textrm{ion}$ model slightly increases $z\sim4\rightarrow10$, due to low mass halos being more prevalent at higher redshifts with a steep LF faint-end slope ($\alpha < -2$). As demonstrated in panel (b), with the $\tau_\mathrm{es}$ and dark fraction constraints, 21\% of our $\dot{\mathbf{N}}_\textbf{ion}$ posterior increases over $z=6\rightarrow8$ ($d\log_{10}\dot{\mathbf{N}}_\textbf{ion}/dz|_{z=6\rightarrow8} > 0$), and only 3\% when the Ly$\alpha$\ damping constraints are included, making an increasing $\dot{N}_\textrm{ion}(z)$ much more unlikely. The \citet{Naidu2019} SFR surface-density $f_\textrm{esc}$ model shifts the burden of reionization onto brighter galaxies, producing a steeper rise in $\dot{N}_\textrm{ion}$ with decreasing redshift \citep[see also,][]{Sharma2018,Seiler2019}. We note that the \citet{Naidu2019} model was fit to most of the same constraints as presented here, so is likely to match by construction, but demonstrates a physically-motivated model which can match the inferred evolution in $\dot{N}_\textrm{ion}$.

The right axes of panels (c) and (d) show the number of ionizing photons per hydrogen atom released in a Hubble time: $\dot{N}_\textrm{ion}/\langle n_\mathrm{H} \rangle/H(z)$, calculated at $z=6$. We infer $\sim5-23$ (68\% range) ionizing photons per hydrogen atom from our $\dot{\mathbf{N}}_\textbf{ion}$ posteriors including the Ly$\alpha$\ damping constraints ($\sim3-36$ with only $\tau_\mathrm{es}$ and the dark fraction), supporting a scenario where ionizing photons are relatively abundant in the late stages of reionization \citep{Becker2013}.

\subsection{Reionization history and CMB optical depth}
\label{sec:results_history}

\begin{figure}
	\includegraphics[width=\columnwidth]{figure2.pdf}
    \vspace{-16pt}
    \caption{The volume-averaged IGM neutral hydrogen fraction as a function of redshift, ${\overline{x}_\textrm{{\textsc{hi}}}}(z)$ from our inferred $\dot{N}_\textrm{ion}(z)$. Grey lines show 1000 draws from the posterior using only the CMB optical depth, $\tau_\mathrm{es}$ and dark fraction constraints, dark purple shows the inferred evolution using the additional observations, which prefers a later reionization. We also plot the constraints used in our likelihood (Section~\ref{sec:nion_like}).}
    \label{fig:history}
\end{figure}

\begin{figure}
	\includegraphics[width=\columnwidth]{figure3.pdf}
	\vspace{-16pt}
    \caption{1D and 2D distributions of reionization mid-point $z({\overline{x}_\textrm{{\textsc{hi}}}} = 0.5)$ and duration $\Delta z_\mathrm{reion}$ from our inferred $\dot{N}_\textrm{ion}(z)$. The grey contours shows the result using only the CMB optical depth, $\tau_\mathrm{es}$ and dark fraction constraints, dark purple shaded contours show the constraints using the additional Ly$\alpha$\ damping observations, which prefer a later reionization. Contours in the 2D plot show 68\% and 95\% confidence regions.}
    \label{fig:zreion}
\end{figure}

\begin{figure}
	\includegraphics[width=\columnwidth]{figure4.pdf}
    \vspace{-16pt}
    \caption{Electron scattering optical depth to CMB photons, $\tau_\mathrm{es}$. The \citet{PlanckCollaboration2015} and \citetalias{PlanckCollaboration2018} measurements are plotted in grey (light dotted, dark dashed respectively) assuming Gaussian errors. We plot the distribution of $\tau_\mathrm{es}(z=15)$ from our inferred $\dot{\mathbf{N}}_\textbf{ion}$. As required to maximise our likelihood, $\tau_\mathrm{es}$ inferred from our model $\dot{\mathbf{N}}_\textbf{ion}$ is fully consistent with the \citetalias{PlanckCollaboration2018} measurement (light grey distribution -- using the \citetalias{PlanckCollaboration2018} $\tau_\mathrm{es}$ and dark fraction, purple distribution -- including Ly$\alpha$\ damping constraints).}
    \label{fig:tau}
\end{figure}

Figure~\ref{fig:history} shows the reionization history obtained in our inference and the observational constraints used in the likelihoods (Section~\ref{sec:nion_like}). As required by the likelihoods our model fits the observations by construction, but it is instructive to observe the allowed reionization histories. While only using the most model-independent constraints allows a broad range of reionization histories \citep[see also][]{Greig2017c}, the Ly$\alpha$\ damping constraints prefer a later reionization.

The reionization timeline we infer is rapid between $z\sim6-8$, driven by the Ly$\alpha$\ damping measurements, and thus does not overlap significantly with an extended timeline, such as the model by \citet{Finkelstein2019}. This is primarily due to our sharply decreasing $\dot{N}_\textrm{ion}(z)$. In the \citet{Finkelstein2019} model low mass galaxies with extremely high escape fractions produce flatter $\dot{N}_\textrm{ion}(z)$ at early times than we infer (see panel (d) of Figure~\ref{fig:nion}), enabling reionization to start early, while more massive galaxies with lower escape fractions dominate the SFR density at later times, as the galaxy LF faint-end slope flattens, extending reionization. Our results imply that the latest reionization constraints favour an alternative balancing of the ionizing photon budget, likely reducing the dependence on very low mass galaxies with high $f_\textrm{esc}$ to complete reionization.

Double reionization (where the neutral fraction dips again at high redshift) is not ruled out \textit{a priori} in our model and can be seen in a few posterior draws. In most physical models the luminosity density of galaxies steadily declines at $z>10$ \citep[e.g.,][]{Mason2015,Tacchella2018a}, making double reionization unlikely. Within our framework, without placing model-dependent priors on the evolution of $\dot{N}_\textrm{ion}$, ruling out this scenario requires higher redshift measurements of the reionization timeline.

Figure~\ref{fig:zreion} shows the distributions of reionization's midpoint and duration obtained from our analysis. We define the midpoint: $z_\textrm{0.5} = z({\overline{x}_\textrm{{\textsc{hi}}}}=0.5)$, and the duration: $\Delta z_\mathrm{reion} = z({\overline{x}_\textrm{{\textsc{hi}}}} = 0.9) - z({\overline{x}_\textrm{{\textsc{hi}}}} = 0.1)$. With just the \citetalias{PlanckCollaboration2018} and dark fraction constraints, we obtain $z_\textrm{0.5} = 7.49_{-0.57}^{+0.71} \, (68\%)$ and $\Delta z_\mathrm{reion} = 2.26_{-0.99}^{+2.00}$, comparable to the \citetalias{PlanckCollaboration2018} analysis. Including the Ly$\alpha$\ damping constraints more tightly constrains the reionization midpoint to later times but does not significantly update its duration: $z_\textrm{0.5} = 6.93\pm0.14$ and $\Delta z_\mathrm{reion} = 2.83_{-0.99}^{+2.40}$. Again, this is due to the current lack of constraints on reionization's early stages ($z>8$).

Figure~\ref{fig:tau} shows the electron scattering optical depth to the CMB obtained in our inference. As required in our likelihood, both models are consistent with the \citetalias{PlanckCollaboration2018} $\tau_\mathrm{es}$ at $z=14$. Using the full ${\overline{x}_\textrm{{\textsc{hi}}}}$ constraints we infer $\tau_\mathrm{es}=0.053\pm0.004$.

\section{Discussion}
\label{sec:disc}

The CMB optical depth and dark fraction provide the most model-independent constraint on $\dot{N}_\textrm{ion}(z)$. Any model of reionization should produce $\dot{N}_\textrm{ion}(z)$ and reionization histories which are consistent with these results. As noted in Section~\ref{sec:results_nion}, the $\dot{N}_\textrm{ion}(z)$ we infer from these constraints is consistent with the decline in the galaxy UV luminosity density, with fixed $f_\textrm{esc}$ and $\xi_\textrm{ion}$, as shown in previous work \citep[e.g.,][]{Bouwens2015,Robertson2015}. As measurements of the reionization history improve, with further Ly$\alpha$\ observations and 21cm experiments, our knowledge of $\dot{N}_\textrm{ion}$ will increase and can be more directly decomposed into its constituent sources.

Figure~\ref{fig:nion} shows the Ly$\alpha$\ absorption constraints prefer a more rapid decline of $\dot{N}_\textrm{ion}$ over $z\sim6\rightarrow8$ to produce later reionization. While a contribution from quasars could be invoked to produce the build-up of $\dot{N}_\textrm{ion}$ from $z\sim8\rightarrow6$, using the model described in Section~\ref{sec:reion} the contribution of quasars to $\dot{N}_\textrm{ion}$ is negligible at $z>5$: the combined galaxies and quasar models shown in Figure~\ref{fig:nion} do not produce an evolution as steep as our inferred $\dot{N}_\textrm{ion}$. This is true even assuming quasars exist to $M_\textsc{uv} < -18$ and have hard ionizing spectra, $\alpha=-0.56$ \citep{Scott2004}. However, the faint-end slope of the quasar LF is still uncertain at $z>7$ \citep[e.g.,][]{Manti2016}: a large population of faint quasars at $z\sim8-10$ could contribute to the increase in $\dot{N}_\textrm{ion}$ at $z<8$.

The inferred rapid decline in $\dot{N}_\textrm{ion}$ is steeper than the evolution in the galaxy luminosity density (if integrated down to $M_\textsc{uv} \geq -13$). While many different effects contribute to $\dot{N}_\textrm{ion}$, using this luminosity density and constant $f_\textrm{esc}$ and $\xi_\mathrm{ion}$ it thus is difficult to describe $\dot{N}_\textrm{ion}$ evolution. Our inferred $\dot{N}_\textrm{ion}$ therefore provides tentative \textit{a posteriori} evidence that galaxies' ionizing photon emission properties evolve with redshift and/or that frequently-solicited, but undetected, faint galaxies contribute less to reionization than previously required \cite[e.g.,][]{Bouwens2015,Robertson2015}.
For example, models where $f_\textrm{esc}$ is higher in more massive galaxies will naturally produce a steeper $\dot{N}_\textrm{ion}$ evolution, due to the more rapid build-up of high mass galaxies at high redshifts \citep{Sharma2018,Seiler2019,Naidu2019}. A steep $\dot{N}_\textrm{ion}$ can also be produced and/or enhanced if the galaxy luminosity density is dominated by `bright' galaxies ($M_\textsc{uv} \simlt -18$), which is the case if the UV LF faint-end slope is shallow, $\alpha > -2$.
Accurate measurements of the faint-end slope at $z>4$ are therefore important for understanding the ionizing photon budget, and will be improved in the next decade with \textit{JWST}~deep imaging observations.
The dominant sources of reionization impact the topology of HII regions during reionization, measurable by future 21cm intensity mapping experiments: e.g. bright, massive sources produce more biased and larger ionized bubbles compared to if extremely faint, low mass sources dominate the ionizing budget \citep[e.g.,][]{McQuinn2007,Seiler2019}.

As demonstrated in Figure~\ref{fig:history} the Ly$\alpha$\ emitter clustering and EW evolution measurements most tightly constrain the $\dot{\mathbf{N}}_\textbf{ion}$ posterior at $z\sim7$, due to their smaller uncertainties compared to the ${\overline{x}_\textrm{{\textsc{hi}}}}$ measurements from individual quasars. As discussed in Section~\ref{sec:nion_like} the Ly$\alpha$\ damping measurements may introduce systematics due to modelling the intrinsic Ly$\alpha$\ emission (both in galaxies and quasars). In particular, an increase in circumgalactic medium absorption at $z>6$ could also play a role in absorbing Ly$\alpha$\ emission from galaxies \citep[e.g., from an increase in self-shielding systems][]{Bolton2013a}, but is likely subdominant to reionization as explaining the Ly$\alpha$\ damping without reionization requires the ionizing background to drop by at least a factor of 20 from $z\sim6\rightarrow7$ \citep{Mesinger2015}.

We find constraints from Ly$\alpha$\ damping provide increasing evidence that the mid-point of reionization was relatively late, $z\sim7$. However, the duration of reionization is still mostly unconstrained due to the lack of observations at $z>8$. Physical models of the $z>10$ galaxy population \citep[e.g.,][]{Mason2015,Tacchella2018a} predict a declining luminosity density, due to the lower abundance of halos to host star formation, which \textit{should} provide continually declining $\dot{N}_\textrm{ion}(z)$, but our work demonstrates that this is not yet confirmed observationally. $\dot{N}_\textrm{ion}$ is best constrained $z\sim6-8$, but future $z>8$ Ly$\alpha$\ surveys and 21cm experiments will measure ${\overline{x}_\textrm{{\textsc{hi}}}}$, and thus $\dot{N}_\textrm{ion}$, at higher redshifts.

\section{Conclusions}
\label{sec:conc}

We have fit for the non-parametric evolution of the hydrogen-ionizing emissivity, $\dot{N}_\textrm{ion}$, at $z=4-14$, being maximally agnostic about the ionizing sources. We use the most recent constraints on hydrogen reionization to constrain our model. This method enables inference of $\dot{N}_\textrm{ion}(z)$ at redshifts higher than is possible with the Ly$\alpha$\ forest, which becomes too heavily absorbed at $z>6$. Our main conclusions can be summarised as follows:

\begin{enumerate}
    \item Current constraints on reionization favour a declining $\dot{N}_\textrm{ion}$ with redshift at $z>6$. This is moderately favoured by model-independent $\tau_\mathrm{es}$ and dark fraction constraints, $d\log_{10}\dot{\mathbf{N}}_\textbf{ion}/dz|_{z=6\rightarrow8} = -0.31\pm0.35$ and more strongly favoured when Ly$\alpha$\ damping constraints are included ($d\log_{10}\dot{\mathbf{N}}_\textbf{ion}/dz|_{z=6\rightarrow8} = -0.44\pm0.22$).
    %
    \item $\dot{N}_\textrm{ion}(z)$ inferred from model-independent reionization constraints is consistent with the declining UV luminosity density of galaxies, with constant $f_\textrm{esc}$ and $\log_{10}\xi_\mathrm{ion}$, as found by previous studies.
    %
    \item When reionization constraints from Ly$\alpha$\ damping in quasars and galaxies are included, a more rapid decline in $\dot{N}_\textrm{ion}(z)$ is inferred, which is less likely to be explained purely by the declining galaxy UV luminosity density (integrated to $M_\textsc{uv} > -13$) with constant ionizing photon output of galaxies, relative to their non-ionizing UV emission, $z=8\rightarrow6$.
    %
    \item Including the Ly$\alpha$\ damping measurements significantly constrains the midpoint of reionization to $z_\textrm{0.5} = 6.93\pm0.14$ (compared with $7.49_{-0.57}^{+0.71}$ from the \citetalias{PlanckCollaboration2018} optical depth and dark fraction alone), but does not provide more information on its duration due to the lack of observational constraints in reionization's earliest stages, $z>8$.
\end{enumerate}

Without making assumptions about galaxy evolution, our analysis demonstrates current reionization measurements are broadly consistent with the evolution in the galaxy luminosity density, but hint at evolution in the ionizing photon output at $z>6$ and/or a lessened requirement for undetected extremely faint galaxies to dominate reionization. Future measurements of the reionization history, from galaxy surveys and 21cm experiments, can be used in this framework to more tightly constrain the evolution of $\dot{N}_\textrm{ion}$, and combined with statistics of the $z>8$ galaxy population, observable with \textit{JWST}, and measurements of the reionization topology, expose the sources of reionization.

\section*{Acknowledgements}

The authors thank Charlie Conroy for useful discussions and for providing comments on a draft of this paper. We thank Steve Finkelstein, Fred Davies and Brad Greig for useful discussions and for providing: the $\dot{N}_\textrm{ion}$ model from \citet{Finkelstein2019}, and ${\overline{x}_\textrm{{\textsc{hi}}}}$ posteriors from \citet{Davies2018} and \citet{Greig2019} respectively. We thank Austin Hoag for sharing the ${\overline{x}_\textrm{{\textsc{hi}}}}$ posterior from \citet{Hoag2019a}.
CAM acknowledges support by NASA Headquarters through the NASA Hubble Fellowship grant HST-HF2-51413.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. RPN gratefully acknowledges an Ashford Fellowship and Peirce Fellowship granted by Harvard University. ST is supported by the Smithsonian Astrophysical Observatory through the CfA Fellowship. JL is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1701487.

\noindent
\textit{Software}: \verb|IPython| \citep{Perez2007a}, \verb|matplotlib| \citep{Hunter2007a}, \verb|NumPy| \citep{VanderWalt2011a}, \verb|SciPy| \citep{Oliphant2007a}, \verb|Astropy| \citep{Robitaille2013}, \verb|dynesty| \citep{Speagle2019}.

\bibliographystyle{mnras}