\section{Introduction}
\label{sec:Introduction}

The \emph{initial (nucleation) stage} of a first-order phase transition, during which the precursors of the new phase appear as a result of random density fluctuations, still remains a challenging problem in the field of condensed matter. Throughout the last century, numerous theoretical models have been developed with the aim of describing such phenomenon. These can be grouped in two categories depending on how the energy needed to form a molecular aggregate is derived: \emph{i)} the \emph{cluster approach}, originally introduced by the pioneers of the nucleation theory \cite{book:gibbs,article:farkas,article:kaischew-stranski,article:becker-doring}; and \emph{ii)} the \emph{density-functional approach} (DFA), which was applied to nucleation first by Cahn and Hilliard.\cite{article:cahn-hilliard}

In the cluster approach, phase transition is modelled  by the formation of unstable molecular aggregates (clusters) whose density is close to that of the new stable phase. This argument is the cornerstone of both the equilibrium and the kinetic view of nucleation; both of them consider that the aggregation of molecules requires energy, which is known as work for cluster formation $W(N)$, and such magnitude can be calculated in terms of the Gibbs Free Energy (e.g. chapter 3 of \citeauthor{book:kashchiev-1}\cite{book:kashchiev-1}).
The mathematical expression for the work of cluster formation is very simple in the case of spherical clusters since the surface tension can be approximated by the one known between the phases involved in the process (e.g. vapor-liquid, liquid-solid) and does not depend on cluster size. Such approximation is also called \emph{capillary approximation}. In the classical nucleation theory (CNT) the work of cluster formation plays a prominent role because the equilibrium cluster size distribution is of Boltzmann type, i.e. it is exponentially proportional to the work of formation, and satisfies the law of mass action. One of the main problems of this theory is that the pre-exponential factor depends on the concentration of potential nucleation sites, an obscure parameter which cannot be determined precisely. Although CNT is easily comprehensible, in many cases it does not accurately predict the nucleation rate.\cite{article:viisanen-strey-reiss}  In the case of kinetic nucleation theory (
KNT) the pre-exponential
factor is derived in terms of the monomer attachment and detachment rates.\cite{book:nielsen} Although the expression for the collision frequency is readily obtainable, the same is not true with the monomer desorption rate. The work of cluster formation enters in KNT via the detailed-balance condition imposed to eliminate the dependency on this unknown quantity. Considering this assumption one obtains a Fokker-Planck equation (FPE) as the time-evolution equation of the cluster size concentration, better known as the Zeldovich-Frenkel equation (e.g. chapter 9 of \citeauthor{book:kashchiev-1}\cite{book:kashchiev-1}). Hence, the equilibrium cluster size distribution can be re-written in terms of the work of formation and the monomer attachment rate yielding the kinetic counterpart of the classical equation. In this framework, the estimated values for nucleation rate are  also far from those obtained experimentally. These deviations led to the development of new theories such
as the Density Functional Approach.

Under DFA, the state of the system is described by the number density of molecules as a function of the space coordinates, $\rho(\mathbf{r})$. No explicit boundary between the phases of higher and lower density is assumed and, hence, it cannot accurately be established whether a molecule belongs to the old or the new phase. This is also known as the diffuse interface theory of nucleation. In this regard, the work of cluster formation is now expressed in terms of the Landau potential which can be calculated by using, for example, the square-gradient approximation\cite{article:lutsko-2011b} and used as a replacement of $W(N)$ in CNT and KNT master equations.

Until recently, the above descriptions of  nucleation  have been the main tools to predict experimental results. Notwithstanding, the underlying mathematical treatment in these theories is completely deterministic, though they describe a random process. For this reason a new formulation for the nucleation theory has been developed within the statistical mechanics framework.\cite{article:ford,article:ford-review} However, this formalism cannot be used in the case of non-stationary conditions and, furthermore, it does not fit to the observed values of (nucleation) induction time.\cite{article:ford} The main motivation of this paper is to overcome both limitations, i.e. the deterministic nature of previous models and the lack of solutions for non-stationary conditions.

\subsection{Outline of the model}

In the middle of the twentieth century, Hendrik Anthony Kramers described thermally activated escape from a metastable state as \emph{Brownian motion} of a fictitious particle along a \emph{reaction coordinate}, which covers all relevant degrees of freedom that define the system state, in a (static) field of force.\cite{article:kramers} Such processes are characterized by the presence of two stable states separated by a potential barrier and can be described in two different but equivalent forms. The first is using stochastic differential equations (SDE) to characterize the time evolution of the reaction coordinate. The conjectured SDE is a Langevin-like equation\cite{article:langevin} which contains both a deterministic and a random part. The second is using the time-evolution equation of the probability density function (PDF) associated with the reaction coordinate. Indeed, this expression can be written as a Fokker-Planck equation.\cite{book:risken-1996a} Thus, as we will demonstrate later,
using the FPE one can derive an analytical expression for the nucleation (or escape) rate in the form proposed by Arrhenius. Moreover,
considering the nucleation process as a rare event one can obtain the PDF for the induction time.  Hence, the proposed model not only allows to reduce the problem from a $N$-dimensional to an unidimensional description but also provides analytical expressions for magnitudes such as nucleation rate under both stationary and non-stationary conditions, which demonstrates the abilities of the model.

The aim of this work is to apply Kramers' reasoning to the problem of nucleation in order to predict the induction time statistics for both stationary and non-stationary conditions, the latter assuming the semiadiabatic limit.\cite{article:talkner-luczka,article:kim-talkner-kyun-Hanggi} We propose cluster radius, $R$, as a reaction coordinate for crystal nucleation and, hence, a stochastic differential equation is postulated for modeling time evolution of the cluster size. In this context, analytical expressions for statistics of induction time can be deduced by employing the theory of stochastic processes. This contribution allows us to better understand the randomness of experimental results and the deviation from the theoretical (classical) nucleation rate.

\section{Equation of motion for the cluster radius}

\subsection{Theoretical background}

The energy required to form (so-called work of formation of) a spherical cluster with the properties of the new phase within a supersaturated mother phase has been extensively studied.\cite{article:lutsko-2011b,article:laaksonen,article:debenedetti,article:nishioka} Although historically the usual convention has been to evaluate the work of cluster formation in terms of the increment of the Gibbs free energy, such a magnitude has to be evaluated depending on the system in terms of the corresponding thermodynamic potential. Thus, if one considers homogeneous nucleation in the context of the Grand Canonical ensemble, the work of formation will be given in terms of the Landau potential, $\Omega$. Nevertheless one can consider the problem in a closed system and, then, one has to evaluate this energy in terms of the Helmholtz free energy, $F$. In any case, Gibbs' interfacial thermodynamics can be used in order to
derive the energetic cost of forming an spherical embryo of  the stable phase inside
the metastable mother phase, taking into account the capillary approximation. These calculations end up in expressions that contains a volume and a surface term.
The former consists of a free energy density term multiplied by the volume of the sphere. The latter contains the surface of the sphere multiplied by a surface tension term which generally depends on the radius. Nevertheless, as has been shown, this dependence is weak and considering the surface free energy as a constant is a very good approximation. However, computing the free energy density inside the cluster is not a straightforward calculation neither in monocomponent nor in multicomponent systems.\cite{article:laaksonen,article:debenedetti,article:nishioka} For the sake of simplicity, we will consider the approximation of such magnitude as the increment of the chemical potential\cite{book:zettlemoyer}, obtaining so the following expression for the reversible work required to form a noncritical cluster\cite{article:nishioka-1992}
\begin{equation}
 W(R)=\Delta \Omega(R)=-\frac{4\pi\,\Delta\mu}{3v_0}R^3+4\pi\sigma R^2,
 \label{eq:1}
\end{equation}
where $\Delta \mu$ is the difference between the chemical potential of the liquid and that of the solid phase, $v_0$ is the volume occupied by a molecule in the cluster of the new phase and $\sigma$ is the surface tension between the old and the new phase (considering a planar interface). At small $R$, the second term of equation (\ref{eq:1}), which is related to the formation of the surface, prevails while the first term plays a decisive role at large values of $R$. The maximum of $W(R)$ defines the critical size,
\begin{equation}
 R^\star=\frac{2v_0\sigma}{\Delta\mu}\label{eq:2},
\end{equation}
with
\begin{equation}
 W^\star=\frac{16\pi}{3}\frac{v_0^2\sigma^3}{\Delta\mu^2}\label{eq:3}.
\end{equation}

As widely known, equation (\ref{eq:1}) characterizes the equilibrium state of the system whereas the cluster radius describes the dynamic of the nucleation process. The latter magnitude is generally strongly coupled with the environmental degrees of freedom and due to such a coupling its dynamic is not deterministic but stochastic. Therefore, one can assume that $R$ plays the role of a reaction coordinate in the Kramers' theory.\cite{book:berne,article:hanggi-talkner-borkovec} Accordingly the variable cluster radius, hereafter denoted as $X_R$, can be understood as a stochastic process, i.e. $X_R(t)$ draws a random trajectory in the reduced phase space\footnote{The mathematics of such a reduction were developed in the framework of statistical mechanics\cite{book:kubo-2}} of the system (figure \ref{fig:2}). Considering these assumptions, the nucleation process is characterized by a stochastic variable which can be interpreted as the instantaneous position of a  \emph{fictitious Brownian particle}. Inspired
by the Kramers' theory one postulates that the equation of motion for such a particle will be given by a Langevin-type equation. In fact, nucleation is nothing but a thermally activated escape process from a potential barrier. Therefore, in order to escape from the well $A=\{X_R\leq X_{R}^\star\}$, the fictitious ``random walker'' must acquire energy to overcome the energy barrier and subsequently it must again lose energy to become trapped by the attractor $B=\{X_R\geq X_R^\star\}$, i.e. a molecular cluster will fluctuate in size until escape from the well $A$ to the attractor $B$ after which the cluster will grow in a deterministic manner.

\subsection{The model}

The stochastic dynamic of the cluster radius, $X_R$, can be phenomenologically postulated by means of an \emph{overdamped} \emph{Langevin} \emph{equation},
\begin{align}
  \textcolor{black}{\eta_{\mathcal{E}}\,\frac{d{X}_R(t)}{dt}=-\frac{\partial W(X_R(t))}{\partial X_R}+\sqrt{2\eta_{\mathcal{E}} k_BT}\xi(t)}
 \label{eq:4}
\end{align}
where $\eta_{\mathcal{E}}$ is the friction coefficient associated with the reduced phase space (so-called \emph{reaction-coordinate viscosity}) and $\xi(t)$ denotes zero-mean, delta-correlated Gaussian white noise (GWN),
  \begin{align}
   \langle\xi(t)\rangle&=0,\nonumber\\
   \langle\xi(t)\xi(t')\rangle&=\delta(t-t')\label{eq:5}.
  \end{align}
