\section{Introduction} \label{sec:intro}
Much of the current interest in cosmological hydrogen
reionisation may be distilled to two related questions: (1) When did it
occur?; and (2) What sources dominated the ionising photon budget?
Theoretical efforts to address the first question have traditionally
considered two observational constraints: the optical depth to Thomson
scattering of cosmic microwave background (CMB) photons, $\tau_{\mathrm{es}}$, and the
opacity of the Lyman-$\alpha$ forest~\citep[][and references therein]{tra09}.
While these works have shed much light on the nature of inhomogeneous
reionisation, their predictions have until recently been uncertain owing
to the unknown relation between ionising luminosity and halo
mass~\citep{mcq07}.  With the recent
installation of the Wide-Field Camera 3 on board the Hubble Space
Telescope, a new generation of observations is now constraining the
abundance and colors of galaxies back to $z=8$ and
beyond~\citep{fin10,bou11,dun12,gon11,gra11,mcl11,oes12}.  These
observations may be used to constrain the relationship between star
formation rate and halo mass via abundance-matching
studies~\citep{tre10,mun11}.  Furthermore, the relationship between ionising
luminosity and star formation rate is now constrained by the direct
detection of escaping Lyman continuum flux from high-redshift
galaxies~\citep{sha06,sia10,nes11}.  Hence there are now tenuous observational
links connecting halo mass with ionising luminosity, a major step toward
understanding the relationship between high-redshift galaxies and their
environment.

By combining observations of galaxies, the CMB, and the Lyman-$\alpha$
forest, it is possible to ask whether the observed ionising
sources and reionisation history are
compatible~\citep[for example,][]{haa12,kuh12}.
Two uncertainties hinder progress toward answering this question.
The first remains the fraction of ionising photons $f_{\mathrm{esc}}$ that escape
into the IGM: direct constraints are only available at lower
redshifts and higher luminosities than are relevant to reionisation.
Furthermore, $f_{\mathrm{esc}}$ is sensitive to parsec-scale physical processes
that occur deep within galaxies' interstellar media and is difficult
to model numerically~\citep[][and references therein]{fer11}.  The second
uncertainty involves the recombination rate of the reionisation-epoch
IGM.  The unknown ratio of the IGM's true recombination rate to its hypothetical
rate under the assumption of uniform density and temperature is often referred
to as the IGM clumping factor $C$.  We will define $C$ explicitly in
Section~\ref{sec:defs}; for the present the important point is that it is
proportional to the volume-averaged hydrogen recombination rate in the IGM.
High values of $C$ require abundant faint galaxies and/or a large $f_{\mathrm{esc}}$
while low values of $C$ are easier to reconcile with existing observations.
$C$ is sensitive to the IGM's density, ionisation, and temperature fields,
hence an accurate estimate requires three-dimensional numerical simulations.

\begin{figure}
\centerline{
\setlength{\epsfxsize}{0.5\textwidth}
\centerline{\epsfbox{fig.clump_literature.ps}}
}
\caption{Our ionisation- and temperature-corrected IGM clumping
factor at $z=6$ in comparison to previous calculations.  Our work
indicates that the clumping factor lies within the range 2.7--3.3,
which is consistent with the low values that have been found
in recent works.
}
\label{fig:clumplit}
\end{figure}

To this end, $C$ has been studied numerous times using numerical
simulations.  We compile predictions of the clumping factor at $z=6$
from a variety of works in Figure~\ref{fig:clumplit}.
The early work of~\citet{gne97} yielded a clumping factor
of 30--50 (their Figure 2) depending on the choice of definition.  However,
this clumping factor included the recombinations that occur in gas whose
ionisation state is dominated by the local field rather than by the global
extragalactic ultraviolet ionising background (EUVB), hence it is an
overestimate of the recombination rate within the gas
that is conventionally associated with the IGM (see
also~\citealt{gne00} and~\citealt{mcq11}).  For the same
reason,~\citet{tra07} also inferred a clumping factor of roughly 30
from simulations that subtended a significantly larger cosmological
volume.

\citet{ili07} applied this insight to calculate an improved clumping
factor using high-resolution N-body simulations whose spatial resolution
was comparable to the Jeans length in the IGM.  By excluding matter within
virialized halos, they derived a redshift-dependent fitting function that
climbs to $\approx10$ at $z=6$, a value which was also found in the recent
N-body simulations of~\citet{rai11}.  This work confirmed that
physically-motivated
density cuts reduce $C$.  However, as the
authors noted, their calculation did not account for the fact that
photoionisation heating (photoheating) tends to smooth the diffuse IGM.
We will confirm
in Figure~\ref{fig:pdf} that this leads to an overestimated clumping factor.

\citet{paw09a} corrected this problem by running cosmological hydrodynamic
simulations that model reionisation using a spatially-homogeneous EUVB.
They used several plausible density cuts to calculate the IGM clumping
factor as a function of redshift.  We show their predictions assuming that
the IGM corresponds to overdensity cuts of
$\Delta \equiv \rho/\langle \rho \rangle < 50$ (lower range)
and $\Delta < 100$ (upper range).  This work generally confirmed the
impression that considering only low-density gas and accounting for
the tendency for the photoheating of diffuse regions
both reduce $C$.  However, it suffered
from two drawbacks.  First, their use of a spatially-homogeneous EUVB
pressurized regions that should in reality be self-shielded, leading
to uncertainty in the gas density distribution from which the
predicted clumping factors were derived.  Second, their analysis
did not account explicitly for the IGM temperature, which could
suppress the recombination rate further.

In a qualitatively similar work,~\citet{shu12} ran a cosmological
hydrodynamic simulation with a spatially-homogeneous EUVB and
computed the temperature-corrected clumping factor within gas
that satisfied cuts on baryon overdensity, temperature, metallicity,
and ionisation state.  Although these extra cuts render direct
comparisons with~\citet{paw09a} uncertain, the inferred range is
consistent.  Both works are also in reasonable agreement
with~\citet{fau08}, who used the gas density distribution
of~\citet{mir00} to find that $C=$3--4 if the clumping factor is
averaged only over gas with overdensity $\Delta$ less than 67.

The uncertainty in the derived clumping factor may be reduced by
using radiation transport simulations to identify directly the gas
that is self-shielded.  ~~\citet{mcq11} use this approach in
post-processing to derive a clumping factor of 2.4--2.9.  Their
calculations do not account for the IGM temperature and,
like the works of~\citet{paw09a} and~\citet{shu12}, begin with
hydrodynamic simulations in which the moderately-overdense gas is
incorrectly exposed to the EUVB (although they confirm directly
that shielding gas that is dense enough to form stars has no
impact).  Nonetheless, their results are clearly consistent with
a relatively low clumping factor.

The goal of the present work is to combine all of these ideas into a
single radiation hydrodynamic framework that accounts realistically
for the time and spatial dependence of the IGM's density, ionisation,
and temperature fields.  Our fiducial simulation reproduces the
observed UV luminosity function of galaxies at
$z\geq6$~\citep{fin11} and is marginally consistent with the effective
optical depth to absorption of Ly$\alpha$\  photons and $\tau_{\mathrm{es}}$.  Hence, our
estimated IGM recombination rate is constrained by a variety of observations.
As Figure~\ref{fig:clumplit} indicates, we will confirm that effects
that have been neglected in recent works based on hydrodynamic
simulations such as incorrectly heating dense regions and ignoring
the gas temperature are indeed small.  Of course, the
purpose of clumping factors is to inform analytical models of
reionisation~\citep[for example,][]{mad99}.
We will therefore verify that analytical reionisation calculations
incorporating clumping factors do, in fact, reproduce the behavior
of numerical simulations.  This check will assess directly the
magnitude of the systematic uncertainties inherent in calculations
that ignore higher-order effects such as shadowing, source
clustering, and the light-travel time of ionising photons.

Having verified that our analytical model reproduces the behavior of
our numerical model reasonably well, we will use it to ask how
many photons are absorbed by filamentary overdensities in the
diffuse IGM.  This will quantify the error introduced in models
that compute reionisation using IGM density distributions in
which the gas is effectively preheated prior to the onset of
reionisation.

The outline of our paper is as follows: We begin by introducing our
suite of simulations in Section~\ref{sec:sims}.  We compare
different definitions of the clumping factor in Section~\ref{sec:defs}
and explore how feedback effects modulate their values in
Section~\ref{sec:feedback}.  In Section~\ref{sec:photCount}, we
incorporate our derived clumping factors into an analytical model
for reionisation and test how well analytical models of
reionisation reproduce the behavior of numerical calculations.
We apply this model to study how many photons are absorbed in
filaments in Section~\ref{sec:photCost}.  Finally, we discuss
our results in Section~\ref{sec:disc} and summarize in
Section~\ref{sec:sum}.  We introduce two improvements to our
treatment of the Eddington tensors in Appendix~\ref{app:fedd}
and quantify resolution limitations in Appendix~\ref{app:res}.

\section{Simulations} \label{sec:sims}

\subsection{Star Formation and Outflows}
Our numerical methods, input physics, and cosmology are similar
to~\citet{fin11} with the exception of our adopted treatments
for an evolving $f_{\mathrm{esc}}$, subgrid self-shielding (see below),
and two numerical optimizations regarding the treatment of
the Eddington tensor field (Appendix~\ref{app:fedd}).
We model the growth of structure and the feedback processes
that couple galaxies with the IGM using our custom version of the
cosmological galaxy formation code {\sc Gadget-2}~\citep{spr05}.
{\sc Gadget-2} implements a formulation of smoothed particle
hydrodynamics (SPH) that simultaneously conserves entropy and
energy and solves for the gravitational potential with a
tree-particle-mesh algorithm.  Dense gas cools radiatively
using the primordial cooling processes in Table 1 of~\citet{kat96}
and the metal-line cooling tables of~\citet{sut93}, which assume
collisional ionisation equilibrium.  In contrast with~\citet{kat96},
however, we evolve the ionisation states of hydrogen and helium and
the electron abundance simultaneously with the cooling equations
using a nested subcycling approach whose timestep is limited by a
chemical Courant condition; see~\citet{fin11} for details.
We initialize the IGM temperature and neutral hydrogen fraction to
the values appropriate for each simulation's initial redshift as
computed by {\sc recfast}~\citep{won08}, and we assume that
helium is initially completely neutral.  We generate the initial
density field using an~\citet{eis99} power spectrum at redshifts
of 249 and 319 for simulations subtending 6 and 3 $h^{-1}{\rm Mpc}$,
respectively.  All simulations assume a cosmology in which
$\Omega_M=0.28$, $\Omega_\Lambda=0.72$, $\Omega_b = 0.046$,
$h=0.7$, $\sigma_8 = 0.82$, and the index of the primordial power
spectrum $n=0.96$.

Our goal of modeling galaxies and reionisation simultaneously
requires a treatment for the ability of cool gas to form stars.
We adopt the subgrid two-phase interstellar medium treatment
of~\citet{spr03}, which can be tuned to reproduce the observed
relation between the surface densities of gas and star
formation~\citep{ken98}.  The physical density threshold for star
formation is 0.13 cm$^{-3}$.  This value is motivated by observations
of a critical density for the onset of star formation and lies within
the range at which the thermo-gravitational instability is expected
to become active~\citep{sch04}.  Varying it has only a minor impact
on the predicted star formation rate density~\citep{sch10}.
We account for metal enrichment owing
to supernovae of Types II and Ia as well as asymptotic giant branch
stars; see~\citet{opp08} for details.

Allowing cold gas to form stars without any feedback invariably
results in overproducing the observed reionisation-epoch star formation
rate density~\citep{dav06,fin11}.  The accepted solution is to allow
feedback from massive stars to expel star-forming gas from galaxies.
As the spatial resolution necessary to form such outflows
self-consistently is beyond the reach of current cosmological
simulations~\citep{mac99,hop12,pow11}, we impose star formation-driven
outflows by stochastically applying kicks to star-forming gas
particles.  The amount of gas kicked per unit stellar mass formed
and the kick velocities are adjusted to reproduce the expected
scalings from momentum-driven outflows~\citep{mur05}.  We temporarily
disable hydrodynamic forces in outflowing gas in order to
mimic the way in which outflows escape through holes in
higher-resolution simulations~\citep{mac99}.  Hydrodynamic forces are
restored once the gas has traveled for $1.95\times10^{4}$ km s$^{-1}/v$
million years (where $v$ is the kick velocity) or its density drops
below 10\% of the threshold density for star formation.  In practice, this
allows outflowing gas to reach galactocentric radii of 50--100
physical kpc and then re-accrete onto the central galaxy on a timescale
of $\sim1$ Gyr~\citep{opp08}.  The distance that outflows travel before
re-accreting is large compared to the typical radius at which hydrodynamic
forces are enabled, hence the decoupling prescription has little impact
on outflows once they escape the ISM.  The long re-accretion timescale enables
expelled gas to modify the IGM recombination rate and opacity depending
on the density of outflowing gas; exploring these possibilities is one
of the goals of the present work.

\subsection{Radiation Transport}\label{ssec:rt}
We evolve the EUVB on a fixed Cartesian grid using the
moments of the cosmological radiative transfer equation.  We close
the moment hierarchy with a variable Eddington tensor field that is
updated periodically via a time-independent ray-casting calculation.
This approach accounts accurately for the inhomogeneous opacity
field as long as the updates occur frequently
enough~\citep{fin09a}.  Each cell's ionising luminosity is given by
a sum over the luminosities of its star-forming gas particles.
Their luminosities, in turn, depend on their star formation rates
and metallicities, where the metallicity dependence follows an
analytical fit to Table 4 of~\citet{sch03}.