The main reason to consider the overdamped limit (avoiding the second time derivative of $X_R$) is because it is consistent with the structure of the classical post-critical growth law. To demonstrate that we only have to consider $X_R>X_R^\star$, where the noise term is negligible in comparison to the driving force, so
$$
\frac{dX_R}{dt} \sim -\eta_{\mathcal{E}}^{-1}\frac{\partial W(X_R)}{\partial X_R}
$$
which shows the same structure as (see Eq. 2.62 of \citeauthor{book:kelton}\cite{book:kelton})
$$
\frac{d N}{dt}=-\frac{f(N)}{k_BT}\frac{\partial W(N)}{\partial N},
$$
with $N$ the number of molecules inside the cluster and $f(N)$ the attachment rate. Although inspired by the latter one may feel tented to consider a size-dependent viscosity, $\eta_{\mathcal{E}}$, we used a constant viscosity since it is in good agreement with the procedure followed in CNT to calculate nucleation rates setting $f(N)$ to $f(N^\star)$ (e.g. p.168 of \citeauthor{book:kashchiev-1}\cite{book:kashchiev-1}), with the advantage that
this simplifies the later calculations. In order to confirm that such a simplification is reasonably good, we can consider for example the case of diffusion-limited kinetics, where the monomer attachment rate is given by (see Eq. 10.18 of \citeauthor{book:kashchiev-1}\cite{book:kashchiev-1})
$$
f(N)=\gamma_N\,4\pi\,D\,\rho_\infty\,X_R(N)
$$
where $\gamma_N\simeq1$ is the sticking coefficient, $D$ is the diffusion constant and $\rho_\infty$ is the monomer number density and where the only dependence on the size enters via the radius $X_R(N)$ with $N=\frac{4\pi}{3v_0}X_R^3$. As can be observed in figure \ref{fig:R_N}, which represents  $X_R(N)$ as a function of $N$, the radius lies in the range $\sim[2,6]$ for a very wide range of $N$. Thus, one can conclude that $\eta_{\mathcal{E}}$ could be ultimately approximated as a constant to make easier the later mathematical treatment. In fact, as we will see later (Eq. \ref{eq:21-new}), if one selects $\eta^{-1}_{\mathcal{E}}=(\frac{\partial N}{\partial X_R})_{X_R^\star}^{-2}f(N^\star)/k_BT$ inspired by the previous reasoning, the nucleation rate derived from CNT is recovered except for a multiplicative factor.

It is worth to note that the friction coefficient $\eta_{\mathcal{E}}$ is an abstraction that provides information about the viscosity of the phase space, $\mathcal{E}$. Hence, there is no trivial expression for $\eta_\mathcal{E}$ in terms of the friction acting on individual molecules in the real space. Nonetheless this magnitude can be estimated, together with the mass associated with the reaction coordinate $m_X$, by using molecular dynamics.\cite{article:huang-attard} One can expect that the former will be related to the mean monomer attachment frequency for a $R$-sized spherical cluster $\overline{f(X_N)}$ since, the first term of the right hand side of equation (\ref{eq:4}) informs about the deterministic behaviour of the cluster size and the second one only can arise from the unpredictable collisions of monomers with the cluster. Indeed, $\eta_{\mathcal{E}}$ should be characterized by the thermodynamic properties of the new phase since at a fixed temperature the monomer attachment frequency depends on
the density of the final state.

Thus the most remarkable characteristic of equation (\ref{eq:4}), comparing with the Ginzburg-Landau equation (e.g. chapter 4 of \citeauthor{book:barrat-hansen-2003-a}\cite{book:barrat-hansen-2003-a}), is the second term on the right hand side which includes not only a random variable but also the temperature of the thermal bath. Hence, this term plays the role of a fluctuating force which comprehends all degrees of freedom associated with the environment. Therefore, the equation of motion of the fictitious Brownian particle presents two contributions: \emph{i)} the deterministic force due to the free energy potential $\textcolor{black}{W}$, and \emph{ii)} the random force \hbox{$\widetilde{\xi}(t):=\sqrt{2\eta_{\mathcal{E}} k_BT}\xi(t)$}.

From the stochastic differential equation (\ref{eq:4}), using the forward Kramers-Moyal expansion, it can be demonstrated that the time-evolution equation of the PDF, $\rho(X_R,t)$, is
\begin{align}
 \frac{\partial \rho}{\partial t}(X_R,t)&=\hat{\mathcal{L}}_{_{SL}}\rho(X_R,t)=-\frac{\partial j(X_R,t)}{\partial X_R}\label{eq:6},
\end{align}
where the linear differential operator $\hat{\mathcal{L}}_{_{SL}}$ is the \emph{Smoluchowski operator},
\begin{equation}
 \hat{\mathcal{L}}_{_{SL}}:=\frac{\partial}{\partial X_R}\left(\frac{1}{\eta_{\mathcal{E}}}\frac{\partial\textcolor{black}{W}(X_R)}{\partial X_R}+\frac{k_BT}{\eta_{\mathcal{E}}}\frac{\partial}{\partial X_R}\right)\label{eq:7},
\end{equation}
a special form of the \emph{Fokker-Planck} operator $\hat{\mathcal{L}}_{FP}$ \citep{book:risken-1996a}, and $j$ is the probability current
\begin{equation}
  j(X_R,t):=-\left(\frac{1}{\eta_{\mathcal{E}}}\frac{\partial\textcolor{black}{W}(X_R)}{\partial X_R}+\frac{k_BT}{\eta_{\mathcal{E}}}\frac{\partial}{\partial X_R}\right)\rho(X_R,t)\label{eq:8}.
\end{equation}
Equation (\ref{eq:6}) is also known as the \emph{Smoluchowski equation} and belongs to the family of Fokker-Planck partial differential equations. It is noteworthy that the above FPE for the proposed continuous random walker should be considered as the continuous counterpart of that derived from a discrete Brownian motion\cite{article:white-1969,book:kashchiev-1} but setting the monomer attachment rate to be the mean monomer attachment frequency, which is directly related with $\eta_{\mathcal{E}}$ by means of equation (\ref{eq:20-new}). The main advantage of having a continuous version of the PDF is that it allows us to use the tools of the continuous calculus to compute in an easy way magnitudes such as the mean first passage time without considering infinite series. Moreover, as will be immediately studied, obtaining the stationary and quasi-stationary distribution functions will be almost straightforward owing to the rules of continuous integral calculus can
be used.

\subsection{Stationary probability density function}

The time evolution of the PDF converges to the stationary solution of equation (\ref{eq:6}) when \hbox{$t\rightarrow\infty$}. With the aid of this solution it is possible to estimate the probability to find a cluster of a given size when the time is much larger than the relaxation time $\tau_S$. This represents the time required to decay within the attractor $A$
\begin{equation}
 \tau_s\sim\left[\frac{1}{m_X}\left(\frac{\partial^2 \textcolor{black}{W} }{\partial^2 X_R}\right)_{X_R=0}\right]^{-1/2},
 \label{eq:9}