The fraction of
ionising photons that escapes into the IGM in order to participate in
reionisation $f_{\mathrm{esc}}$ is unknown.  Recent work suggests that it must
vary steeply with redshift in order to bring observations of galaxies,
the Lyman-$\alpha$ forest, and the CMB optical depth to Thomson scattering
into agreement~\citep{ino06,fin11,haa12,kuh12}.  We adopt a strongly
redshift-dependent $f_{\mathrm{esc}}$ in our simulations:
\begin{eqnarray}
f_{\mathrm{esc}} = \left\{
\begin{array}{lc}
  f_{\mathrm{esc,5}} \left(\frac{1+z}{6}\right)^\kappa & z < 10 \\
  1.0 & z \geq 10
\end{array} \right.
\end{eqnarray}
Here, the normalization $f_{\mathrm{esc,5}}$ sets the escape fraction at $z=5$
and $\kappa$ controls
how strongly $f_{\mathrm{esc}}$ varies with redshift.  For simulations with outflows,
we use $f_{\mathrm{esc,5}}=0.0519$ and $\kappa=4.8$.  The normalization is consistent
with observations of Lyman break galaxies~\citep{nes11} at $z\sim3$, while the
power-law index lies within
the range (1--6) that~\citet{kuh12} indicate is required by observations
if the UV luminosity function of galaxies extends to absolute magnitudes
of -15 (as required by observations of gamma-ray bursts;~\citealt{tre12}).
With this parametrization, the effective optical depth to Lyman-$\alpha$
absorption at $z=5$ is $\tau_\alpha=3.1$, marginally consistent with the observed
range of 2--3~\citep{fan06}.  Meanwhile, the predicted ionising emissivities
are (4.6, 7.0) $\times10^{50}$s$^{-1}$Mpc$^{-3}$ at $z=$(5, 6).  At $z=5$,
this is consistent with the observed range of $4.3\pm2.6$~\citep{kuh12}.
At $z=6$, however,  it overshoots the observed limit $<2.6$, suggesting
that the predicted emissivity should strengthen rather than declining
toward the end of reionisation.  Note that our model is not unique in
failing to reproduce the observed emissivities at $z=5$ and $z=6$
simultaneously; this evolution is evidently quite strong compared to
expectations from models (see, for example, Figure 10 of~\citealt{kuh12}
or Figures 8 and 15 of~\citealt{haa12}).  The resulting reionisation
history corresponds to an optical depth to Thomson scattering of 0.071,
only slightly below the observed $1\sigma$ confidence intervals of
$\tau_{\mathrm{es}}=0.088\pm0.015$~\citep{kom11}.

Simulations without outflows require a steeper power-law index and a lower
normalization because they overproduce the observed galaxy abundance~\citep{dav06}.
We adopt $f_{\mathrm{esc,5}}=0.0126$ and $\kappa=7.21$.  With this dependence,
they predict emissivities of 5.4 and 7.8 $\times10^{50}$s$^{-1}$Mpc$^{-3}$ at $z=5$
and 6, respectively, once again indicating good agreement at $z=5$
but incorrect evolution to higher redshift.  The predicted Lyman-$\alpha$
optical depth at $z=5$ is 2.2, well within the observed range.
The predicted Thomson scattering optical depth is 0.074, within
the observed $1\sigma$ confidence intervals.

The physical interpretation of an evolving $f_{\mathrm{esc}}$ could take many
forms.  For example, direct measurements of ionising continuum flux
from galaxies at $z\approx3$ suggest that current stellar
population models underestimate the ionising luminosity of
low-metallicity Population II stars~\citep{nes11}.  If confirmed,
this observation could imply an evolutionary trend to higher ionising
luminosity at high redshift that mimics an evolving $f_{\mathrm{esc}}$.
Alternatively, observations indicate that galaxies grow bluer
and by inference less dusty beyond $z=4$~\citep{fin10}; this could
drive up $f_{\mathrm{esc}}$ if dust dominates the absorbing column
(although~\citealt{gne08} argue that $f_{\mathrm{esc}}$ is dominated by
the ISM's geometry rather than its dust content).  Finally,
galactic outflows can fill a halo with gas, potentially
creating a significant optical depth to ionising photons.
If these screens grow weaker with increasing
lookback time, then the fraction of ionising photons that travel
unimpeded from the edge of a galaxy's ISM to the virial radius
could grow accordingly; this would mimic an evolving $f_{\mathrm{esc}}$.
We defer detailed investigation of these possibilities to future
work; for the present, we simply choose a parametrization that
will improve the realism of our simulated density field.
We discuss drawbacks to this approach in
Section~\ref{sec:disc}.

We assume that each photoionisation deposits 4.08 eV of latent
heat into the IGM.  This heats newly-reionised gas to
$\approx15,000$ K, which lies within the range that is expected
for a Population II stellar spectrum.  Each radiation
transport cell's opacity is given by its volume-averaged
neutral hydrogen density (counting only the gas that is not
self-shielded; see below) multiplied by the cross-section
for ionisation at an energy of 17.68 eV.

\subsection{Self-shielding Within Overdense Systems}

Two central challenges in using numerical models to compute the
IGM recombination rate involve deciding which gas belongs to the
IGM and ensuring that overdense gas is allowed to self-shield
against the EUVB.

The need to define the IGM in a physically-motivated way is
illustrated by the results of the pathbreaking radiation
hydrodynamic simulations of~\citet{gne97} and~\citet{gne00}.
These simulations faithfully modelled the growth of ionised
regions as well as the hydrodynamic response of the heated gas.
However, the volume-averaged recombination
rates were computed over all the gas in the simulation, leading to
recombination rates that were much larger than what is plausible for
the moderately overdense ($\rho/\langle \rho \rangle<100$) gas that
is conventionally associated with the IGM.  It also led to a
strong dependence of the recombination rate (and the reionisation
photon budget) on the simulation dynamic range, with higher-resolution
simulations absorbing more photons~\citep{gne00}.

Improvement results from using simple density cuts to isolate
low-density gas in post-processing from precomputed simulations.
This approach has been applied both to N-body
simulations~\citep{ili07,rai11} and to cosmological hydrodynamic
simulations~\citep{paw09a,shu12}.  It reduces $C$ as expected,
but the results remain dependent on the precise cuts used to
define the IGM.  In the case of hydrodynamic simulations, additional
uncertainty results from assuming a homogeneous EUVB, which may
pressurize regions that are in reality dense enough to self-shield
and does not account for the spatial inhomogeneity of reionisation.
Both of these effects could modify the IGM's density and temperature
fields~\citep{fur09}, which in turn impact the recombination rate.

In this work, we combine these two approaches by using radiation
hydrodynamic simulations of the reionisation epoch that incorporate
a physically-motivated treatment for self-shielding within overdense
regions.  Effectively, we incorporate directly into our simulations
the definition that \emph{the IGM consists only of that gas whose
ionisation state is not determined by local sources}; denser systems
are allowed to remain neutral.  By isolating overdense systems from
the EUVB, we (1) improve the accuracy of our predicted density and
temperature fields; (2) prevent our simulations from overestimating
the reionisation photon budget; and (3) obviate the need to choose
a self-shielding threshold overdensity (such as 50 or 100) when
deriving the IGM recombination rate in post-processing: the
threshold evolves on-the-fly in a way that follows the (local)
gas temperature and (global) EUVB amplitude.

We model the self-shielding of overdense regions using a subgrid
recipe because our radiation transport solver does not resolve
them spatially even though our hydrodynamic solver does.  For
example, Lyman limit systems (with neutral column densities of
$\sim10^{17}$cm$^{-2}$) have a characteristic size of 10 physical
kpc~\citep{sch01,mcq11} and must be treated with 5--10 resolution
elements~\citep{gne06}.  In our
higher-resolution simulation (r6n256wWwRT32), the gravitational softening
length is $0.1$ kpc (Plummer equivalent) while the radiation transport solver achieves
a resolution of 38 kpc (where both numbers are in physical units
at $z=6$).  Hence while our hydrodynamic solver certainly resolves
Lyman limit systems, our radiation solver's resolution is a factor
of $\approx20$ too coarse.

We overcome this limitation by following
the argument of~\citet{sch01}.  This work
showed that, if Lyman limit systems are in hydrostatic equilibrium,
then they can be identified directly with gas above a critical
overdensity~\citep[see also][]{mcq11}.  The threshold, $\Delta_{\rm{lls}}$,
is given by:
\begin{equation}\label{eqn:dlls_typ}
\Delta_{\rm{lls}} = 25
\left(\frac{T}{10^4 \rm{K}}\right)^{0.3}
\left(\frac{1+z}{7}\right)^{-3}
\left(\frac{\Gamma_{\rm{HI}}}{2\times10^{-13}}\right)^\frac{2}{3}
\end{equation}
at the Lyman limit.  Here we have assumed that hydrogen is fully
ionised and helium is singly-ionised, as expected for soft
ionising sources~\citep{cia12,fri12}.  $\Gamma_{\rm{HI}}$ represents the volume-averaged
EUVB and does not include the influence of local sources.
This estimate is within a factor of two of what others have
found~\citep[for example,][]{bol07}.  We do not expect slight
differences owing, for example, to our assumed recombination
rate to affect our results.  We enforce $\Delta_{\rm{lls}}\geq2$ in order to
preserve self-shielding at $z\gg10$; this floor is not expected to
affect our results.  We use this approximation to model subgrid
self-shielding as follows: For each gas particle whose density lies
below the physical threshold for star formation (0.13 cm$^{-3}$),
we compute $\Delta_{\rm{lls}}$ using the local gas temperature, the case-A
recombination rates, and the global mean ionisation rate per
hydrogen $\Gamma_{\rm{HI}}$.  We then assume that the optical depth to the
diffuse background varies as a power-law of the overdensity
\begin{eqnarray}\label{eqn:tau_delta}
\tau_\Gamma = \left\{
\begin{array}{lc}
  0 & \Delta < \Delta_{\rm{lls}} \\
  \tau_0\left(\frac{\Delta}{\Delta_{\rm{lls}}}\right)^b & \Delta \geq \Delta_{\rm{lls}}
\end{array} \right.
\end{eqnarray}
In the optically-thin limit, $\tau_0 = 1$ and $b=1.5$~\citep{sch01}.
In reality, $\tau_\Gamma$ increases more rapidly with density because
the neutral
fraction in the partially shielded outskirts of an absorber is larger
than in the optically-thin case.  We have experimented with both
$b=1.5$ and $b=3$ and found little difference in practice.  We adopt
the fiducial choices $\tau_0 = 1.0$ and $b = 3$ in this
work.

We implement Equations~\ref{eqn:dlls_typ}--\ref{eqn:tau_delta} directly
into our simulations so that overdense gas is shielded from the
EUVB on-the-fly.  Our use of different discretizations for the
radiation and hydrodynamic solvers requires that we do so in two
ways.  When updating the ionisation state of a given particle, we
attenuate the radiation field of its host cell by the factor
$e^{-\tau_\Gamma}$, where $\tau_\Gamma$ is evaluated using the global
EUVB and the particle's temperature.  Conversely, when assembling the
opacity field on the radiation solver's grid, we reduce each particle's
contribution to its host cell's opacity by the same factor.  We include
the opacity of each gas particle's self-shielded region by treating it
as an opaque sphere with volume equal to $fx_{\mathrm{H\,I}}(1-e^{-\tau_\Gamma})\Delta V$,
where $\Delta V$ is the ratio of the gas particle's mass to its
density, $x_{\mathrm{H\,I}}$ is its neutral fraction, and $f$ is a parameter
that is tuned via high-resolution calculations to $1/8$.  For example,
if $e^{-\tau_\Gamma}=0.5$ for a particle, then its ionisation rate is
50\% of that of its host cell's, and only 50\% of its neutral
hydrogens contribute to the cell's opacity.  We also augment the cell's
opacity by the amount
$4.84(0.5x_{\mathrm{H\,I}}\Delta V)^{-2/3}/A_{\rm{cell}}$, where $A_{\rm{cell}}$
is the cell's area.  Note that the opacity from
partially-neutral gas particles whose overdensity falls below the
self-shielding threshold is included directly in the radiation
transport solver as in our previous work.

Intuitively, this treatment divides self-shielded regions into an optically
thin skin and an optically-thick core.  The overdensity range over which
a region transitions from optically-thin to optically-thick follows
from Equation~\ref{eqn:tau_delta}.  The opacity owing to the
optically-thin region is distributed uniformly throughout the
radiation transport cell while the core is treated as a photon sink.
This approximation preserves the ability of dense gas
to remain self-shielded within a coarse radiative transfer grid.
It will enable radiative hydrodynamic simulations to subtend
cosmological volumes even as they treat the radiation field's small
scale structure with sufficient detail to model observables such
as the abundance of low-ionisation metal absorbers and the
post-reionisation IGM opacity.

In reality, regions that are dense enough for star formation to occur
($\Delta>1000$) contain a mixture of ionised and neutral
gas owing to the local radiation field from massive stars.
They remain entirely neutral in our simulations owing to our self-shielding
prescription.   Accounting for their recombinations would lead to much
larger clumping factors~\citep{gne97,gne00}.  However, the corresponding
absorptions should be identified with an ionising escape fraction $f_{\mathrm{esc}}$
that is less than unity rather than an increased clumping factor
because they are not dominated by the EUVB.  We encapsulate their
impact via an assumed $f_{\mathrm{esc}}\leq1$ that is tuned using observational
insight (see above) and defer the challenging task of modeling $f_{\mathrm{esc}}$
directly to future work.

\begin{figure}
\centerline{
\setlength{\epsfxsize}{0.5\textwidth}
\centerline{\epsfbox{fig.dlls.ps}}
}
\caption{The threshold overdensity for self-shielding $\Delta_{\rm{lls}}$
as a function of redshift within our fiducial simulation at a
characteristic temperature of $10^4$K.  $\Delta_{\rm{lls}}$ peaks at over
200 during the heyday of reionisation before shrinking to
roughly 150 by $z=5$ owing to declining escape fractions.
}
\label{fig:dlls}
\end{figure}

Figure~\ref{fig:dlls} shows how $\Delta_{\rm{lls}}$ evolves in our fiducial
simulation for a characteristic temperature of $10^4$K.  Note
that this figure is illustrative; in reality our simulations
compute this threshold on-the-fly using the local temperature.
$\Delta_{\rm{lls}}$ grows from values near unity at the onset of reionisation
to values near 200 by the time that reionisation
has completed.  Gas that is denser than this threshold sees an
attenuated EUVB, suppressing both the reionisation photon budget
and the volume-averaged recombination rate.

In detail, the
values in Figure~\ref{fig:dlls} may be somewhat overestimated
(see, for example,~\citealt{mcq11}).  This occurs because our
monochromatic radiation transport solver uses a smaller
photoionisation cross section $\sigma$ than the value at the
Lyman limit $\sigma_{\rm{LL}}$ in order to match the assumed
photoheating rate of 4.08 eV per photoionisation.  In particular,
$\Delta_{\rm{lls}}$ varies with the photoionisation cross section $\sigma$
as $(\sigma/\sigma_{\rm{LL}})^{-2/3}$~\citep{sch01}, hence our
simulations may overestimate $\Delta_{\rm{lls}}$ by a factor of $\approx1.7$.
This is roughly the difference between our value at $z=6$ and
that of~\citet{mcq11}.  The overall impact of this limitation on
our inferred value of $C$ is weak; we will quantify it in
Section~\ref{sec:disc}.

Table~\ref{table:sims} introduces our simulation suite.  The simulation
names encode the physics and numerical resolutions.  For example,
the r6n256wWwRT16 simulates a $6h^{-1}{\rm Mpc}$ volume (r6) with $2\times256^3$
particles (n256), includes outflows (wW) and evolves the EUVB on
a Cartesian grid with $16^3$ cells (wRT16).  The first four simulations
explore the sensitivity of gas clumping to photoheating
and outflows.  The next three explore convergence with respect to our
radiation and gas solvers as well as the cosmological volume.  The
final simulation assumes an optically-thin EUVB for comparison with
previous work.  Throughout this work, we will use the
r6n256wWwRT16 simulation as our fiducial case.

Our simulations allow us to explore how the IGM's temperature impacts
its recombination rate within the context of a model that treats
photoheating from Population II stars self-consistently
(Section~\ref{ssec:feedbackAll}).  As a check on whether our predictions
are realistic, we compare the temperature at mean
density to recent observations.~~\citet{bol12} used measurements
of the Doppler widths of Ly$\alpha$ absorption lines along the line
of sight to seven quasars to infer that the IGM temperature at mean
density is $\log(T)=3.85\pm0.08$ at $z=6.08\pm0.33$ (their Table 3,
fiducial model).  Our fiducial simulation predicts
$\log(T,z=6)=3.93\pm0.10$.  Here, we report the median $\log(T)$ over
particles with overdensity $-0.2<\log(\Delta)<0.2$; switching
from the median to the mean changes results by $\sim0.1\%$
because the scatter in the temperature distribution is not large.
This comparison is of course incomplete because it does not account
for a variety of observational systematics.  Nonetheless, the fact
that the simulation reproduces the observed value within the errors
indicates that the predicted temperature fields are plausible.

\section{What is gas clumping?} \label{sec:defs}
The motivation for defining clumping factors comes from a desire to take the
volume-weighted mean of the ionisation rate equations.  For example, the
rate of change of the neutral hydrogen abundance $n_{\mathrm{H\,I}}$ is given by
\begin{equation} \label{eqn:dnHIdt}
\frac{dn_{\mathrm{H\,I}}}{dt} = \Gamma_{\rm{HI}} n_{\mathrm{H\,I}} + k_2 n_{\mathrm{H\,II}} n_{\mathrm{e}} - k_1 n_{\mathrm{H\,I}} n_{\mathrm{e}},
\end{equation}
where $\Gamma_{\rm{HI}}$ is the photoionisation rate per hydrogen atom, $k_1$ is the
collisional ionisation rate, and $k_2$ is the (case A) recombination rate.
Each term
on the right side is nonlinear, hence, in principle, it is not possible to compute their
spatial averages unless the cross-correlation of the abundances of $n_{\mathrm{H\,I}}$ and
$n_{\mathrm{H\,II}}$ with each other as well as with the radiation and temperature fields
are known (since $k_1$ and $k_2$ are temperature-dependent).  It is convenient
to encapsulate this subgrid information with clumping factors.  For
example, the spatially averaged recombination rate is approximated as

\begin{equation} \label{eqn:Cdef}
\langle k_2 n_{\mathrm{H\,II}} n_{\mathrm{e}}\rangle \equiv C_{{\rm H\,II,T_b}}\langlen_{\mathrm{H\,II}}\rangle\langlen_{\mathrm{e}}\rangle\langle k_2 \rangle,
\end{equation}
where angle-brackets indicate averages over the entire simulation volume.
Similar clumping factors could be defined for each of the terms on the
right-hand side of Equation~\ref{eqn:dnHIdt}.  In practice, however, analytical
calculations bypass the need for the first term's clumping factor by
assuming that the volume-averaged $\Gamma_{\rm{HI}}$ equals the ionising emissivity
divided by the mean hydrogen abundance.  Meanwhile, the third term can
be neglected because $\Gamma_{\rm{HI}} \gg k_1n_{\mathrm{e}}$ except within shocks.  This means
that, along with the ionising escape fraction, $C_{{\rm H\,II,T_b}}$ is one of the major
uncertainties that hamper efforts to connect the observed abundance of
ionising sources (such as Lyman Break galaxies) with the ionisation state
of the IGM (see~\citealt{koh07} for an expanded discussion
of the other clumping factors).  In Section~\ref{sec:photCount}, we will
show that the numerical
results may be reproduced by an analytical model that considers only the
recombination clumping factor even though the full numerical calculation
accounts for all three terms.

\begin{figure}
\centerline{
\setlength{\epsfxsize}{0.5\textwidth}
\centerline{\epsfbox{fig.clumpdefs.ps}}
}
\caption{The evolution of four different clumping factors in our
fiducial simulation volume.  The clumping factors of all gas whose
density exceeds the star formation threshold reaches 15 by $z=6$.
Accounting for the fact that overdense gas is self-shielded and that
ionised gas is somewhat warmer than $10^4$K leads to much lower clumping
factors (see text for discussion).
}
\label{fig:clumpdefs}
\end{figure}

Equation~\ref{eqn:Cdef} shows that the clumping factor depends on
the density, ionisation, and temperature fields.  Partial clumping
factors are often used when one or more of these fields is
unavailable.  In Figure~\ref{fig:clumpdefs}, we show how the choice
of definition impacts the clumping factor's value in our fiducial
simulation.  The black solid curve shows the
evolution of the clumping factor averaged over all gas whose density
lies below the adopted threshold for star formation (0.13 cm$^{-3}$),
$C_{\mathrm{b}}\equiv\langle\rho^2\rangle/\langle\rho\rangle^2$.  It increases
monotonically as overdense regions collapse, with a brief enhancement around
the reionisation redshift.  The enhancement occurs when the EUVB grows strong
enough to photoevaporate gas from photosensitive
halos~\citep{fin11}, temporarily boosting the mass of gas whose
density lies just below the threshold for star formation.  The
clumping factor of all gas (including star forming gas) decreases at
this epoch as expected, but the clumping factor of gas that lies
outside of galaxies---that is, $C_{\mathrm{b}}$---increases until the
newly-pressurized gas has expanded into hydrostatic equilibrium.

The clumping factor within ionised regions---which is the important
quantity for reionisation---is expected to be different from $C_{\mathrm{b}}$
because reionisation begins in overdense regions and progresses
into voids~\citep{fur04,mcq07} in such a way that the last regions to
remain neutral are overdense~\citep{mir00}.  Furthermore, gas that is
self-shielded does not contribute to recombinations.
To explore these factors, we compute the clumping factor of electrons
and protons
\begin{equation}\label{eqn:chii}
C_{{\rm H\,II}}\equiv \langlen_{\mathrm{e}}n_{\mathrm{H\,II}}\rangle/\langlen_{\mathrm{e}}\rangle\langlen_{\mathrm{H\,II}}\rangle.
\end{equation}

$C_{{\rm H\,II}}$ can be computed in either of two ways.  One way is to compute the
averages over all ionised gas.  Unfortunately,
this introduces uncertainty owing to the choice of ionisation threshold.
For example, averaging only over gas particles that are more than 50\%
ionised gives a slightly different answer than averaging only over
particles that are more than 95\% ionised because the simulated IGM is
not a strictly two-phase medium.  The alternative is to take averages
over the entire simulation volume and then weight by the ionised volume
fraction $x_{\mathrm{H\,II,V}}$.  This can be understood by recalling that, in a
partially-ionised universe with a uniform density, $C_{{\rm H\,II}}=1/x_{\mathrm{H\,II,V}}$.
Therefore, $\xhiivC_{{\rm H\,II}}$ indicates the recombination rate in units of
the recombination rate of a homogeneous universe that has the same
ionisation state, which is the goal of the clumping factor.

We show this quantity with a red dotted line.
The tendency for overdense regions
to remain neutral suppresses $\xhiivC_{{\rm H\,II}}$ below $C_{\mathrm{b}}$ at
all times, with reionisation-epoch values that lie below $\approx6$.
This difference immediately reveals the importance of self-shielding,
which creates a boundary between the ionised IGM and the neutral or
locally-ionised
regions near halos.  It can be modeled directly by spatially resolving
physical scales of $\sim1$ kpc~\citep{gne06}.  However, this incurs
significant computational expense, limiting the simulation's
cosmological volume.  Our approach of coupling a somewhat coarse
grid for the radiation solver with a subgrid prescription for
self-shielding may open up the possibility of simulating cosmological
reionisation within large volumes while treating the IGM's thermal
and ionisation states faithfully.

$C_{{\rm H\,II}}$ is an improved description of the volume-averaged
recombination rate over $C_{\mathrm{b}}$, but it remains incomplete because
the recombination rate $k_2$ also depends on the gas temperature.
To illustrate this, we use a blue dashed curve to show the
temperature-corrected clumping factor $x_{\mathrm{H\,II,V}} C_{{\rm H\,II,T_b}}$, computed
following Equation~\ref{eqn:Cdef} (referred to by~\citealt{shu12}
as $C_{\rm{RR}}$).  This curve remains near zero until
reionisation is well under way because the recombination rate per
proton is spatially anti-correlated with the abundance of ionised
gas.  Once $x_{\mathrm{H\,II,V}}$ exceeds $\approx0.5$, it begins to climb
because $\langle k_2\rangle$ is no longer suppressed by cold neutral
regions.  However, it does not reach $C_{{\rm H\,II}}$ owing to the lingering
presence of self-shielded regions. Following reionisation,
$x_{\mathrm{H\,II,V}} C_{{\rm H\,II,T_b}}$ climbs slowly because the volume-weighted mean
temperature is supported by the slow photoheating of filaments,
which suppresses $\langle k_2 \rangle$.  It never exceeds 3,
reflecting the tendency of photoheating to suppress
recombination rates.

$\xhiivC_{{\rm H\,II,T_b}}$ provides a reasonable description of the clumping
factor within ionised regions once reionisation is well under
way, but it is less informative at early times ($x_{\mathrm{H\,II,V}} < 0.1$)
when the high recombination rate in the predominantly cold, neutral
IGM suppresses $C_{{\rm H\,II,T_b}}$ despite the high clumping factor of
ionised regions.  We may compute a more informative clumping factor
by re-evaluating $C_{{\rm H\,II,T_b}}$ using only those regions (that is, SPH
particles) whose ionised mass fraction is greater than 0.95.
This clumping factor (green dot-dashed) is somewhat larger at earlier
times because it is not suppressed by the high recombation rates
per proton in regions that have not yet been heated.
It merges with $\xhiivC_{{\rm H\,II,T_b}}$ once
$x_{\mathrm{H\,II,V}}\sim1$.  A drawback of $C_{{\rm H\,II,T_b}}$ is that its value depends
on the threshold ionisation fraction that is used to compute it.

While $C_{{\rm H\,II,T_b}}$ emerges naturally from a spatial average of
Equation~\ref{eqn:dnHIdt}, it is difficult to compare with observations
because the mean IGM temperature (and hence $\langle k_2 \rangle$)
at $z=6$ remains poorly-constrained; this is of course equally true
even if $C_{{\rm H\,II,T_b}}$ is computed only over regions whose ionisation fraction
exceeds some threshold.  To resolve these problems, we also compute
the ``observational temperature-corrected clumping factor of ionised
gas"
\begin{equation}\label{eqn:chiiten}
C_{{\rm H\,II,}10^4{\rm K}} \equiv \frac{\langle n_{\mathrm{e}} n_{\mathrm{H\,II}} k_2(T)\rangle}
	{\langle n_{\mathrm{e}} \rangle \langle n_{\mathrm{H\,II}} \rangle k_2(10^4\mathrm{K})}
\end{equation}
in which we replace the volume-averaged recombination rate
$\langle k_2\rangle$ in the denominator with the recombination rate at
$10^4$K.  This clumping factor (referred to by~\citealt{gne00} as $C_{\rm{HII}}$)
expresses the mean recombination rate
without requiring knowledge of the topology of reionisation or the
temperature of ionised gas.  It is slightly higher than $C_{{\rm H\,II,T_b}}$
following reionisation, reflecting the tendency for gas to cool below
$10^4$ K once reionisation is complete (that is,
$\langle k_2 \rangle > k_2(10^4\mathrm{K})$ for $z<6$).  The function
\begin{equation} \label{eqn:cfit}
\xhiivC_{{\rm H\,II,}10^4{\rm K}} = 9.25 - 7.21 \log_{10}(1+z)
\end{equation}
fits its evolution to within 30\% for $z\leq15$ and 10\% for $z\leq 10$.

Note that uncertainties in the recombination rate translate directly
into uncertainties in the clumping factor.  For example, our adopted
(case A) recombination rate is $\sim20\%$ higher than that of~\citet{hui97}
in the temperature range $10^3$--$10^5$K.  This means that ionisation
fronts do not penetrate as far into overdense regions in our simulation
as they would if we assumed the~\citet{hui97} rates, reducing the
inferred clumping factor of ionised gas by the same ratio.  The
resulting uncertainty does not affect our results and is trivial to
adjust for: Simply divide our inferred clumping factors by the ratio
of the alternative recombination coefficient to our
adopted coefficient at $10^4$ K.

In summary, Figure~\ref{fig:clumpdefs} indicates that the clumping
factor's value depends rather strongly on its definition.  Taking
self-shielding into account suppresses it from $\approx30$
to $\approx6$ at $z=5$ because the most overdense regions are self-shielded
and do not contribute to the IGM recombination rate.  Taking the IGM's
temperature into account further suppresses clumping to 2--4
because photoionised regions are generally warmer than $10^4$ K.

\section{The Role of Feedback}\label{sec:feedback}

Galactic outflows and photoheating both affect
gas clumping.  Heating suppresses gas clumping by smoothing density
fluctuations on scales smaller than the local Jeans length~\citep{efs92}.
Outflows increase IGM clumping partly by boosting the amount of
dense gas that resides within halos but outside of galaxies and
partly (for a given $f_{\mathrm{esc}}(z)$) by delaying reionisation and heating.
The latter effect is a straightforward consequence of suppressing
star formation.  The former is expected based on our
previous finding that outflows boost the number of ionising photons
required to achieve reionisation~\citep{fin11} although it may in
reality be weaker if ejected gas self-shields.  In this section, we
quantify how these processes modulate gas clumping.

\subsection{Averaging over the IGM} \label{ssec:feedbackAll}
\begin{figure}
\centerline{
\setlength{\epsfxsize}{0.5\textwidth}
\centerline{\epsfbox{fig.clumpfeedback.ps}}
}
\caption{The impact of photoheating and outflows on the
overall IGM baryonic clumping factor $C_{\mathrm{b}}$.  All simulations subtend
a $6h^{-1}{\rm Mpc}$ volume discretized with $2\times256^3$ particles (r6n256);
see Table~\ref{table:sims}.  Outflows have negligible impact on $C_{\mathrm{b}}$
in the absence of an EUVB, but when acting in concert with an EUVB
their impact is stronger than uncertainties owing to limitations in
our particular treatment of the EUVB.
}
\label{fig:cb_feedback}
\end{figure}

We begin by discussing the separate impacts of outflows and
photoheating on the IGM clumping factor.  In order to
facilitate this discussion, we revert to the definition $C_{\mathrm{b}}$,
which uses all gas that is too diffuse to form stars
(Section~\ref{sec:sims}) without reference to
its ionisation state or temperature.  While this definition
includes gas that is dense enough to self-shield and is therefore
an imperfect estimate of the IGM recombination rate, it allows us
to compare simulations with and without an EUVB.

We use solid black and dotted blue curves in
Figure~\ref{fig:cb_feedback} to show how $C_{\mathrm{b}}$ evolves in simulations
without and with outflows in the absence of an EUVB.  They are
nearly coincident, indicating that outflows affect the gas density
distribution only weakly if they do not couple to an EUVB.  Adopting
a spatially-homogeneous~\citet[][hereafter HM01]{haa01} EUVB dramatically
suppresses the clumping factor following $z=9$ (long-dashed green;
see also~\citealt{paw09a}).

Our self-consistent simulation omitting outflows (r6n256nWwRT16, short-dashed
red) predicts a reionisation history in which the neutral hydrogen fraction
drops to 50\% at $z=9.2$, quite similar to the redshift at which the HM01
EUVB reionises the Universe ($z=9$).  Prior to $z=5$, its $C_{\mathrm{b}}$ exceeds
the HM01 case owing to its more extended reionisation history.  By $z=5$,
its $C_{\mathrm{b}}$ converges to the HM01 case although it shows several spikes
that are likely associated with photoevaporation of star-forming gas
from low-mass ($M_h < 10^9M_{\odot}$) haloes in regions that are just being
reionised.  This suggests that $C_{\mathrm{b}}$ is weakly sensitive to the
reionisation history at early times but converges to expectations
from calculations using a spatially-uniform EUVB fairly rapidly
(probably on the sound-crossing time of the filaments).

Turning to our fiducial simulation (magenta dot-dashed), which does include
outflows, we find that its $C_{\mathrm{b}}$ is higher than the no-wind case at all
redshifts.  This model has significantly less star formation than the
no-wind model, but it also assumes a higher ionising escape fraction in
order to compensate (Section~\ref{ssec:rt}).  Combining these factors leads
to a reionisation history in which the neutral fraction drops below 50\%
at $z=8.9$, only slightly later than the other models.  Despite
their similar reionisation histories, the wind model's $C_{\mathrm{b}}$ is roughly
30\% higher at $z=5$.  Hence although the reionisation histories of the
wind and no-wind models are not exactly the same,
Figure~\ref{fig:cb_feedback} confirms that outflows enhance $C_{\mathrm{b}}$ by
boosting the amount of gas that is overdense but not star-forming.

We draw several conclusions from Figure~\ref{fig:cb_feedback}.  First,
outflows do not affect $C_{\mathrm{b}}$ unless they are heated by an EUVB.  Intuitively,
outflows consist largely of material entrained from their host galaxies'
interstellar media, hence they remain relatively cold and dense unless they
are heated by an EUVB.  Second, an EUVB heats and evaporates gas out of
shallow potential wells such as filaments and minihalos,
suppressing $C_{\mathrm{b}}$.  Finally, a better understanding of the nature of
galactic outflows is {\it at least} as important to our understanding
of the IGM density structure as the nature of the EUVB.  In
Section~\ref{ssec:feedbackRhomax}, we will ask whether this
uncertainty impacts the recombination rate of diffuse gas.

\subsection{The recombination rate as a function of density} \label{ssec:feedbackRhomax}
\begin{figure}
\centerline{
\setlength{\epsfxsize}{0.5\textwidth}
\centerline{\epsfbox{fig.pdf.ps}}
}
\caption{(top)
The gas density PDF in units of mass-weighted probability
per unit $\log_{10}(\Delta)$ for six different simulations as indicated by
the legend in the middle panel at $z=6$; all simulations subtend volumes
with side
length $6h^{-1}{\rm Mpc}$ and resolve halos of total mass $1.4\times10^8M_{\odot}$.
The vertical dashed line indicates the threshold density for star and
wind formation (0.13cm$^{-3}$).  Photoionisation heating pressurizes
gas with densities $\Delta>10$, and outflows transfer gas from
galaxies into the IGM.
(middle)
The cumulative baryonic clumping factor as a function of maximum
overdensity compared to the fiducial simulation of~\citet{paw09a}
(black dotted).  An EUVB suppresses the clumping of all gas,
with minor differences owing to the details of outflows and
reionisation.
(bottom)
The neutral hydrogen fraction as a function of overdensity.
Gas above a threshold overdensity (which is related to $\Delta_{\rm{lls}}$) is
neutral owing to our self-shielding model.
}
\label{fig:pdf}
\end{figure}

Figure~\ref{fig:cb_feedback} suggests that outflows and photoheating
modulate the IGM's density structure, but their impact on the IGM
recombination rate remains unclear because $C_{\mathrm{b}}$ averages over gas
that is in reality self-shielded.  In order to understand how feedback
processes modulate the abundance of gas at different densities, we
show the mass-weighted
probability density function of gas (PDF) at $z=6$ as a function
of overdensity $\Delta\equiv\rho/\langle\rho\rangle$ in the top
panel of Figure~\ref{fig:pdf}.  Comparing the blue and black
solid curves to the lower curves shows how photoheating
smooths over filaments and removes gas from low-mass haloes
(see also~\citealt{paw09a}).  At densities higher than the
adopted threshold for star formation (vertical dashed line), the
PDF is depleted by star formation and outflows.

Above $\Delta\sim240$, both reionisation simulations (short-dashed
red and dot-dashed magenta) show an increasing gas abundance because
the EUVB is attenuated in these regions owing to self-shielding.
Outflows modestly increase the gas abundance throughout the region
$10 < \Delta < 1000$.  The long-dashed green curve shows the PDF
from the r6n256nWHM01 simulation, which omits outflows and assumes
a spatially-homogeneous EUVB (with no self-shielding).  Its gas
abundance is suppressed with respect to the self-shielding
calculations in regions where $\Delta>1000$ and enhanced near
$\Delta\sim200$ because the EUVB evaporates dense gas that should
in reality be shelf-shielded.

Differences in the gas density PDFs translate into different amounts
of gas clumping.  In the middle panel of Figure~\ref{fig:pdf}, we show
the cumulative clumping factor as a function of the threshold
overdensity $C(<\Delta_{\rm{thr}})$, defined following~\citet{paw09a}:
\begin{equation}\label{eqn:dtr}
C(<\Delta_{\rm{thr}}) \equiv \frac{\int_{0}^{\Delta_{\rm{thr}}} d\Delta \Delta^2 \mathcal{P}(\Delta)}{\int_{0}^{\Delta_{\rm{thr}}} d\Delta \mathcal{P}(\Delta)}
\end{equation}
Here, $\mathcal{P}(\Delta)$ is the volume-weighted probability
distribution function of gas density.  Note that each curve's value at
the star formation threshold equals the corresponding simulation's
$C_{\mathrm{b}}$ at $z=6$ in Figure~\ref{fig:cb_feedback}.
The solid blue and black curves
confirm that outflows have little impact on the clumping factor without
an EUVB.  Adopting a spatially-homogeneous EUVB (green long-dashed)
reproduces the~\citet{paw09a} result (black dotted) as expected, with
slight differences likely owing to their addition of extra heating
around the reionisation redshift.  Our radiation-hydrodynamic
simulations again predict very similar clumping factors.  Outflows
boost clumping as expected from the top panel, but the change is slight.
In particular, we confirm the finding by~\citet{paw09a} that outflows change
$C(<\Delta_{\rm{thr}})$ by $<10\%$ when averaged over gas with $\Delta<100$.  This is
somewhat nontrivial given that our outflow prescription ejects roughly twice
as much gas per unit stellar mass formed as theirs at this redshift.
Overall, the impact of changing the EUVB and the outflow model on
$C(<\Delta_{\rm{thr}})$ is quite modest.  The fact that the gas density PDF at low
densities ($\Delta < 200$) is relatively robust to these factors
indicates that the recombination rate in the diffuse IGM---which is the
important quantity for understanding the progress of reionisation---is
only weakly sensitive to the uncertain details of star formation and
feedback.  At higher densities, the impact of feedback and self-shielding
is much stronger (top panel), but recombinations in these regions do not
count toward the IGM recombination rate as their ionisation state is
dominated by local sources.

The ability of outflows to boost the abundance of neutral gas at moderate
overdensities ($\Delta > 240$) rather than only at the outskirts of
galaxies ($\Delta > 1000$) suggests that they could function as
an additional population of Lyman limit systems residing near or
within halos.  If so, then they would have to be self-shielded against
the EUVB.  We show in the bottom panel the neutral hydrogen fraction as a
function of $\Delta$ for our fiducial simulation.  The neutral fraction
jumps from 0 to 1 at $\Delta\approx240$, indicating that outflows could
indeed be partially self-shielded.  This possibility is consistent
with direct observations of cold gas
in galactic outflows~\citep[for example,][]{mar05,wei09} as well as
recent theoretical results~\citep{van12}.

In summary, an EUVB significantly suppresses IGM clumping by
pressurizing moderately overdense regions while outflows modestly
boost the recombination rate by
returning gas from galaxies into the IGM.  Meanwhile, self-shielding
hides the most overdense gas from the EUVB.  This boosts the mass
fraction of neutral gas while suppressing the clumping factor in
ionised regions.  Outflows are at least partially self-shielded,
suggesting that they could contribute to the IGM opacity.

\section{Application I: Photon-Counting}\label{sec:photCount}
Analytical models of reionisation use clumping factors to account
for IGM inhomogeneities~\citep{mad99}, prompting the need for
theoretical insight into the nature of clumping.  We have discussed the
various effects that modulate the value of the clumping factor that we
derive from our simulations.  In order to close the loop, we now ask
whether these clumping factors do indeed constitute a good model for
the IGM recombination rate during reionisation.  We do this by
using an analytical model to explore how well several different
clumping factors fare in reproducing the reionisation history of our
r6n256wWwRT16 simulation.  We will show that clumping factors perform
reasonably well, and that they fare better if they incorporate information
from the simulated temperature and ionisation fields.

We use as our model Equation 20 of~\citet{mad99} (this can also be
derived by taking the spatial average of Equation~\ref{eqn:dnHIdt} and
normalizing by the mean hydrogen number density):
\begin{equation} \label{eqn:dnHIdtmean}
\frac{dx_{\mathrm{H\,I,V}}}{dt} = -\dot{n}_\gamma + \mathrm{recombinations}.
\end{equation}
Here, $x_{\mathrm{H\,I,V}}\equiv\langle n_{\mathrm{H\,I}}/n_{\mathrm{H}} \rangle$ is the volume-averaged
neutral hydrogen fraction while $\dot{n}_\gamma$ is the ionising
luminosity per hydrogen atom into the IGM.  We will use the value
for $\dot{n}_\gamma$ that is predicted
directly by our simulation.  The form of the recombination term varies
depending on the clumping factor definition.

There are three caveats to this widely-used formalism.  The first is that
the second term represents a \emph{mass-average} while the third term is
generally computed in a \emph{volume-averaged} way (see, for example,
Equation~\ref{eqn:Cdef}).  To see this, recall that the second term is
meant to model the rate of growth of the ionised volume
fraction~\citep{mad99}.  The rate at which an
ionised volume $V_I$ grows depends on the ratio of the ionising luminosity
$\dot{N}_\gamma$ to the gas density at the position of the ionisation front
$n_H$, $\dot{V}_I = \dot{N}_\gamma/n_H$.  In practice, however,~\citet{mad99}
substituted the volume-averaged gas density $\langle n_H \rangle$ into the
denominator.  This approximation is only accurate if the density is
homogeneous or if the sources are widely separated (as in the case of
helium reionisation;~\citealt{mcq09}), hence it is equivalent to
assuming that the neutral mass and volume fractions are equal,
$x_{\mathrm{H\,I,M}}=x_{\mathrm{H\,I,V}}$.  It is possible to use knowledge of the
distribution of gas densities and the topology of reionisation (for
example,~\citealt{mir00}) to compute the $\dot{V}_I$ from the ionising
emissivity (or, alternatively, to compute the recombination rate from
the ionised mass fraction), but these steps are omitted in analytical
calculations.  Owing to this approximation, we will refer only to the
neutral fraction $x_{\mathrm{H\,I}}$ rather than to $x_{\mathrm{H\,I,M}}$ or $x_{\mathrm{H\,I,V}}$ within the
context of our analytical model.  The second caveat is that
Equation~\ref{eqn:dnHIdtmean} omits collisional
ionisations because they are nontrivial to model
analytically.  Finally, the third term requires knowledge of the IGM's
density and temperature distributions that is generally encapsulated via
clumping factors.  All of these simplifications are relaxed in our numerical
calculations.  In this section, we will use our simulations to ask
whether they are indeed valid.
\begin{figure}
\centerline{
\setlength{\epsfxsize}{0.5\textwidth}
\centerline{\epsfbox{fig.photCount.ps}}
}
\caption{The dependence of the reionisation history on the definition of the
clumping factor, as compared to our fiducial radiation hydrodynamic simulations.
In photon-starved reionisation scenarios, a rapidly-growing clumping factor
can lead to a double-reionisation history although this does not happen with
our preferred definition of the clumping factor.
}
\label{fig:photCount}
\end{figure}

We compare the resulting reionisation histories with our numerical simulation
in Figure~\ref{fig:photCount}.  As a baseline, the short-dashed blue curve
illustrates the reionisation history if there are no recombinations ($C=0$).
Comparing it with the dot-dashed green curve indicates that recombinations
delay the completion of reionisation from $z\approx10$ to $z\approx7.5$.
The delay is expected given that this simulation requires 3.5 ionising
photons per hydrogen atom to reach a volume-averaged neutral fraction
of $x_{\mathrm{H\,I}}=0.01$.  The other colored
curves represent three different ways of treating recombinations.  The
solid black curve uses the clumping factor of all IGM gas and computes
the recombination term in Equation~\ref{eqn:dnHIdtmean} as:
$C_{\mathrm{b}} k_1(10^4\mathrm{K})x_{\mathrm{H\,II}}\langlen_{\mathrm{H}}\rangle$, where $x_{\mathrm{H\,II}} = 1 - x_{\mathrm{H\,I}}$
represents the ionised fraction.  Down to $z=10$, this
crude treatment is already a significant improvement over the $C=0$
curve.  Below $z=10$, recombinations begin to win over the declining
ionising emissivity, which is in turn driven by the assumed $f_{\mathrm{esc}}(z)$;
the result is a multimodal reionisation history.  The fact that the
simulation does not, in fact, yield a multimodal reionisation history
emphasizes the importance of allowing overdense regions to self-shield
against the ionising background.  The dotted red and long-dashed
magenta curves correspond to the definitions $C_{{\rm H\,II}}$ and $C_{{\rm H\,II,}10^4{\rm K}}$
from Section~\ref{ssec:feedbackAll}.  In these models, the
recombination term is
$C k_1(10^4\mathrm{K}) x_{\mathrm{H\,II}}^2 \langlen_{\mathrm{H}}\rangle$.  They are each a
clear improvement over using $C_{\mathrm{b}}$.  Recall that the dotted red curve
accounts for the order in which regions of different overdensities
reionise while the long-dashed magenta curve additionally accounts
for the impact of temperature fluctuations.  Including this information
accelerates reionisation by suppressing recombinations
(Figure~\ref{fig:clumpdefs}).

In detail, reionisation is more extended in the numerical model even
when the clumping factor accounts for both the temperature and the
ionisation fields (long-dashed magenta).  We attribute this to the fact
that the analytical model effectively assumes an infinite speed of light
$c$.  We have previously shown that artificially boosting $c$ causes
reionisation to complete earlier~\citep{fin09b} although the effect
was slightly weaker ($\Delta z = 0.1$) than what is seen here
($\Delta z\approx0.9$).  This delay is only a few times larger than
the light travel time across the photon mean free path at this epoch
($\sim10$ Mpc;~\citealt{fur05}), consistent with the possibility
that light-travel time effects delay the completion of reionisation.

There are two additional effect that extend the final stages of
reionisation in numerical models.  The first results from the fact
that the timescale for ionising an absorber varies inversely with its
geometrical cross section.  For example, an EUVB with amplitude $\Gamma$
at the frequency corresponding to an absorption cross-section $\sigma$
ionises a spherical absorber of uniform density $n$ and radius $r$
on a timescale $\frac{2}{3} n\sigma r/\Gamma$ compared to $1/\Gamma$
in the optically-thin case.  Approximating Lyman limit systems as
uniform spheres in hydrostatic equilibrium at a temperature of
$10^4 K$, we find that systems with overdensities of 10 and 100
are ionised by a background $\Gamma=2\times10^{-13}$s$^{-1}$ on
timescales of $\approx$ 8 and 24 Myr at $z=6$, compared to $\approx0.2$
Myr in the optically thin case.  This effect is omitted in analytical
models but present in numerical simulations, although it probably
delays the completion of reionisation by only $\Delta z\approx 0.1$.

The second delay results from the fact that photons with high energy
and small absorption cross section take longer to be absorbed than
photons near the Lyman limit once the ionised volume fraction is large.
This effect is missing from our simulations owing
to our monochromatic radiation transport solver.  However, it is not
large: for $z\geq6$, energies less than 54.4 eV, and neutral hydrogen
fractions greater than 0.01, the delay is less than 20 Myr.

We may quantify the error in the reconstructed reionisation history using
the optical depth to Thomson scattering $\tau_{\mathrm{es}}$.  We compute this for the
simulations and the analytical models assuming that helium is
singly-ionised with the same ionisation fraction as hydrogen down to
$z=3$~\citep{cia12,fri12}, and doubly-ionised thereafter.  In our
simulation, the predicted value is 0.071.  Ignoring
recombinations yields $\tau_{\mathrm{es}}=0.095$, which confirms that
galaxies could have dominated reionisation modulo uncertainties regarding
the star formation efficiency of low-mass halos and the IGM recombination
rate.  Adopting realistic clumping factors measured directly from our
simulation yields values for $\tau_{\mathrm{es}}$ that are between 0.062 (solid black)
and 0.082 (long-dashed magenta).  The discrepancy between the long-dashed
magenta and dot-dashed green curves can be regarded as an estimate of the
uncertainty in analytical reionisation calculations owing to complications
such as light-travel time effects, shadowing, and source clustering;
it is $\approx10\%$.  In short, analytical models of reionisation
can reproduce the behavior of more realistic models reasonably well, with
improvement possible if the adopted recombination rate accounts
for the inhomogeneous temperature and ionisation fields.  Our analytical
modeling also confirms that collisional ionisations and the other
complications that we have mentioned are subdominant.

We note that~\citet{gne00} previously performed an extensive comparison
between his numerical radiation hydrodynamic simulations and analytical models,
finding broad support for the basic assumptions underlying analytical models.  Our
results agree with his.  The principal difference is that, in computing
the volume-averaged recombination rate, their analysis did not distinguish
between gas within galaxies and the IGM. This led to a much higher value
for the clumping factor (100--200) than is appropriate for densities
that are conventionally associated with the IGM (for example,
$\Delta < 1000$).  Our simulations sidestep this complexity by using
$\Delta_{\rm{lls}}$ as a physically-motivated model for the boundary between the IGM
and the condensed gas whose absorptions are more appropriately
attributed to $f_{\mathrm{esc}}$.

\section{Application II: The Photon Cost of Reionising Filaments}\label{sec:photCost}
\begin{figure}
\centerline{
\setlength{\epsfxsize}{0.5\textwidth}
\centerline{\epsfbox{fig.photCost.ps}}
}
\caption{(bottom) The dependence of the overlap redshift on heating redshift:
Delaying photoheating from $z=19.5$ to $z=7.5$ delays overlap by $\Delta z=3$.
(top) The reionisation photon budget per hydrogen atom as a function of heating redshift:
Delaying photoheating can more than double the photon cost of reionisation.
The crosses correspond to a homogeneous IGM ($C=1$).
}
\label{fig:photCost}
\end{figure}
While much has been learned about the photon cost of removing
gas from minihalos~\citep{hai01,bar02,sha04,ili05,cia06}, the ability
of moderately-overdense gas to delay reionisation has received
less attention.  This is potentially important for two reasons.  First,
analytical models may underestimate the reionisation photon budget by deriving
the clumping factor from hydrodynamic simulations in which the IGM is
effectively pre-heated.  For example, it is possible to compute the clumping
factor at high redshift by applying Equation~\ref{eqn:dtr} to the baryonic
density PDF of~\citet{mir00}~\citep{wyi08,pri10}.  However, that PDF
derives from simulations in which the diffuse IGM is already
smoothed by an EUVB.  In reality, there is a photon cost associated
with smoothing the IGM because the initially cold IGM has a higher
clumping factor (Figure~\ref{fig:pdf}).  Ignoring this cost may
cause models to underestimate the photon budget of reionisation.
Alternatively, it is possible that the IGM is heated before
reionisation is well under way.  For example, an early X-ray
background could suppress gas clumping prior to the onset of
reionisation~\citep{oh03}.  In this case, filaments would absorb fewer
ionising photons than in models (such as ours) in which the primary
heating source is UV photons.  In this section, we estimate the
magnitude of this uncertainty by asking how many ionising photons are
required to reionise filaments, or regions with overdensities of
$\Delta=$1--50.

Our approach involves decoupling reionisation from heating and
determining how reionisation changes if the universe is heated at
different redshifts.  If the universe is heated at high redshift
(such as $z=19.5$), then the photon cost of reionisation is smaller
because the universe remains smooth after being heated.  By contrast,
if the universe is heated at a lower redshift (such as $z=7.5$), then
reionisation is more expensive because more of its gas has condensed into
filaments by the time that reionisation is under way.  By comparing these
histories, we gain some insight (albeit not completely self-consistently) into
the photon price that the universe pays for waiting to ionise its filaments;
this is probably comparable to the actual photon cost of reionising
moderately-overdense regions.

\citet{paw09a} provide fitting functions for the baryonic clumping factor
in simulations where an instantaneously imposed spatially-homogeneous
EUVB is used to heat the IGM at a variety of redshifts.  We have
combined their fitting functions (choosing their $C_{50}$) with the
predicted ionising emissivity from our r6n256wWwRT16 simulation
within our analytical model and evaluated how reionisation
depends on the smoothing redshift.  We show in the bottom
panel of Figure~\ref{fig:photCost} how the reionisation redshift (where
$x_{\mathrm{H\,I}}\rightarrow 0.01$) varies.  Delaying smoothing from $z=19.5$
to $z=7.5$ delays reionisation by $\Delta z \approx 4.5$.  This relatively
strong dependence owes to the extended reionisation history that is
enforced by our choice of $f_{\mathrm{esc}}(z)$.  We have verified that extracting
our emissivity history from simulations that assume a constant
$f_{\mathrm{esc}}=0.5$ yields a smaller delay of $\Delta z = 1.1$ (not shown)
because the growth rate of collapsed matter (such as stars) is much
faster than the growth rate of moderately overdense structures
(such as filaments).

The top panel maps the reionisation redshift into the number of ionising
photons consumed per hydrogen atom.  If the IGM is
heated at $z=19.5$, reionisation consumes 1.8 photons per
hydrogen.  Observations suggest that reionisation was well under way
by $z=$9--10~\citep{pri10,kuh12}.  Adopting this as the photoheating
redshift yields a total photon cost of 2.5--2.9 photons per hydrogen.
Subtracting, we find that the universe pays a price of 0.7--1.1 photons
per hydrogen for waiting to heat its filaments.  Put
differently, analytical models that assume that the universe is smoothed
at a high redshift underestimate the total photon cost of reionisation
by 0.7--1.1 photons per hydrogen atom.

These calculations may underestimate the total photon cost of reionisation
because they do not account for gas that is condensed into minihalos,
or halos with virial temperatures below $10^4$ K ($\approx10^8M_{\odot}$ for
redshifts of 6--10).  Prior to reionisation, much of the gas may have
been condensed into halos with virial masses between $10^6$ and $10^8M_{\odot}$.
The clumping factors of~\citet{paw09a} were derived from simulations that
resolved halos at the upper end of this mass range with 100 particles,
hence lower-mass systems are effectively missing.  These halos can
consume up to 5 photons per hydrogen atom~\citep{sha04}, with the
implication that the amount of gas clumping could have been much
higher if the IGM temperature did not exceed 1000 K prior to reionisation.

On the other hand, an early X-ray background could deposit an early
entropy floor due to the large mean free path of
X-rays~\citep{oh01,ven01,ric04,mad04}.  Such a background is widely
predicted in many models of reionisation and could arise from X-ray
binaries, inverse Compton scattering in supernova remnants, or early
mini-quasars.  This process significantly reduces gas
clumping~\citep{oh03}, and if it occurs well before the reionisation
epoch then it could decrease the reionisation photon budget.

In short, while filaments may have absorbed $\sim1$ additional photon
per hydrogen, the total number of photons required to complete reionisation
remains uncertain owing to the unknown abundance of small absorbers.
Future work incorporating X-ray heating and a wider dynamic range will
be required in order to explore these processes.

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

Could galaxies have produced enough ionising photons to reionise the
Universe? Could they have maintained the IGM's ionisation state at $z=6$?
Two major uncertainties hamper efforts to connect observations of
galaxies with measurements of the neutral hydrogen fraction in order
to address these questions.  The first is the fraction of ionising
photons that escape galaxies' ISMs, $f_{\mathrm{esc}}$, and the second is the
recombination rate of the IGM, which is often parametrized by a
clumping factor $C$.  The observed population of galaxies can
maintain an ionised IGM if the ratio $C/f_{\mathrm{esc}}$ lies within the range
1--10~\citep[][Figure 15]{fin10}.

We previously showed
that our simulations reproduce the observed UV luminosity function of
galaxies and self-consistently complete reionisation by $z=6$ when we
assume $f_{\mathrm{esc}}=0.5$~\citep{fin11}.  While encouraging, those simulations
suffered from two drawbacks.  First, they overproduced the observed
EUVB amplitude at $z\leq6$ and underproduced $\tau_{\mathrm{es}}$.  Second, the
radiation transport solver's limited spatial resolution did not shield
moderately-overdense regions from the EUVB, with the result that the
predicted clumping factor of ionised gas was artificially high
($\approx12$ at $z=6$).  Consequently, completing reionisation
($x_{\mathrm{H\,I,V}}<0.01$) required up to 5 ionising photons per hydrogen atom.

In this work, we have improved on those calculations by adopting a
subgrid self-shielding treatment and a redshift-dependent $f_{\mathrm{esc}}$.
Self-shielding limits the clumping factor and the reionisation photon
budget by isolating overdense gas from the EUVB.  Meanwhile, an
evolving $f_{\mathrm{esc}}$ (1) reconciles the model with the observed $\tau_{\mathrm{es}}$ by
increasing the electron fraction at early times; (2) heats
the IGM at an earlier redshift, which suppresses gas clumping at
later times; and (3) raises the Gunn-Peterson optical depth at $z=5$
into improved agreement with observations.  Together, we find that
these
modifications suppress $C$ to values in the range 2.7--3.3 at $z=6$
(the uncertainty comes from possibly overestimating
the overdensity at which gas fully self-shields; see below).
For our adopted $f_{\mathrm{esc}}(z)$, this leads to $C/f_{\mathrm{esc}}\approx$25--30
at $z=6$, slightly too high for currently-observed galaxy
populations to maintain a completely ionised IGM.  Our
simulations overcome this barrier by resolving galaxies down
to absolute magnitudes of $M_{\mathrm{UV}}=-15$, which is 2--3
magnitudes fainter than current limits for observing galaxies
directly but consistent with requirements inferred from other
star formation tracers~\citep{tre12}.

Despite these improvements, there are a number of ways in which we
can build on this work in the future.  First, we use the global EUVB
to compute $\Delta_{\rm{lls}}$.  Near sources, this underestimates the EUVB,
leading to an overdensity threshold that is too low.  In voids, it
overestimates the EUVB, leading to too little
self-shielding~\citep[see also][]{cro11}.  It would be preferable
to estimate the EUVB by averaging
over regions that are a few times larger than the expected length scale of
self-shielding systems in order to incorporate the EUVB's small-scale
fluctuations into the ionisation field.  This improvement would result in
more-ionised overdensities and more-neutral voids (or in other words a
more strongly inside-out reionisation topology), but the overall impact
on $C$ is difficult to predict.

Second, our results would benefit from an improved understanding
of the IGM's temperature $T$.  The IGM recombination rate is
sensitive to $T$ through the recombination coefficient $k_2$ and
through the minimum overdensity $\Delta$ at which gas is neutral, which
in turn is proportional to the self-shielding threshold $\Delta_{\rm{lls}}$.
For temperatures between
$10^{3.5}$--$10^{5.5}$, $k_2\propto T^{-0.77}$.  Meanwhile,
$C(<\Delta_{\rm{thr}}) \propto \Delta^{0.25}$ (Figure~\ref{fig:pdf}) and
$\Delta_{\rm{lls}}\propto T^{0.3}$ (Equation~\ref{eqn:dlls_typ}), leading
to a weak overall dependence
of $C(<\Delta_{\rm{thr}})\propto T^{0.075}$.  Combining these, we find
that the overall IGM recombination rate varies as $T^{-0.7}$.

Uncertainty in $T$ stems from the unknown latent heat of photoionisation
and from resolution limitations.  We have effectively assumed
that reionisation heats gas to $\approx15,000$ K.  While this assumption
yields a mean IGM temperature that is consistent with observations
(Section~\ref{sec:sims}), the observations are themselves uncertain
owing to the necessity to correct the measured IGM temperature for
the influence of nearby quasars~\citep{bol12}.  Assuming a (harder,
softer) ionising background would lead to a (hotter, colder)
post-reionisation IGM.  Additionally, the limited spatial
resolution of our radiation solver may cause our
simulations to underestimate the post-reionisation gas temperature
because ionisation fronts are artificially broadened, increasing the
amount of time during which partially-ionised gas can cool through
collisional excitation of neutral hydrogen~\citep{mir94}.  The magnitude
of this effect is likely small compared to the uncertain latent heat of
photoionisation.  Nevertheless, improved measurements of the IGM
temperature at $z>5$ would constrain our model further.

Third, our simulations expose somewhat too much overdense gas to the EUVB,
boosting the IGM recombination rate.  For example, the
threshold overdensity for self-shielding at $10^4$K and $z=6$ should
be $\Delta_{\rm{lls}}=25$ (Equation~\ref{eqn:dlls_typ}).  If the neutral column
density of self-shielded systems varies strongly with $\Delta$, then
systems with overdensity more than a few times greater than $\Delta_{\rm{lls}}$ should
be essentially neutral.  For example, the high-resolution calculations
of~\citet{mcq11} indicate that the neutral fraction should increase
rapidly above $\Delta\sim100$ (their Figure 3).  By contrast, the gas
in our simulations remains ionised until $\Delta\approx200$
(Figure~\ref{fig:pdf}), which may be slightly too high.  If so,
then it owes largely to the fact that our simulations assume that
all photons have an energy of 1.3 times the threshold energy for
photoionisation $h\nu_{\rm{LL}}$ in order to photoheat the IGM.  As we
argued in the discussion of Figure~\ref{fig:dlls}, this assumption may
artificially boost $\Delta_{\rm{lls}}$ by a factor of 1.7.  Invoking once
again the scaling $C\propto\Delta_{\rm{lls}}^{0.25}$ (Figure~\ref{fig:pdf}), we
find that,
if the amplitude of the ionising background were unchanged but the
radiation field were dominated by photons with energy nearer to
$h\nu_{\rm{LL}}$, then the clumping factor of ionised gas would be
lower by $\approx12\%$.  Combining this uncertainty with the range of
values inferred from our fiducial and high-resolution simulations
(r6wWwRT16 and r6wWwRT32), we conclude that the clumping factor at
$z=6$ lies between 2.7--3.3.

Finally, our limited cosmological volume forces us to adopt a rather
dramatic evolution for $f_{\mathrm{esc}}(z)$.  We do this because our predicted $\tau_{\mathrm{es}}$
and $\tau_\alpha$ are not converged with respect to simulation volume or mass
resolution, as is well known for simulation volumes of this
size~\citep{gne06,bol09}.  This is not to say that the need for a
decreasing $f_{\mathrm{esc}}(z)$ is purely a resolution limitation; empirical
arguments strongly support such a trend and in fact our adopted
relation lies within the observationally-inferred range~\citep{kuh12}.
However, with a larger volume, structure formation
would begin sooner, allowing us to
adopt a weaker redshift-dependence for $f_{\mathrm{esc}}(z)$ without compromising
the agreement between the observed and predicted $\tau_{\mathrm{es}}$ and $\tau_\alpha$.
This would change the predicted IGM temperature-density relationship
as well as the overall gas density distribution.  Both of these changes
could modify the predicted IGM recombination rate, although the sign
of the effect is difficult to predict.

\section{Summary} \label{sec:sum}
We have used a suite of cosmological radiation hydrodynamic simulations
to explore the impact of outflows and an EUVB on the recombination rate
of the reionisation-epoch IGM.  Feedback processes modify
the IGM's gas density, ionisation, and temperature fields.  We may
illustrate their impact on the gas density distribution by considering
the clumping factor
$C_{100}\equiv\Delta^2 \mathcal{P}(\Delta)$ at redshift $z=6$ and
overdensity $\Delta=100$; the overall clumping factor $C_{\mathrm{b}}$ (which neglects
the ionisation and temperature fields) scales with this value.
\begin{itemize}
\item In the absence of feedback, baryons follow dark matter,
and $C_{100}=10.7$ as previously found in N-body simulations~\citep{ili07,rai11}.
\item Galactic outflows return star-forming gas into the IGM, boosting
$C_{100}$.  The increase is less than 30\% for any EUVB amplitude
(including zero).
\item An ionising background pressurizes the diffuse IGM, suppressing
$C_{100}$ to 2.0.
\end{itemize}
In reality, the IGM's volume-averaged recombination rate depends on
its density, ionisation, and temperature fields.  Taking this information
into account one piece at a time within the context of our fiducial
simulation, we find (still at $z=6$):
\begin{itemize}
\item The clumping factor of all baryons below our adopted threshold
density for star formation is $C_{\mathrm{b}}=16.4$.
\item Averaging only over the ionised gas yields a much lower value
of $C_{{\rm H\,II}}=4.9$ because overdense gas is self-shielded.
\item Accounting for the temperature field further suppresses the
recombination rate to $C_{{\rm H\,II,}10^4{\rm K}}=3.3$ because gas in the recent aftermath
of reionisation is somewhat hotter than the canonical $10^4$ K.
\item The temperature-corrected clumping factor of ionised regions
increases from $0.79$ at $z=15$ to 3.3 at $z=6$.  Equation~\ref{eqn:cfit}
reproduces this evolution reasonably well throughout $z=15\rightarrow5$.
\item Our fiducial simulation may overestimate the minimum density of
neutral gas.  Correcting for
this may suppress $C_{{\rm H\,II,}10^4{\rm K}}$ by an additional factor of 12\%.  Hence
our most realistic estimate of the ionisation- and temperature-corrected
clumping factor at $z=6$ is that it lies within the range 2.7--3.3.
\end{itemize}

We have constructed an analytical reionisation model and tested how well
different definitions of the clumping factor reproduce our numerical
simulation's reionisation history.  We find that the clumping factor
averaged over all IGM baryons $C_{\mathrm{b}}$ significantly overestimates the
IGM recombination rate because it treats self-shielded gas as if it
were optically thin.  Accounting for the ionisation field through
$C_{{\rm H\,II}}$ improves agreement with the numerical simulation, although
the recombination rate is still overestimated because the
post-reionisation IGM is in reality slightly hotter than the canonical
$10^4$ K.  By contrast, using $C_{{\rm H\,II,}10^4{\rm K}}$ to
account additionally for the IGM temperature field yields a monotonic
reionisation history whose $\tau_{\mathrm{es}}$ overestimates the numerical result by
only 10\%.  We speculate that this small remaining difference owes at
least partly to light-travel time effects that are missing from
the analytical calculation although other effects may also contribute.

We have used our analytical model to estimate the photon
cost of reionising filaments, or regions with $\Delta=$1--50.  With
our simulated emissivity history,
if some process (other than reionisation) heats the IGM at $z=19.5$
then the photon cost of reionisation is $\gamma/H=1.8$ photons per
hydrogen atom.  Delaying the redshift at which the IGM is heated
from $z=19.5$ to $z=9.0$ increases $\gamma/H$ to 2.9 because the
clumping factor increases.  The difference between these numbers
indicates that the Universe pays a price of $\approx1$ photon
per hydrogen in order to reionise its filaments.

The ionisation state of dense gas has observational implications
that merit further study.  Our simulations show that, in the presence
of self-shielding, outflows can boost the amount of gas at densities
corresponding to the circumgalactic medium ($10 < \Delta < 1000$).
If a significant fraction of this gas remains neutral, then it could
play two important roles in modulating the growth of the EUVB.  Within
halos, it could constitute a significant absorbing column, modifying
the fraction of ionising photons that escape the host halo.  This effect
could mimic a redshift-dependent $f_{\mathrm{esc}}$, which seems to be required by
observations~\citep{ino06,fin11,haa12,kuh12}.  Additionally, gas that
travels past the virial radius while retaining a self-shielded component
could boost the abundance of Lyman limit systems, as suggested
by~\citet{van12}.

A complementary motivation for improving our understanding of
self-shielded gas is that it can be observed directly in the form of
damped Lyman $\alpha$ systems and low-ionisation metal absorbers.
The IGM's ionisation field is a major theoretical uncertainty hampering
efforts to interpret observations of low-ionisation metal ions.  For
example,~\citet{opp09} studied low-ionisation metal absorbers
through the use of homogeneous ionising backgrounds.  While this
approach may be appropriate for ions such as SiIV and CIV, it is
less appropriate for low-ionisation species such as OI and SiII
because it artificially ionises, heats, and smooths gas that should
in reality be self-shielded.  Our current simulations allow dense gas to
self-shield, hence they will yield improved predictions for the abundance
of low-ionisation metal absorbers.  This will reduce the theoretical
parameter space and increase the power of current observations.

\section*{Acknowledgements}
We thank the referee for a thoughtful and detailed report that
improved the draft.  We thank M.\ McQuinn and M.\ Prescott for helpful
conversations.  KF thanks D.\ Godinez and J.\ Sivapalan for exploring
the ability of galactic outflows to self-shield within our simulations.
Our simulations were run on the University of Arizona's Xeon cluster.
Support for this work was provided by the National Science Foundation
through grant AST-0907998, the NASA Astrophysics Theory
Program through grant NNG06GH98G, and grant number
HST-AR-10647 from the SPACE TELESCOPE SCIENCE INSTITUTE, which is
operated by AURA, Inc. under NASA contract NAS5-26555.  Support for
this work, part of the Spitzer Space Telescope Theoretical Research
Program, was also provided by NASA through a contract issued by the
Jet Propulsion Laboratory, California Institute of Technology under
a contract with NASA.  KF gratefully acknowledges support from NASA
through Hubble Fellowship grant HF-51254.01 awarded by the Space
Telescope Science Institute, which is operated by the Association
of Universities for Research in Astronomy, Inc., for NASA, under
contract NAS 5-26555.  SPO acknowledges support from NSF grant
AST-0908480.