\end{equation}
with $m_X$ the effective mass corresponding to the fictitious particle. In such a case, i.e. when the time is larger than the time required for the spontaneous decomposition of a subcritical cluster, the non-stationary PDF converges to the stationary one
\begin{align}
\rho_{\text{st}}(X_R)&=\frac{1}{\zeta}
  \exp\left\{-\int^{X_R}\frac{1}{k_BT}\left(\frac{\partial\textcolor{black}{W} (X_R')}{\partial X_R'}\right)dX_R'\right\}\nonumber\\
  &=\frac{1}{\zeta}
  \exp\left\{-\frac{\textcolor{black}{W} (X_R)}{k_BT}\right\},\label{eq:10}
  \end{align}
where $\zeta$ is the normalization constant (see \ref{ap:1})
\begin{align}
\zeta&=\int_{0}^{\infty}\rho_{st}(y)dy\simeq
\sqrt{\frac{\pi k_BT}{2}}\frac{\exp\left\{
-\frac{\textcolor{black}{W}(X_R=0)}{k_BT}\right\}}{\sqrt{\left(\frac{\partial^2\textcolor{black}{W} }{\partial X_R^2}\right)_{X_R=0}}}\textcolor{black}{=\sqrt{\frac{\pi k_BT}{8\pi\sigma}}}
\label{eq:11}
\end{align}
The identity (\ref{eq:10}) is completely in accordance with the cluster size distribution of CNT\cite{article:izmailov-myerson-arnold} with $\zeta$ playing the role of the pre-exponential factor.

Nevertheless, the stationary PDF is not useful to obtain the nucleation rate since the boundary condition $X_R\geq0$ implies that the probability current must be zero at $X_R=0$, i.e. $j(0)=j(X_R)=0$, and hence at any value of $X_R$ including $X_R^\star$. Such restriction implies that the escape rate must be equal to zero due to the definition of this magnitude,
\begin{equation}
 k^+:=\frac{j}{n},\label{eq:12}
\end{equation}
with $n$ the stationary probability that the particle has not crossed the boundary $X_R^->X_R^\star$,
\begin{equation}
 n=\int_0^{X_R^-}\rho_{st}(X_R)dX_R,
  \label{eq:13}
\end{equation}
which is usually approximated by $n\sim1$. This fact shows that the stationary distribution is a good but unrealistic approximation. For this reason we shall introduce the quasi-stationary solution of the Smoluchowski equation.

\subsection{Quasi-stationary probability density function}

Now we shall assume that $\textcolor{black}{W/k_BT\,>\,1}$ and $T$ is constant. Under these conditions, before reaching the stationary state  (i.e. $\tau_s\ll t<\infty$)  the Smoluchowski solution remains in a quasi-stationary state and, therefore, the current probability is almost time-independent. In such a quasi-stationary state, the probability current over the top of the potential is very small near to $X_R^\star$ and the time change of the PDF is also very small. Therefore, the small value of the probability current is almost independent of $X_R$, i.e. $j(X_R,t)\sim j$. Hence, equation (\ref{eq:8}) can be written as
\begin{align}
j&\simeq-\frac{k_BT}{\eta_{\mathcal{E}}}e^{-\textcolor{black}{W} (X_R)/k_BT}\frac{\partial}{\partial X_R}\left(e^{\textcolor{black}{W} (X_R)/k_BT}\rho(X_R,t)\right)
\label{eq:14}.
\end{align}
Considering now the boundary condition $\rho(X_R^-,t)=0$, one readily gets
\begin{align}
 \rho_{\text{q-st}}(X_R,t)&=\frac{\eta_{\mathcal{E}} j}{k_BT}e^{-\textcolor{black}{W} (X_R)/k_BT}\int_{X_R}^{X_R^-} dy\,e^{\textcolor{black}{W} (y)/k_BT}
  \label{eq:15}.
\end{align}

\section{Nucleation rates and induction time statistics}

We want to calculate the mean time that the system needs to produce a supercritical cluster, i.e. when $X_R>X_R^{\star}$. This magnitude is also called \emph{mean first-passage time} (MFPT) (or Kramers' time) and can be easily related to the induction time. In this section we present the expressions obtained for the induction time and the nucleation rate under both stationary and non-stationary conditions.

\subsection{Nucleation rate under steady conditions}

From equations (\ref{eq:12},\ref{eq:13},\ref{eq:15}) the following expression for the escape rate can be derived,
\begin{align}
 \frac{1}{k^+}&=\int_{0}^{
X_R^-} dz\,\frac{\eta_{\mathcal{E}}}{k_BT}e^{-W(z)/k_BT}\int_{X_R}^{X_R^-} dy\,e^{W(y)/k_BT}\label{eq:16}.
\end{align}
Accordingly, applying the Laplace (or Gaussian steepest-descent) method (e.g. page 124 of cite \citeauthor{book:risken-1996a}\cite{book:risken-1996a})  the following approximation can be obtained (see \ref{ap:2})
\begin{align}
k^+&=\sqrt{\frac{\partial^2W(0)}{\partial X_R^2}}\sqrt{\left|\frac{\partial^2W(X_R^\star)}{\partial X_R^2}\right|}\frac{e^{-W^\star/k_BT}}{\pi\eta_{\mathcal{E}} }\textcolor{black}{=4\pi\sigma\frac{e^{-W^\star/k_BT}}{\pi\eta_{\mathcal{E}} }}
\label{eq:17}.
\end{align}
The inverse of this escape rate is also known as MFPT,
\begin{equation}
 \tau_K=\frac{\pi\eta_{\mathcal{E}} }{\sqrt{\frac{\partial^2W(0)}{\partial X_R^2}}\sqrt{\left|\frac{\partial ^2W(X_R^\star)}{\partial X_R^2}\right|}}e^{W^\star/k_BT}.
  \label{eq:18}
\end{equation}
The nucleation rate and the induction time are given by the expressions\cite{book:barrat-hansen-2003-a}
\begin{align}
J=\frac{\rho_\infty}{\tau_K}&= \rho_\infty\,k^+=\rho_\infty\sqrt{\frac{\partial^2W(0)}{\partial X_R^2}}\sqrt{\left|\frac{\partial^2W(X_R^\star)}{\partial X_R^2}\right|}\frac{e^{-W^\star/k_BT}}{\pi\eta_{\mathcal{E}} }\label{eq:19-2},\\
t_{\text{ind}}&=\frac{1}{JV}=\frac{\tau_K}{\mathcal{N}_1}=\frac{\pi \eta^{\text{eff}}}{\sqrt{\frac{\partial^2W(0)}{\partial X_R^2}}\sqrt{\left|\frac{\partial ^2W(X_R^\star)}{\partial X_R^2}\right|}}e^{W^\star/k_BT},\label{eq:indtime}
\end{align}
where $\rho_\infty$ is the monomer number density, i.e. the equilibrium number of monomer per unit volume $\mathcal{N}_1/V$, and $\eta^{\text{eff}}=\eta_{\mathcal{E}}/\mathcal{N}_1$. As can be observed, equation (\ref{eq:19-2}) looks like the classical expression for the nucleation rate. Indeed, it can be rewritten as,
\begin{align}
J&=\rho_\infty\frac{\textcolor{black}{2\,k_BT}\,\left(\frac{\partial X_N}{\partial X_R}(X_R^\star)\right)}{\eta_{\mathcal{E}}}\,Z_0\,Z_D\,e^{-W^\star/k_BT},\label{eq:19-new}
\end{align}
with
\begin{align}
 X_N&=\frac{4\pi}{3 v_0}X_R^3,\nonumber\\
  Z_0&=\left(\frac{1}{\textcolor{black}{2}\pi\textcolor{black}{k_BT}}\frac{\partial^2W(0)}{\partial X_R^2}\right)^{\frac{1}{2}}\textcolor{black}{=\sqrt{\frac{4\pi\sigma}{2\pi k_BT}}},\nonumber\\
 Z_D&=\left(\frac{1}{2\pi k_BT}\left|\frac{\partial^2W(X_N^\star)}{\partial X_N^2}\right|\right)^{\frac{1}{2}},\nonumber
\end{align}
being the number of molecules inside the $R$-sized cluster, the curvature of the energy landscape at the basin and the Zeldovich's factor (typically\cite{book:kashchiev-1}\ $10^{-2}\leq Z_D\leq1$), respectively. Therefore, the term $\textcolor{black}{Z_0}\frac{\textcolor{black}{2k_BT}}{\eta_{\mathcal{E}}}\left(\frac{\partial X_N}{\partial X_R}(X_R^\star)\right)$ sets the time scale of the phase transition. In fact,  following the classical reasoning (that the characteristic time for nucleation is determined by the attachment rate of monomers to the critical cluster) one could postulate that
\begin{equation}
Z_0\frac{\textcolor{black}{2k_BT}}{\eta_{\mathcal{E}}}\left(\frac{\partial X_N}{\partial X_R}(X_R^\star)\right)\sim  f(X_R^\star) \label{eq:22-22}
\end{equation}
so that
\begin{align}
  \eta_{\mathcal{E}}^{-1}&\sim Z_0^{-1}\frac{f(X_R^\star)}{\textcolor{black}{2k_BT}}\left(\frac{\partial X_N}{\partial X_R}(X_R^\star)\right)^{-1}\label{eq:20-new}
\end{align}
Thus, substituting equation (\ref{eq:20-new}) into (\ref{eq:19-new}) one gets,
\begin{equation}
J\sim \rho_\infty\,Z_D\,f(X_N^\star)\,e^{-W^\star/k_BT}\textcolor{black}{=J_{\text{CNT}}}.\label{eq:21-new}
\end{equation}
Therefore, making such an interpretation of the viscosity parameter, one recovers the CNT expression of the nucleation rate\cite{book:kashchiev-1}, $J_{_\text{CNT}}$. Note that this expression of $\eta_\mathcal{E}$ is slightly different from that one we expected, i.e. $\eta^{-1}_{\mathcal{E}}= \left(\frac{\partial X_N}{\partial X_R}\right)_{X_R^\star}^{-2}f(N^\star)/k_BT$, inspired by the post-critical growth rate. In such a case, the nucleation rate would be
$$
J\sim \rho_\infty\,2\,Z_0\,Z_D\,\left(\frac{\partial X_N}{\partial X_R}\right)_{X_R^\star}^{-1}\,f(X_N^\star)\,e^{-W^\star/k_BT}=2\,Z_0\,\left(\frac{\partial X_N}{\partial X_R}\right)_{X_R^\star}^{-1}\,J_{\text{CNT}}
$$
which differs from the CNT expression in a multiplicative pre-exponential factor which is of the order of $Z_D$.

Nonetheless, in this work we propose a similar expression to (\ref{eq:20-new}) but considering the average value of the collision rate that an individual cluster feels, i.e. $\overline{f(X_N)}$. Although this magnitude is completely unknown we can assert that it should be lower than $f(X_N^\star)$. Therefore, our proposed expression for the nucleation rate would be as,
\begin{equation}
 J\sim  \rho_\infty\,Z_D\,\overline{f(X_N)}\,e^{-W^\star/k_BT},\label{eq:22-new}
\end{equation}
which should predict lower values than the classical one since, $\overline{f(X_N)}\leq f(X_N^\star)$. Nonetheless, under this assumption the nucleation rate depends on an unknown magnitude, $\overline{f(X_N)}$. This is the reason by which $\eta_{\mathcal{E}}$ must be fitted to experimental values for $t_{\text{ind}}$. Actually, one could use the fitted value of $\eta_{\mathcal{E}}$ in order to estimate the effective value of\ $\overline{f(X_N)}$\ and, then, better understand the kinetics of the phase transition. With the aid of such a value information on the underlying mechanism of cluster formation could be obtained by comparing with theoretical values obtained using the different expressions for the monomer attachment frequencies (chapter 10 of cite \citeauthor{book:kashchiev-1}\cite{book:kashchiev-1}).

\subsection{Induction time statistics under steady conditions.}

Once an estimation for the induction time has been obtained, it seems interesting to derive an analytical equation for the induction time statistics which would be helpful for a better description of the random nature of such magnitude. Thus, if nucleation is considered as a  homogeneous Poisson process\cite{article:white-1969,book:nucleation-1996,article:peters-2011,article:jiang-2011,article:limay} characterized by the escape rate, $k^+=1/\tau_K$, the induction time can be considered as a Gamma-distributed random variable with a PDF given by
\begin{equation}
 \varrho(t)=\frac{\rho_\infty}{\tau_K}\exp\left\{-\frac{\rho_\infty}{\tau_K}t\right\},
\label{eq:19}
\end{equation}
so that,
\begin{equation}
 P(t_{\text{ind}}\leq t)=\int_0^t\varrho(s)ds=1-\exp\left\{-\frac{\rho_\infty}{\tau_K}t\right\}.
\label{eq:20}
\end{equation}
The latter equation is also called the Kramers law of the escape time statistics. According to equation (\ref{eq:18}), if we know the viscosity of the reduced phase space $\eta_{\mathcal{E}}$ then the statistics of the nucleation process can be estimated from equation (\ref{eq:20}).

Although equations (\ref{eq:18}, \ref{eq:19}, \ref{eq:20}) are good approximations, one can compare the accuracy of these expressions by numerical integration. As in the case of ordinary differential equations, there exists a huge number of techniques (stochastic integrators, SINT) for integrating a SDE as the \emph{ Euler-Maruyama method}.\cite{book:kloeden-1992-a}  In such a case, the upper limit of $X_R$ (i.e. the absorbing wall $X_R^-$) must be fixed at a value greater than $X_R^\star$ and obeying
\begin{equation}
 |\textcolor{black}{W}  (X_R^-)-\textcolor{black}{W} (X_R^\star)|\geq2k_BT.
\end{equation}

\subsection{Time-dependent nucleation rate and induction time statistics under the semiadiabatic limit}

Let us consider the case when the potential barrier $\textcolor{black}{W} $ changes over time (e.g. due to changes in bulk concentration, temperature, pH, etc.\cite{article:fermin-2009}) but this change is slow compared to the relaxation time of the system $\tau_S$. According to the reasoning of  Talkner and co-workers\cite{article:talkner-luczka} the system reaches a quasi-stationary state instantaneously (semiadiabatic approximation) and therefore an analogous deduction to the developed for equation (\ref{eq:20}) can be made for unsteady conditions of $\textcolor{black}{W} $
\begin{equation}
 \tau_K(t)=\frac{\pi\eta_{\mathcal{E}} }{\sqrt{\frac{\partial^2\textcolor{black}{W} (0,t)}{\partial X_R^2}\left|\frac{\partial^2\textcolor{black}{W} (X_R^\star(t),t)}{\partial X_R^2}\right|}}e^{\textcolor{black}{W} ^\star(t)/k_BT},
  \label{eq:22}
\end{equation}
or equivalently,
\begin{equation}
 k^+(t)=\frac{\sqrt{\frac{\partial^2\textcolor{black}{W} (0,t)}{\partial X_R^2}\left|\frac{\partial^2\textcolor{black}{W} (X_R^\star(t),t)}{\partial X_R^2}\right|}}{\pi\eta_{\mathcal{E}} }e^{-\textcolor{black}{W} ^\star(t)/k_BT},
\label{eq:23}
\end{equation}
also called \emph{instantaneous escape rate}. Proceeding in a similar manner to the previous section, one calculates the expression for the \emph{instantaneous nucleation rate} and, hence, for the \emph{instantaneous induction time}
\begin{align}
J(t)&=\rho_\infty(t)\frac{\sqrt{2\,k_BT}\,\left(\frac{\partial X_N}{\partial X_R}(X_R^\star(t))\right)}{\eta_{\mathcal{E}}}\,Z_0(t)\,Z_D(t)\,e^{-W^\star(t)/k_BT},\label{eq:24-new}\\
t_{\text{ind}}(t)&=\frac{1}{J(t)V}=\frac{\tau_K(t)}{\mathcal{N}_1(t)}\label{eq:24-2-new}.
\end{align}

Equation (\ref{eq:22}) constitutes a good approximation only when the topology of the potential barrier $\textcolor{black}{W} (X_R,t)$ does not change, i.e. the character of its local maxima and minima must be invariant. In fact, Talkner and co-workers \cite{article:talkner-luczka} derived a more accurate expression which contains expression (\ref{eq:22}) and a second-order correction term. They called the latter as \emph{geometric correction term} because it is related to the geometric change of the barrier shape. We will not consider this second-order correction in this work.

Following a similar procedure to that of the previous section, the nucleation process can be understood now as a non-homogeneous Poisson process\cite{article:pellerey-franco-shaked-moshe-zinn-joel} and, hence, the first-passage time statistics is governed by the equation below
\begin{equation}
 \varrho(t)=\frac{\rho_\infty(t)}{\tau_K(t)}\exp\left\{-R_K(t)\right\}\label{eq:24},
\end{equation}
with
\begin{equation}
 R_K(t)=\int_0^t \frac{\rho_\infty(s)}{\tau_K(s)}ds\label{eq:25},
\end{equation}
and assuming
\begin{equation}
 R_K(\infty)=\int_0^\infty \frac{\rho_\infty(s)}{\tau_K(s)}ds=\infty\label{eq:26}.
\end{equation}
Accordingly, the distribution function of the induction time is given by
\begin{equation}
 P(t_{\text{ind}}\leq t)=\int_0^t \,\varrho(s)ds=1-e^{-R(t)}\label{eq:27}.
\end{equation}
The accuracy of our theoretical predictions is strongly related to the assumptions  made above. Consequently, the results predicted with the aid of equations (\ref{eq:22}) and (\ref{eq:27}) can be improved using a SINT. In fact, the stochastic integration is indispensable in order to obtain more realistic predictions when the hypotheses of the semiadiabatic limit are not fulfilled. Indeed, SINTs are the only tool to predict escape rates and, then, induction times under strongly unsteady conditions. As mentioned in the previous section, the instantaneous position of the absorbing wall $X_R^-(t)$ must be fixed at a value greater than $X_R^\star(t)$ and obeying the following relation
\begin{equation}
 |\textcolor{black}{W}  (X_R^-,t)-\textcolor{black}{W} (X_R^\star,t)|\geq2k_BT\label{eq:28}.
\end{equation}

\section{A qualitative analysis of the model}

At this point, a short break should be taken in order to summarize and qualitatively analyze the theoretical results obtained so far, before making use of them in the next section.

Thus far the presented model has shown the ability of reproducing the main theoretical results, such as the stationary size distribution, the nucleation rate or induction time equations, of CNT and another previous works in case of setting the viscosity parameter to be the monomer attachment rate of the critical cluster, as we discussed in section 2. Nevertheless, the major difference with previous theories is indeed that the Brownian model does not consider this as the only possibility.
In fact, we have emphasized that there is no apparent growth mechanism to be the same as that which governs the growth of post-critical clusters. Although that could seem a disadvantage because the magnitudes mentioned above will depend on an \emph{a priori} undetermined attachment frequency, it allows us to make better predictions by deducing such a parameter experimentally under conditions when the experiments are highly reproducible and, then, using that fitted value of viscosity, $\eta$, into equations one can predict nucleation rates (or induction times) for different conditions, e.g. large values of supersaturation ratio. Therefore, following this line of reasoning, the gap between the experimental and classical predictions can be overcome, as will be shown in the next section.

Besides, another advantage of this heuristic model is to describe the randomness underlying the experiments in order to determine how reproducible they are. The induction times ineluctably have a standard deviation that cannot be reproduced by the CNT procedures. However, the Brownian model is based on a SDE and hence, the stochasticity is considered. Therefore, the present work endows the classical description with a mathematical apparatus which covers the inhere experimental deviations. That allows us to know whether or not the experimental values lie in the theoretically predicted statistics and, ultimately, knowing whether or not the experiments fulfill
our predictions.

In order to highlight and verify these claims, we will compare the results computed by using the Brownian model against those predicted by using CNT expression (i.e. Eq.(\ref{eq:21-new})). Yet more, the statistics will be tested showing a slight deviation from the experimental one, but this can be explained based on the finite number of assays.

\section{The model at work}

The Brownian model proposed in sections 2 and 3 has been tested {under fixed experimental conditions} by fitting the viscosity parameter through the expression for the induction time (\ref{eq:indtime}) to experimental results of such a magnitude measured for  hydrated calcium sulphate CaSO$_4\cdot$2H$_2$O (gypsum) in a volume of $200\,\mu L$\footnote{I. Rodr\'{i}guez-Ruiz, A.E.S. Van Driessche and J.M. Garc\'{i}a-Ruiz, data to be published.}.
This kind of applications represents one of the main contributions of this paper because neither classical nor non-classical nucleation theories allow to follow a nucleation event in such a simple manner. As a summary, one must fit $\eta_{\mathcal{E}}$ to experimental data using equation (\ref{eq:indtime}) and apply both the SDE (Eq. (\ref{eq:4})) and equations (\ref{eq:19-2}-\ref{eq:20}, \ref{eq:24-new}-\ref{eq:27}) to predict both nucleation rates and induction time statistics under both steady and unsteady conditions.
Once we have the estimation of $\eta_{\mathcal{E}}$, one can apply equation (\ref{eq:20}) in order to calculate the induction time statistics and verify the ability to predict $P(t_{ind}\leq t)$. Moreover, with the aim of considering a finite number of assays the time evolution of the reaction coordinate $X_R$ was simulated using the Euler-Mauryama method. Furthermore, using the fitted value of $\eta_{\mathcal{E}}$ one could simulate more complicated experimental conditions where classical expressions cannot be used.

Twelve assays measuring induction time were considered at a fixed supersaturation and temperature to estimate $\eta_{\mathcal{E}}$, e.g. $\textcolor{black}{\overline{S}}=C/C_e=1.90$\footnote{In solution crystal growth, the supersaturation (ratio between the concetrarion $C$ and the solubility $C_e$) is used as a measure of the driving force for the phase change,
$$
\Delta \mu=k_BT\ln\left(\frac{a}{a_e}\right)\simeq k_BT\ln\left(\frac{C}{C_e}\right),
$$
where $a$ and $a_e$ are the activity and the equilibrium activity of the solute, respectively.} at $T=328K$ (the point marked with an arrow in figure \ref{fig:4}), yielding the value $\eta^{\text{eff}}\sim 2.25$. This value was computed by equaling equation (\ref{eq:indtime}) to the experimental average value at this supersaturation, i.e.
\begin{equation}
\eta^{\text{eff}}=\frac{\frac{1}{M}\sum_{k=1}^{M}t_k^{\,\text{exp}}(\overline{S})}{\frac{1}{4\sigma}\exp\left({W^\star(\overline{S})}/k_BT\right)}
\label{eq:etaeff}
\end{equation}
with $M=12$ in our case.

\noindent{}Substituting the estimated value of $\eta^{\text{eff}}$ into equation (\ref{eq:indtime}) one can calculate a theoretical curve of predicted induction times as a function of supersaturation. The predicted curves for induction times are in very good agreement with the measured values (figure \ref{fig:4}).
To carry out the calculation of these curves, as well as the classical predictions, we need to know the surface tension. In this work
we used\cite{book:volmer,book:kashchiev-1}
\begin{align}
\sigma&=\Psi^{1/3}(\Theta_w)\,\sigma_{_{\text{HON}}},\\
\Psi(\Theta_w)&=\frac{1}{4}(2+\cos(\Theta_w))(1-\cos(\Theta_w))^2,
\end{align}
with $\sigma_{_{\text{HON}}}=14\times10^{-3}\,$J/m$^2$ the value estimated by \citeauthor{article:alimi}\cite{article:alimi} for homogeneous nucleation (HON) in the same range of temperatures and, the contact angle $\Theta_w\simeq \pi/2$, so that $\sigma_{_{\text{HEN}}}=\sigma_{_{\text{HON}}}/2^{^{1/3}}$. This is often employed to model heterogeneous nucleation.

Moreover, as can be observed in figure \ref{fig:4}, the ratio $J/J_{_{\text{CNT}}}$ is of order $10^{-2}$  which means the effective attachment rate is $\overline{f(X_N)}\simeq 10^{-2}f(X_N^\star)$ (dividing Eq. (\ref{eq:22-new}) by (\ref{eq:21-new})). Hence, the methodology followed in CNT of approximating $f(X_N)\simeq f(X_N^\star)$  (e.g. p.168 of \citeauthor{book:kashchiev-1}\cite{book:kashchiev-1}) to obtain $J_{_\text{CNT}}$ lead to an overestimation of the nucleation rate by several orders of magnitudes.

The model also provides the necessary tools for calculating the induction time statistics by considering the Kramers law (Eqs. (\ref{eq:20}) and (\ref{eq:27})). The same set of twelve experimental values for induction time were used in order to calculate the experimental cumulative distribution function (blue triangles in figure \ref{fig:5}). Using the estimated value of $\eta^{\text{eff}}$ into equation (\ref{eq:20}), $P(t_{\text{ind}}\leq t)$ (solid orange line in figure \ref{fig:5}) was obtained. As can be observed, the theoretical prediction does not fit as expected to the experimental statistics. Such a disagreement between the predicted and the observed curves is not due to the assumptions made to integrate the escape rate but to the fact that in experiments we truncated the statistics to the longest time observed, $t_{\text{max}}^{\text{exp}}$, i.e. no nucleation events were recorded after an arbitrary time $t_{\text{max}}^{\text{exp}}$ corresponding to the duration of the experiment.
 So the experimental $t_{\text{ind}}$ values were biassed towards small values. Using a SINT (Euler-Maruyama) to simulate with equation (\ref{eq:4}) a finite number of assays with an upper limit equal to the experimental observation time (i.e. if a simulation exceeds the upper limit $t_{\text{max}}$, then such assay is not considered and another one starts) the simulated results (green squares in figure \ref{fig:5}) are much closer to the experimental curve.

\section{Conclusions}
In this work a stochastic differential equation was presented as the equation of motion of the classical order parameter in the classical nucleation theory, i.e. the cluster radius. This SDE for modelling the time evolution of the radius of an individual cluster was used to obtain theoretical equations to predict (nucleation) induction times and its statistics. These equations are applicable to the often used crystallization setups in which supersaturation changes over time and can be used to identify the underlying mechanism of cluster formation by fitting the measured nucleation rates to equations (\ref{eq:19-new})-(\ref{eq:20-new}) and (\ref{eq:22-new}), as well as to predict induction times (Eqs. (\ref{eq:indtime}) and (\ref{eq:24-2-new})) and induction time distributions (Eqs. (\ref{eq:19})-(\ref{eq:20}) and (\ref{eq:24})-(\ref{eq:27})).
We present here not only an application of the model but also a method to obtain
theoretical
and simulated predictions of both induction times (and hence of nucleation rates) and cumulative distribution functions at different concentrations and temperatures.
The first test of the stochastic model against experimental data reveals its potential ability for calculating nucleation rates and induction time statistics. Both the analytical and the numerical results predicted by the Brownian model seems to be in good agreement with the experimental data. Hence, the interpretation of the nucleation stage as an escape process could be an optimal tool to study deeper
problems of the first-order phase transitions. However the theoretical approximations usually offers worse results than the simulations of the SDE by using  stochastic integrators due to the experimental truncation of the statistics to the longest time observed. Additional work is in progress to deduce analytical expressions for nucleation rate, induction times and induction time distributions when neither stationary nor semiadiabatic limit can be considered.

\section*{Acknowledgements}

This research was supported by Ministerio de Ciencia e Innovaci\'{o}n,
FPI grant BES-2010-038422 (project AYA2009-10655). The authors wish to thank I. Rodr\'{i}guez-Ruiz for providing us unpublished results. We are also immensely grateful to Prof. J.F. Lutsko, Dr. J.M. Delgado-L\'{o}pez, Dr. A.E.S Van Driessche and Dr. J.A. Gavira for their comments on an earlier version of the manuscript.


\section{Introduction}
\label{Mot}

The first generation of galaxies is expected to have reionized cosmic
hydrogen by a redshift $z\approx 7$ \citep{Loeb13}. One of the most
important unknown parameters regulating the way in which reionization
proceeds is the fraction of ionizing photons, $\fesc$, that escape
outside the virial radius of the galaxies in which they were produced
and into the surrounding intergalactic medium (IGM).  The growth of
ionized regions with cosmic time depends on the average value of this
parameter over viewing angle per galaxy and over the galaxy population
\citep{Wyithe03, Trac07, Bouwens12, Finkelstein12}.

Observations of galaxies at $z= 1-3$ indicate a broad set of escape
fraction values ranging from a few percent to tens of percent
\citep{Steidel01, Fernandez03, Shapley06, Inoue06, Siana07,
Giallongo08, Iwata09, Boutsia11, Grazian11, Vanzella12a,
Vanzella12b}. This potentially reflects strong variations of the
escape fraction with viewing angle and evolutionary time.

A calculation of the escape fraction of ionizing radiation from first
principles is difficult, as it depends on the spatial distribution of
ionizing sources relative to the neutral hydrogen in galaxies and is
therefore sensitive to the small scale clustering of young stars and
the interstellar medium.  Due to their higher density, high-redshift
galactic disks are expected to have small escape fractions during
reionization \citep{Wood00,Gnedin08, Razoumov06, Wise09}. Even for
local disk galaxies, \citet{Dove00} concluded that $\fesc$ should not
exceed a few percent in most cases, as a result of efficient radiation
trapping by the shells of the expanding superbubbles around OB star
associations.

Since the mean density of the Universe (and hence galaxies) scales as
$(1+z)^3$, theoretical calculations tend to conclude that $f_{\rm
esc}$ should decrease with increasing redshift.  This expectation
appears to be in conflict with recent reionization models and data
\citep{Bolton07, Salvaterra11, Mitra12a, Mitra12b, Kuhlen12,
Finkelstein12} that point towards the need for an increasing UV
emissivity towards high redshifts. This can be achieved in two
ways. The first is to postulate an increasingly top-heavy initial mass
function (IMF) of stars that would increase the number of photons
emitted per baryon incorporated in stars.  Although possible, this
explanation seems to be disfavored by the lack of observational
evidence for an early population of substantially more massive stars
\citep{Cayrel04, Caffau11} and by recent numerical simulations
\citep{Greif12}. A more appealing solution is represented by a
possible increase in the average escape fraction $\langle f_{\rm
esc}\rangle$ during reionization, recently suggested by \cite{Mitra12b}.
Here we explore a novel physical explanation for this unexpected trend.

We show that the lowest-mass galaxies near or below (so-called
mini-halos) the hydrogen cooling threshold have high escape fractions
but their contribution to the UV production rate decreases with cosmic
time due to large-scale radiative feedback processes during
reionization that either photo-heat or sterilize them by dissociating
their H$_2$.  However, internal feedback is of foremost importance as
well: as these small systems achieved large $\fesc$ by rapidly (within
$\sim 10^4$ yr) ionizing their entire interstellar medium, this gas
will become loosely bound to the host galaxy. Under these conditions
it is difficult to sustain a continuous mode of star formation (SF),
particularly because following the death of massive stars powerful
supernova explosions will clear the gas out of the potential well
before a stable disk structure is established. More massive galaxies
will form a disk with the net result of trapping most of their
ionizing photons.

If high-redshift dwarf galaxies form stars predominantly in episodic
bursts, as suggested by the high incidence of galaxy mergers at
increasing redshift \citep{Munoz11}, the conventional argument that
their star formation efficiency is suppressed by supernova (SN)
feedback may not be valid. In particular, SN feedback is inhibited if
the duration of the starburst is shorter than $\sim 10~{\rm Myr}$, the
lifetime of SN progenitors. In this regime, the global radiative
feedback sets the amount of gas initially available for making stars
and hence the overall efficiency of star formation, $f_\star$.  In our
calculations, we will calibrate the value of $f_\star$ based on recent
measurements of the cosmic stellar mass density at high redshifts
\citep{Gonzalez11}.

For an instanteneous starburst, $f_{\rm esc}$ can be large only if an
escape route for ionizing photons is opened within a few Myr, prior to
the death of the massive stars that produce these photons. This
process cannot be mediated by supernova explosions, which occur after
the emission of ionizing photons has already started to decline.  The
key question is therefore whether the surrounding blanket of absorbing
hydrogen atoms can be ionized before the massive stars in the
starburst end their life. We address this question within the context
of a simple model in which we populate dark matter halos with
primordial gas that cooled to the temperature floor of atomic
hydrogen, $\sim 10^4$K.

Our basic point is simple: massive galaxies exist at all redshifts and
the escape fraction of ionizing radiation from them is calculated and
observed to be small; however, the lowest-mass galaxies ($\Tvir \simlt
10^4$ K) exist only before and during reionization and have a high
$\fesc$.  Our assumption of a short-lived star formation phase without
SN feedback is self-consistent, as we will show that photoionization
heating of the gas suppresses star formation in dwarf galaxies after a
period of time much shorter than the lifetime of massive stars.

In \S 2 we describe the details of our calculation and in \S 3 we
present our numerical results. Finally, we discuss our main
conclusions in \S 4.  Throughout the paper, we adopt the WMAP7 set of
cosmological parameters with $\Omega_m = 0.27$, $\Omega_{\Lambda} = 1
- \Omega_m$, $\Omega_b h^2 = 0.023$, $h=0.71$, $\sigma_8=0.81$, $
n_s=0.97$, and ${\rm d} n_s/{\rm d} \ln k =0$ \citep{Larson11}.

\section{Method of Calculation}
\label{Met}

We derive the escape faction of the ionizing photons, $\fesc$, based
on a set of simplifying assumptions.  We consider a spherical dark matter
halo of mass, $M$, virializing at redshift $z$ and characterized by an
internal density (spherically averaged) profile following the
Navarro-Frenk-White (\cite{Navarro97}, NFW) form,
\be
\rho(r) = \frac{\rho_c
\delta_c }{cx(1+cx)^2},
\label{eq1}
\ee
where $x \equiv r/r_{vir}$, $r_{vir}$ is the virial radius of the
system, $c$ is the halo concentration parameter, $\delta_c =
(18\pi^2/3)c^3/F(c)$ is a characteristic overdensity, and \be F(t) =
\ln(1+t) - \frac{t}{1+t}.
\label{eq2}
\ee For the concentration parameter we use the results of
\cite{Prada12}, extrapolating them towards lower masses and higher
redshift when necessary. The circular velocity, $v_c^2(r) = GM(r)/r$,
can be expressed as \be v_c^2 = V_c^2 \frac{F(cx)}{xF(c)},
\label{eq3}
\ee with $V_c^2 \equiv GM/r_{vir}$. Baryons fall into the dark matter
halo potential well are shock-heated to the virial temperature $\Tvir
= (\mu m_p/2k) V_c^2$, and then relax to a hydrostatic configuration
\citep{Makino98}, \be \ln\rho(r) = \ln \rho_0 -
\frac{V^2}{V_c^2}[v_e^2(r=0)-v_e^2(r)].
\label{eq4}
\ee
\begin{table}
\begin{minipage}{85mm}
\caption{Values of central number density, $n_0$ (in units of cm$^{-3}$)  for dark matter halos of different
total mass, $M$, at $z=9$ and for different values of the gas temperature, $T$, expressed in terms of
the virial temperature $\Tvir$ through the parameter $V^2=\Tvir/T$.}
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline\hline
$\log_{10}\Tvir$ [K]  &             4.31     &         4.98        &         5.65             &          6.32                  \\ \hline
$M$ [$\msun$]      &      $10^8$ &      $10^9$ &     $10^{10}$   & $10^{11}$            \\ \hline
$V^2=1.0$                &   20.6  &         23.9         &         36.0        &          167.8             \\ \hline
$T =2\times 10^4$ K   &   24.3          &        $1.9\times 10^4$         &       $4.3\times 10^6$                 &       ---                    \\ \hline
\end{tabular}
\end{center}
\label{Tab1}
\end{minipage}
\end{table}
Here, $v_e$ is the halo escape velocity from the halo,
\be v_e^2(r) = 2 V_c^2 \frac{F(cx)+
cx/(1+cx)}{xF(c)};
\label{eq5}
\ee $\mu=1.22$ is the mean molecular weight of a neutral primordial
H-He gas; $k_B$ is the Boltzmann constant and $m_p$ the proton
mass. We have also inserted in equation (\ref{eq4}) the extra factor
$V^2=V_c^2/c_s^2 = \Tvir/T > 1$, where $c_s^2 = 2k_BT/\mu m_p$ is the
effective sound speed at a generic temperature $T\neq \Tvir$. The
standard case in which $T=\Tvir$ is obtained by setting $V=1$. By
manipulating equation (\ref{eq4}) we then get the final expression for
the gas density profile as a function of radius, \be \rho(x) =\rho_0
e^{-\gamma^2 V^2} (1+cx)^{\gamma^2 V^2/cx} \equiv \mu m_p n_0 \eta(x);
\label{eq6}
\ee where we have defined $\gamma^2=2c/F(c)$, and the ancillary
function $\eta(x) = e^{-\gamma^2 V^2} (1+cx)^{\gamma^2 V^2/cx}$.  The
central density, $\rho_0$, is obtained by requiring that the total gas
mass in the halo is equal to the cosmological value, i.e. $M_b =
(\Omega_b/\Omega_m) M = f_b M$. This procedure yields \be \rho_0 =
(18\pi^2/3) f_b c^3 e^{\gamma^2 V^2} \left[ \int_0^c dt
(1+t)^{\gamma^2 V^2/t} t^2\right]^{-1} ,
\label{eq7}
\ee Table 1 provides the values of the central number density, $n_0$,
for different halo masses and $V$ values at a fixed redshift,
$z=9$. As expected, $n_0$ increases in larger masses even in the case
$V=1$ in which the gas is at the virial temperature of the halo. This
reflects the strong dependence of $n_0$ on the concentration
parameter, which increases from $c\approx 4$ for $M=10^8 \msun$ to
$c\approx 6.5$ for $M=10^{11} \msun$. This increasing trend is
amplified if the gas is cooler than the virial temperature ($V>1$) or
if it is thermostated at a fixed temperature (such as $T=2\times 10^4$
K in Table 1), leading to unrealistic densities for the most massive
halos with $M=10^{11} \msun$ (which represent rare $3\sigma$
fluctuations of the density field at $z=9$). For massive halos we
obtain artificially high densities because we ignore angular
momentum. Once the gas condenses to a sufficiently small scale, its
rotation will halt contraction and it will settle to a disk. Disks
were considered by \cite{Wood00} and \cite{Dove00} who independently
concluded that the escape fraction was negligible for a smooth gas
distribution in high-redshift disks and also local spirals. We
therefore argue that once the gas contracts enough to make a disk, its
density will not increase by as much as our spherical model predicts
but UV photons will not be able to escape from it anyway, based on
earlier studies.

\subsection{Ionization front evolution}
Having specified the gas density distribution inside halos, we may now
derive the evolution of an ionization front (IF) driven by the
emission of photons with energy $>13.6$eV from a burst of star
formation at the halo center.  We assume that a fraction $f_\star$ of
the available gas, $f_b M$, of a given halo is instantaneously
converted into stars distributed in mass according to a Salpeter IMF
in the range ($m_{low}, m_{up}$) = ($1 \msun, 100 \msun$) and
with absolute metallicity $Z=10^{-3}$.
The ionizing photon production rate, $Q(t)$, by the stellar cluster
can be computed exactly from population synthesis models: we use here
\texttt{Starburst99}\footnote{http://www.stsci.edu/science/starburst99/}
\cite{Leitherer99}. The time dependence of the production rate of
ionizing photons rate under these conditions is,
\be Q(t) =
\frac{Q_0}{1+(t/t_0)^4} ,
\label{eq8}
\ee with ($Q_0, t_0$) = ($10^{47} \rm s^{-1} \msun^{-1}, 10^{6.6}$
yr). Equation (\ref{eq8}) illustrates the important point that after
$\sim 4$ Myr, the production rate of ionizing photons rapidly drops
as a result of the death of short-lived massive stars. This has
important implications for $\fesc$ as discussed hereafter. Note that
equation (\ref{eq8}) implies that the number of ionizing photons
emitted per baryon incorporated into stars is $N_\gamma \approx 0.5 Q
t_0 m_p/\msun =5\times 10^3$.

The time evolution of the IF radius, $r_I$, is described by an
ordinary differential equation that expresses a detailed balance between
the ionization and recombination rates within the volume enclosed by the
IF, \be \frac{dr_I}{dt} = \frac{1}{4\pi n_H r_I^2} \left[ Q(t) -
\frac{4\pi}{3}r_I^3 n_H^2 \alpha^{(2)}\right],
\label{eq9}
\ee where $n_H=0.92 n$ is the hydrogen density for a primordial gas
and $\alpha^{(2)} = 2.6\times 10^{-13} (T/10^4{\rm K})^{-1/2}$ is the
Case B recombination rate of hydrogen \citep{Maselli03}. Normalizing
by an effective recombination rate in the halo, $ {\cal{\dot N}}_{rec}
= (4\pi/3) \alpha^{(2)} n_{H0}^2 r_{vir}^3$, and adopting a
dimensionless time variable $\tau=t/t_{rec}=\alpha^{(2)} n_{H0} t$,
equation (\ref{eq9}) can be written in a dimensionless form,
\be \frac{dx}{d\tau} = \frac{Q(t)}{3 {\cal{\dot N}}_{rec}
\eta(x) x^2} - \frac{1}{3} x \eta(x).
\label{eq10}
\ee Equation (\ref{eq10}) describes the expansion of the \HII region
into the stratified gas distribution within the halo. As the density
decreases outwards, the IF accelerates and eventually exits the virial
radius at a time $t_{out} \equiv t(x=1)$.  Until $t_{out}$ all the
ionizing photons are absorbed inside the halo\footnote{Here we ignore
the possibility of a highly inhomogeneous medium in which low density
channels for escape exist.}; hence, $\fesc\approx 0$.  However, for $t
> t_{out}$, ionizing photons will be only used \emph{to keep} the halo
ionized by balancing recombinations within it, but a large fraction of
them will be able to escape into the IGM, thus making $\fesc>0$. We
then express the net escape fraction as follows: \be \fesc=
\frac{\int_{\tau(x=1)}^\infty d\tau Q(t)}{\int_{0}^\infty d\tau Q(t)}
.
\label{eq11}
\ee As we will see shortly, though, either the IF is efficiently
confined by the density, ${dx}/{d\tau} \rightarrow 0$, or the blow out
will occur on a very short time scale thus making $\fesc \approx 1$.
We will refer to these two different situations as a ``confined'' or
``unconfined'' IF in the rest of the paper.

\section{Results}
\label{Res}
We solved numerically equation (\ref{eq10}) for a number of halo
masses with $V=1$ in the range $M=10^{8}$--$10^{11} \msun$ at two
selected redshifts, $z=9$ and $z=14$.  Figures \ref{Fig1} and
\ref{Fig2} show the results for $f_\star=0.2$\%, a reasonable value for
for nearly primordial star-forming halos \citep{Barkana01, Ciardi05, Pawlik12};
note that \cite{Wise09} found that $\langle \fesc f_\star \rangle =0.02$ averaged over all
atomic cooling ($\Tvir \ge 8000$ K) galaxies assuming a Salpeter IMF.
\begin{figure}
\includegraphics[width=90mm]{Figure_01.png}
\caption{Time evolution of the ionization front radius, $r_I$, at
$z=9$ within halos of different total mass $M=10^{8.+j/3} \msun$, with
$j=0,..,9$ from the top curve ($M=10^8 \msun$) to the bottom one
($M=10^{11} \msun$).  The assumed star formation efficiency is
$f_\star=0.002$ and the gas temperature factor $V=({\Tvir/T
})^{1/2}=1$.  }
\label{Fig1}
\end{figure}

\begin{figure}
\includegraphics[width=90mm]{Figure_02.png}
\caption{Same as Fig. \ref{Fig1}, but for $z=14$.
}
\label{Fig2}
\end{figure}
IFs in these early halos evolve rapidly to the strong radiative input
by the stars and reach a radius $r_I = r_{vir}$ within $\sim 1$ Myr,
yielding $\fesc = 0.98-0.90$ for halo masses up to $M\approx 2\times
10^{10}\msun$, beyond which the IF is confined within the halo, and
hence $\fesc=0$. It requires fine tuning for a single object to achieve
intermediate $\fesc$ values; rather, the process favors on-off
states implying that the escape fraction of photons from galaxies can be
intermittent: depending on variations of $f_\star$ in different bursts
(say, due to the increasing gas metallicity or decreasing gas fraction
due to mass loss from supernova-driven winds) the same halo might
switch from being dark to bright above the Lyman limit.  This is consistent
with the results of \cite{Shapley06} who showed that the PDF of the escape fraction
from $z=3$ star-forming galaxies is bimodal (see their Fig. 5).  Note also
that for the largest IF-confined halos, the \HII region starts to
shrink at $t\simgt 3$ Myr as a result of recombinations and the
decline in the supply rate of ionizing photons.

The situation is dramatically different at $z=14$ (for the same
$f_\star$): at that epoch all halos confine their IF and have
$\fesc=0$. This is due to the increase in the gas density and the
larger concentration factor at a fixed mass. A more global view of the
dependence of $\fesc$ on halo mass and star formation efficiency at
$z=9$ is given by Fig. \ref{Fig3}. The abrupt transition from $\fesc
\approx 1$ to zero is very weakly dependent on halo mass, except for
the narrow range of $-2.8 < \log f_\star < -2.2$ in which larger halos
tend to be more opaque to ionizing photons. Above this range, $\fesc =
1$ independently of halo mass and star formation efficiency, and below
this range, $\fesc=0$.
\begin{figure}
\includegraphics[width=90mm]{Figure_03.png}
\caption{Map of $\fesc$ values as a function of halo mass and star
formation efficiency at $z=9$ assuming $V=({\Tvir/T })^{1/2}=1$. The
color scheme defined by the vertical bar on the right hand side.}
\label{Fig3}
\end{figure}

Figure \ref{Fig4} shows the examples with $V >1$ introduced in Table
1, for which $T=2 \times 10^4$ K, corresponding to $V^2 = (1.04, 4.82,
22.39, 103.96)$ for halos of mass $M/\msun=(10^8, 10^9, 10^{10},
10^{11})$ and hence resulting in considerably denser and more compact
gas configurations. The IF front is now largely confined within all
halos with $M>10^8 \msun$ even if the star formation efficiency has
been increased by a factor of 50 with respect to the case with $V=1$
shown in Fig.  \ref{Fig1}. In addition, the IF escapes the virial
radius in a very short time\footnote{We have checked that the IF
speed is subluminal, $dr_I/dt < c$}, $t_{out}\sim 10^4$ yr.

The different behavior relative to the case of $V=1$ is caused by the
associated high gas density (see Table 1). The small scale height of
the gas distribution, correspondingly reduced by a analogous factor $
\propto V^2$, forces the steep acceleration seen just prior to blow
out.  In analogy with the $V=1$ case, we show the dependence of
$\fesc$ on halo mass and star formation efficiency at $z=9$
(Fig. \ref{Fig5}).  In this case the $\fesc=1$ area is much narrower
and limited to small halos with $M< 10^{8.7} \msun$ and $f_\star
\simgt 10^{-3}$.  Thus, according to our model, it is only the
smallest halos that are truly able to inject ionizing photons in the
IGM and therefore substantially contribute to reionization. Larger
halos are instead able to confine their \HII regions very
effectively. At $z=14$ the $\fesc = 1$ region is qualitatively the
same, but its area shrinks because of the larger background density
and the increased concentration parameter at a fixed halo mass.

 To check the effects of the central density we have performed two additional experiments: (a) we have introduced a central density core of size $x_{core} = r_{core}/r_{vir}$ such that the value the density in the core is equal to $\rho(x < x_{core}) = \rho(x_{core})$; (b) in a more extreme case, we have completely removed the gas, i.e. $\rho(x < x_{core}) = 0$ within the core. The two cases lead to similar conclusions: appreciable deviations from the solutions presented in Fig 4 of the paper are seen only for values of $x_{core} \simgt 0.1$. Stated differently, it would be necessary to completely remove all the gas within 10\% of the virial radius to produce appreciable deviations from what we have obtained with our density profile.

The existence of a minimum value of $f_\star \approx 10^{-3}$ to allow
the escape of photons can be understood based on a simple argument.
This minimum number of ionizing photons provided to each hydrogen atom
in the halo must balance the number of recombinations it undergoes
during the lifetime of the starburst, approximately equal to $t_0 =
10^{6.6}$ yr (see equation \ref{eq8}). This can be expressed by the
following inequality: \ba f_\star &\ge & \frac{C}{{\cal
N}_\gamma}\frac{\alpha^{(2)} t_0}{\mu m_p} (18\pi^2)\Omega_b \rho_c
(z) \\ & = &0.7\times 10^{-3} \left(\frac{C}{3}\right)
\left(\frac{5\times 10^3}{{\cal N}_\gamma}\right)
\left(\frac{1+z}{10}\right)^3 ,
\label{eq11a}
\ea where $C=\langle n_H^2\rangle/\langle n_H\rangle^2\approx 3$ is
the clumping factor accounting for the gas density structure inside a
NFW halo.

The results for halos with $M\simlt 10^8 \msun$ in which $\Tvir \simlt
2 \times 10^4$ K at $z=9$ need some extra attention.  In fact, by
fixing the gas temperature to $2\times 10^4$K, the gas might have
never been accreted in the first place; this is the case if for
example the halo is located within an already ionized region. To
correct this problem, we have set $V=1$ in these mini-halos; of
course, this requires that some cooling agent, such as \HH, HD or LiH,
would cool the gas in these systems to make stars.
\begin{figure}
\includegraphics[width=90mm]{Figure_04.png}
\caption{Same as Fig. \ref{Fig1}, but for halos with mass in the range
$M=10^{8.+2j/30} \msun$, with $j=0,..,9$ from the top curve ($M=10^8
\msun$) to the bottom one ($M=3\times 10^{8} \msun$). The star
formation efficiency is $f_\star=0.1$ and the gas temperature has been
fixed to $T=2\times 10^4$ K, giving a variable $V$ temperature factor
(see text for details).  }
\label{Fig4}
\end{figure}
\begin{figure}
\includegraphics[width=90mm]{Figure_05.png}
\caption{Same as Fig. \ref{Fig3} but for a fixed gas temperature
$T=2\times 10^4$ K.  }
\label{Fig5}
\end{figure}

\begin{table*}
\begin{minipage}{170mm}
\caption{Values of various quantities related to the global evolution of the mass-averaged escape fraction. See text for more details.}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline\hline
&\multicolumn{7}{c}{Mass-averaged $\fesc$ --  Constant $f_\star$} \\
             & $M_{sf}$ & $M_{f1}$               & $M_{f0}$      &  $f_{coll}(> M_{sf})$   & $f_{coll} (>M_{f0})$   & $\langle\fesc\rangle$  &  $f_\star$  \\
            & [$\msun$] &  [$\msun$]  &[$\msun$]  &    &                                    &                \\\hline
$z=6$  &2.977e+08   & 2.977e+08                 &     4.852e+08      &   1.379e-01           &    1.262e-01      & 8.478e-02  & 3.000e-02          \\
$z=9$  & 1.169e+07 & 1.169e+07 & 2.190e+08 & 1.130e-01 & 5.936e-02 & 4.747e-01 & 3.000e-02 \\
$z=9$  & 1.169e+07 & 1.169e+07 & 1.253e+08 & 1.130e-01 & 6.897e-02 & 3.897e-01 & 3.000e-03\\
$z=9$  & 1.169e+07 & 1.169e+07 & 5.818e+08 & 1.130e-01 & 4.374e-02 & 6.130e-01 & 1.000e+00 \\
$z=9$  & 1.169e+07 & 1.169e+07 & 1.169e+07 & 1.130e-01 & 1.130e-01 & 0.000e+00 &  1.000e-04\\
$z=14$ & 1.634e+06 & 1.634e+06 & 8.131e+07 & 5.373e-02 & 1.471e-02 & 7.263e-01 & 3.000e-02 \\\hline\hline
&\multicolumn{7}{c}{Mass-averaged $\fesc$ -- Variable $f_\star(z)$}\\
             & $M_{sf}$ & $M_{f1}$               & $M_{f0}$      &  $f_{coll}(> M_{sf})$   & $f_{coll} (>M_{f0})$   & $\langle\fesc\rangle$  &  $f_\star$  \\
            & [$\msun$] &  [$\msun$]  &[$\msun$]  &    &                                    &                \\\hline
$z=8$   & 5.361e+07 & 5.361e+07 & 2.019e+08 & 1.098e-01 & 8.315e-02 & 2.426e-01 & 8.119e-03\\
$z=10$ & 3.664e+06 & 3.664e+06 & 1.119e+08 & 1.102e-01 & 5.217e-02 & 5.267e-01 & 4.799e-03\\
$z=12$ & 2.030e+06 & 2.030e+06 & 7.128e+07 & 8.003e-02 & 3.148e-02 & 6.066e-01 & 4.281e-03\\\hline
\end{tabular}
\end{center}
\label{Tab2}
\end{minipage}
\end{table*}

\subsection{Global redshift evolution}
Next we would like to constrain the redshift evolution of $\fesc(z)$
when averaged over the entire galaxy population.

The dark matter halo mass function, $n(M,z)$, is well described by the
\cite{Sheth02} form: \be
n(M,z)dM = A\left(1+\frac{1}{\nu'^{2q}}\right) \sqrt{\frac{2}{\pi}}
\frac{\bar\rho}{M} \frac{d\nu'}{dM}
\exp\left(-\frac{\nu'^{2}}{2}\right) dM,
\label{eq12}
\ee with $A=0.322$, $q=0.3$, $\nu'=\sqrt{a}\nu$, $a=0.707$,
$\nu=\delta_c/D(z)\sigma(M)$, with $\delta_c=1.686$, $D(z)$ being the
linear growth factor and $\sigma(M)$ being the r.m.s.  amplitude of
density fluctuations on a mass scale $M$.

We define the mass-averaged escape fraction, $\langle\fesc\rangle$ at
a given $f_\star$ and redshift $z$ as follows: \be
\langle\fesc\rangle(z, f_\star) =1 - \frac {F\left
[>M_{f0}(z,f_\star)\right ]}{F\left[>M_{sf}(z)\right]} ,
\label{eq13}
\ee where $F(>M,z)=\int_M^\infty nM dM/\int_0^\infty nM dM$ is the
collapse fraction of all matter in haloes with mass exceeding $M$ at redshift
$z$.

Equation (\ref{eq13}) contains two important characteristic
masses. The first, $M_{f0}(z,f_\star)$, is defined as the halo mass
above which, for a given redshift $z$ and star formation efficiency
$f_\star$, no ionizing photons can escape, i.e. $\fesc=0$. The results
presented above allow us to precisely determine $M_{f0}(z,f_\star)$.

The second mass scale, $M_{sf}(z)$, denotes the redshift-dependent
minimum mass of star-forming haloes which is set by radiative
feedback. The evolution of $M_{sf}$ is determined by two distinct
radiative feedback processes.  The first one is related to the
increase of the cosmological Jeans mass in progressively ionized
cosmic regions; as a result, the infall of gas in haloes below a given
virial temperature is quenched. The evolution of $M_{sf}(z)$
depends on the details of the reionization history \citep{Gnedin00,
Okamoto08, Schneider08}. Guided by these findings, we adopt the the
value $M_{sf}=M(\Tvir=T_2)$ with $T_2=3\times 10^4$ K after the end of
reionization, here assumed to occur at $z_{rei}=6$.

Before reionization (at least in the mostly neutral regions) a second
type of feedback, involving the photodissociation of \HH molecules by
the Lyman-Werner (LW) background photons, becomes important. As H$_2$
is the primary coolant for minihalos ($\Tvir< 10^4$ K), fragmentation
of the gas to stars could be suppressed in these objects through H$_2$
dissociation by the UV background \citep{Haiman96, Ciardi00,
Kitayama00, Machacek01}.  According to \cite{Dijkstra04}, $z\sim 10$
mini-halos with $\Tvir = T_1 \approx 2\times 10^3$ K can self-shield
H$_2$ and cool, hence we use this value as the lower threshold for
star formation.  During reionization the interplay between the
photo-ionization heating and LW-feedback is complicated (see the
discussion in \cite{Mesinger08} and Fig. 25 of \cite{Ciardi05}) and
and the evolution of $M_{sf}$ during this epochs is uncertain. To
circumvent this uncertainty, we follow the phenomenological approach
of \cite{Salvadori09} and \cite{Salvadori12} who suggested that the
metallicity distribution and iron-luminosity relation of the Ultra
Faint dSphs in the Milky Way halo can be explained if these objects
are relic mini-halos that formed their stars prior to
reionization. With this perspective, one can reconstruct $M_{sf}(z)$ by
fitting local data. This approach yields the functional form, \be M_{sf}(z) =
M[\Tvir = (T_2-T_1) e^{(z-z_{rei})^2/\Delta z}+T_1] ,
\label{eq14}
\ee with $\Delta z=4$, which we adopt here.

In general $M_{f0}(z,f_\star) > M_{sf}(z)$. However, for sufficiently
low values of $f_\star$, $M_{f0}(z,f_\star) = M_{sf}(z)$ and
$\langle\fesc\rangle =0$.  To illustrate this point (see Table
\ref{Tab2}), we examine the results at $z=9$ in Figure \ref{Fig5}.
For $f_\star=0.03$, stars are formed only in halos more massive than
$M_{sf} = 1.17\times 10^7 \msun$. Ionizing photons can only escape from
halos less massive than $M_{f0} = 2.19\times 10^8 \msun$. As a result,
there is a narrow range of objects that contribute to the ionizing
photons budget in the IGM. Below that range stars cannot form and
above it the gas is sufficiently dense to confine the \HII
region. When weighted with the halo mass function, though, the
collapsed mass fraction above $(M_{sf}, M_{f0})$ is $(0.113, 0.059)$,
yielding the large value $\langle\fesc\rangle(z=9, f_\star=0.01) =
0.847$. Since \HII regions in abundant low mass halos are unbound, the
average escape fraction is high. The width of the range progressively
shrinks as star formation becomes less efficient.  This is evident
from Table \ref{Tab2}, where for $z=9$ and $f_\star = 10^{-4}$,
$M_{f0}(z,f_\star) = M_{sf}(z)$ and the escape fraction drops to
zero. These results underline the importance of the possible evolution
in the star formation efficiency, which we address below.
\begin{figure}
\includegraphics[width=90mm]{Figure_06.png}
\caption{Redshift evolution of the mass-averaged escape fraction for
different assumptions concerning the star formation efficiency. Solid
lines refer to the redshift-independent values $f_\star=0.03, 0.003$
as indicated by labels; dotted lines are the same cases in which a fixed
$\fesc=0.05$ has been assumed for large halos with $M>M_{f0}$ (see Sec. \ref{Con}); the shaded area refers to $f_\star(z)$ (shown
by the shaded cyan area and multiplied by a factor 100) derived from
the stellar mass density evolution by \cite{Gonzalez11} based on equation
(\ref{eq15}).}
\label{Fig6}
\end{figure}

The redshift evolution of the results can be understood in terms of a
balance between two opposite trends. On the one hand, the halo central
densities tend to become larger with increasing redshift, providing a
more efficient confinement of the IF and acting to decrease
$\langle\fesc\rangle$. On the other hand, more dwarf galaxies and
mini-halos are able to produce ionizing photons as the critical mass
$M_{sf}$ decreases with increasing $z$. The overall evolution for the
two specific values $f_\star = (0.03, 0.003)$ can be seen in Figure
\ref{Fig6} (solid curves). For the high $f_\star$ case, the escape
fraction continues to increase with redshift, reaching
$\langle\fesc\rangle = 0.83$ at $z=20$. If instead $f_\star = 0.003$,
the escape fraction reaches a peak value of $0.59$ at $z=11.5$ where
it starts to decrease until it reaches zero at $z=17.5$ when the
density effect dominates over the enhanced photon production rate.
Interestingly, a nearly redshift-independent $f_\star$ above a
threshold of $\approx 0.1$\% (see Fig. \ref{Fig5}) always leads to
$d\langle\fesc\rangle/dz>0$ during the final phase of reionization.
\begin{figure}
\includegraphics[width=90mm]{Figure_07.png}
\caption{Cumulative fraction of ionizing photons produced by halos
below a given halo mass $M$. Lines refer to a constant star formation
efficiency $f_\star=0.003$; shaded regions refer to the upper limit of
the redshift-dependent $f_\star$ (equation \ref{eq15}) based on
\cite{Gonzalez11} SMDs. For both sets of curves, solid (dashed) lines
refer to $z=8$ ($z=10$); for $f_\star=0.003$ the case for $z=12$
(dotted) is also shown.  }
\label{Fig7}
\end{figure}
Although a rather weak mass dependence of $f_\star$ for galaxies of
stellar mass $M_\star > 3\times 10^7 \msun$ is inferred from local
SDSS data \citep{Dekel03}, there are indications of a decreasing star
formation efficiency for high-$z$ galaxies. This is in line also with the
conclusion (based on numerical simulations) by \cite{Munshi12} that low
mass ($M < 10^9 M_\odot$) halos have generally lower (0.1\%) values
of $f_\star$. Recently, different
groups \citep{Stark09, Labbe10, Gonzalez11} have obtained stellar
masses of galaxies in the redshift range $4 < z < 7$ from spectral
energy distribution fitting to rest-frame optical and UV fluxes from
Hubble-WFC3/IR camera observations of the Early Release Science field
combined with the deep GOODS-S Spitzer/IRAC data. From this data it
has been possible to reconstruct (see Fig. 4 of \cite{Gonzalez11}) the
evolution of the comoving Stellar Mass Density (SMD), which we define
here as $\rho_\star(z)$. The main result is that $\rho_\star(z=6)=10^7
\msun $ Mpc$^{-3}$; in addition, the SMD grows with decreasing
redshift as $\propto (1 + z)^{-3.4 \pm 0.8}$. We can exploit this
observational result to obtain a robust estimate for $f_\star(z)$
based on the relation,
\be
f_\star(z) = \frac{\rho_\star(z)}{F\left[>M_{sf}(z)\right] \Omega_b\rho_c} .
\label{eq15}
\ee
In other words, we tune the star formation efficiency so that
baryons in collapsed halos account for the abundance of stars observed at
$z<7$. Note that the \cite{Gonzalez11} SMD has also been corrected for
incompleteness at low masses.  The range of allowed values of
$f_\star$ implied by the data and obtained through equation
(\ref{eq15}) are plotted as a shaded area in Fig. \ref{Fig6}. The data
suggest $f_\star = 0.003$ at $z=6$ (for display reasons we have
multiplied the values by a factor 100 in the Figure) decreasing
thereafter by a factor of about three above $z=10$.

The main consequence of a suppressed star formation efficiency in low
mass objects is to introduce a new intermediate mass scale
$M_{f1}$. Because of their lower star formation efficiency, halos above $M_{sf}$
might not produce enough stars to ensure $\fesc=1$. IFs become unconfined only for
somewhat larger halos of characteristic mass  $M_{f1}$, such that $M_{sf} <
M_{f1} < M_{f0} $. Thus  $\fesc > 0$ only in the interval
[$M_{f1},M_{f0}$]. Stated differently, even objects that are above $M_{sf}$ and can
form stars, are now more efficient at confining their IFs; therefore
their escape fraction becomes zero. Table \ref{Tab2} shows how
$M_{f1}$ gradually approaches $M_{f0}$ towards high redshifts, closing
the non-vanishing escape fraction window.

Another interesting question is how the global production of ionizing
photons is distributed among halos of different masses. We quantify
this through the cumulative fraction of ionizing photons produced by
halos below a given mass, \be f_\gamma(< M, z) \propto \int_{M_{sf}}^M
\fesc(M', z, f_\star ) f_\star(M_{sf},z ) M' dM',
\label{eq16}
\ee as shown in Figure \ref{Fig7}.  Due to their large escape
fraction, halos with masses $M\simlt 10^8 \msun$ are the dominant
contributors to the cosmic ionizing photon budget. This holds for the
more realistic redshift-dependent $f_\star$, but also for the constant
$f_\star$ case. Thus, the low-mass galaxies and mini-halos might have
dominated reionization. In particular, at $z=10$ the mass
corresponding to $T=10^4$ K limiting the mini-halo mass range from
above, is $\sim 3\times 10^7 \msun$. Figure \ref{Fig7} shows that
independently of the star formation efficiency prescription, we
predict that mini-halos contribute $\sim 40\%$ of the total
emissivity.

Reionization models suggest an increasing emissivity towards high
redshifts \citep{Bolton07, Mitra12a, Mitra12b, Kuhlen12}, whose most
natural explanation is in terms of a $\fesc$ increasing with redshift
as we find here. For example, \cite{Mitra12b} find that: {\it (i)} the
escape fraction increases from $\langle\fesc\rangle =
0.068_{-0.047}^{+0.054}$ at z=6 to $0.179_{-0.132}^{+0.331}$ at $z=8$;
and {\it (ii)} at $z=10$ $\langle\fesc\rangle > 0.146$.

\section{Discussion}
\label{Con}

We have calculated the escape fraction of ionizing photons from
starburst galaxies during reionization, assuming that galactic gas
cools down to a temperature floor ($\sim 10^4$K for atomic hydrogen
and $\sim 2\times 10^3$K for self-shielded molecular hydrogen).  We
have found that most of the escaping ionizing photons originates in
the lowest mass galaxies with $T_{\rm vir} \simlt 10^4~{\rm K}$.  Due
to the shallowness of the gravitational-potential well in these
galactic halos, the gas maintains a low density and its recombination
rate is low. This allows the ionization front to break through the
virial radius rapidly, resulting in a high escape fraction of the
emitted ionizing radiation. As reionization progresses, the
intergalactic medium is photo-heated and the abundance of galaxies
with $T_{\rm vir}\sim 10^4$K declines, leading to a decrease in
$\langle f_{\rm esc}\rangle$.  This physical mechanism may account for
the increase of $\langle f_{\rm esc}\rangle$ with redshift inferred by
\cite{Mitra12b}.

Our study is consistent with previous studies of low mass galaxies by
\cite{Wise09}, who found a high escape fraction ($> 0.1$). Low mass
galaxies do not form a disk \citep{Pawlik11}, which could trap
effectively ionizing photons. The relative importance of minihalos for
the global production rate of ionizing photons is determined by the
product $\fesc \times f_\star$. Here we have used a phenomenological
calibration of $f_\star$. In order to produce the observed stellar
density at $z=6$--$7$ by minihalos, $f_\star$ has to be $\simgt
0.001$. This value lies in the middle range resulting from
\cite{Wise09} idealized halo simulations, but they find the star
formation efficiency in halos with $M< 10^7 \msun$ to be smaller.  We
find that the required increase of $\fesc$ with redshift cannot be
achieved if minihalos are sterile or form stars inefficiently.

Our simplified calculation ignored angular momentum which drives cold
gas to a rotating disk configuration. Our model allowed cold gas to
condense to arbitrarily small scales without making a disk, and so one
might worry that it underestimates $f_{\rm esc}$. However, previous
calculations have demonstrated that on the scale of galactic disks the
escape fraction is very small \citep{Wood00,Gnedin08, Razoumov06,
Wise09}, and so our inclusion of angular momentum could not increase
$f_{\rm esc}$ significantly. In particular, the conclusion that
$f_{\rm esc}$ is high in low-mass halos must be robust.  Although our
simplified model predicts zero escape fractions from brief starbursts
in massive halos, the actual value of $f_{\rm esc}$ would be finite in
reality due to persistent star formation histories and winds driven by
supernova feedback or Ly$\alpha$ radiation pressure
\citep{Dijkstra08}.  Our results are in qualitative agreement with
\cite{Gnedin08} and \cite{Razoumov06} who simulated more massive
galaxies ($M \simgt 10^{10} \msun$) and inferred an extremely small
escape fraction values $\fesc\simlt 1$\%), even in situations where
the star formation rate is high ($\sim 10 \msun$ yr$^{-1}$). This
results from the high gas densities in massive galaxies due the
presence of a rotationally-supported, geometrically-thin disk.

To get a simple estimate of the impact of these larger galaxies,
we have also run cases allowing a certain ionizing photon fraction
(here assumed to be 5\%, guided by local observations) to escape from
halos of mass $M > M_{f0}$.
The results are shown for two cases in Fig. 6 (dotted lines).
The reported trend of an increasing $\fesc$ for $z<11$ remains clearly visible.

Finally, we note that the trend we find for the increase of the escape
fraction with redshift should be enhanced if the initial mass function
of stars was tilted towards more massive stars. This is because a
higher production rate of ionizing photons would allow the
ionization front to ionize the surrounding hydrogen blanket even
faster than we calculated.

\section*{Acknowledgments}
Our computations were based on the Scipy open source scientific tools
for Python \texttt{http://www.scipy.org/}. This work was supported in
part by NSF grant AST-0907890 and NASA grants NNX08AL43G and
NNA09DB30A.

\bibliographystyle{apj}