Rabacus: A Python Package for Analytic Cosmological Radiative Transfer Calculations

We describe Rabacus, a Python package for calculating the transfer of hydrogen ionizing radiation in simplified geometries relevant to astronomy and cosmology. We present example solutions for three specific cases: 1) a semi-infinite slab gas distribution in a homogeneous isotropic background, 2) a spherically symmetric gas distribution with a point source at the center, and 3) a spherically symmetric gas distribution in a homogeneous isotropic background. All problems can accommodate arbitrary spectra and density profiles as input. The solutions include a treatment of both hydrogen and helium, a self-consistent calculation of equilibrium temperatures, and the transfer of recombination radiation. The core routines are written in Fortran 90 and then wrapped in Python leading to execution speeds thousands of times faster than equivalent routines written in pure Python. In addition, all variables have associated units for ease of analysis. The software is part of the Python Package Index and the source code is available on Bitbucket at https://bitbucket.org/galtay/rabacus . In addition, installation instructions and a detailed users guide are available at http://pythonhosted.org//rabacus .

The Impact of Different Physical Processes on the Statistics of Lyman-limit and Damped Lyman-alpha Absorbers

We compute the z = 3 neutral hydrogen column density distribution function f(NHI) for 19 simulations drawn from the OWLS project using a post-processing correction for self-shielding calculated with full radiative transfer of the ionising background radiation. We investigate how different physical processes and parameters affect the abundance of Lyman-limit systems (LLSs) and damped Lyman-alpha absorbers (DLAs) including: i) metal-line cooling; ii) the efficiency of feedback from SNe and AGN; iii) the effective equation of state for the ISM; iv) cosmological parameters; v) the assumed star formation law and; vi) the timing of hydrogen reionization . We find that the normalisation and slope, D = d log10 f /d log10 NHI, of f(NHI) in the LLS regime are robust to changes in these physical processes. Among physically plausible models, f(NHI) varies by less than 0.2 dex and D varies by less than 0.18 for LLSs. This is primarily due to the fact that these uncertain physical processes mostly affect star-forming gas which contributes less than 10% to f(NHI) in the the LLS column density range. At higher column densities, variations in f(NHI) become larger (approximately 0.5 dex at NHI = 10^22 cm^-2 and 1.0 dex at NHI = 10^23 cm^-2) and molecular hydrogen formation also becomes important. Many of these changes can be explained in the context of self-regulated star formation in which the amount of star forming gas in a galaxy will adjust such that outflows driven by feedback balance inflows due to accretion. Data and code to reproduce all figures can be found at the following url: https://bitbucket.org/galtay/hi-cddf-owls-1

Urchin: A Reverse Ray Tracer for Astrophysical Applications

We describe URCHIN, a reverse ray tracing radiative transfer scheme optimised to model self-shielding from the post-reionisation ultraviolet (UV) background in cosmological simulations. The reverse ray tracing strategy provides several benefits over forward ray tracing codes including: (1) the preservation of adaptive density field resolution (2) completely uniform sampling of gas elements by rays; (3) the preservation of galilean invariance; (4) the ability to sample the UV background spectrum with hundreds of frequency bins; and (5) exact preservation of the input UV background spectrum and amplitude in optically thin gas. The implementation described here focuses on Smoothed Particle Hydrodynamics (SPH). However, the method can be applied to any density field representation in which resolution elements admit ray intersection tests and can be associated with optical depths. We characterise the errors in our implementation in stages beginning with comparison to known analytic solutions and ending with a realistic model of the z = 3 cosmological UV background incident onto a suite of spherically symmetric models of gaseous galactic halos.

Cold accretion flows and the nature of high column density HI absorption at redshift 3

Simulations predict that galaxies grow primarily through the accretion of gas that has not gone through an accretion shock near the virial radius and that this cold gas flows towards the central galaxy along dense filaments and streams. There is, however, little observational evidence for the existence of these cold flows. We use a large, cosmological, hydrodynamical simulation that has been post-processed with radiative transfer to study the contribution of cold flows to the observed z=3 column density distribution of neutral hydrogen, which our simulation reproduces. We find that nearly all of the HI absorption arises in gas that has remained colder than 10^5.5 K, at least while it was extragalactic. In addition, the majority of the HI is rapidly falling towards a nearby galaxy, with non-negligible contributions from outflowing and static gas. Above a column density of N_HI = 10^17 cm^-2, most of the absorbers reside inside haloes, but the interstellar medium only dominates for N_HI > 10^21 cm^-2. Haloes with total mass below 10^10 Msun dominate the absorption for 10^17<N_HI < 10^21 cm^-2, but the average halo mass increases sharply for higher column densities. Although very little of the HI in absorbers with N_HI <~ 10^20 cm^-2 resides inside galaxies, systems with N_HI > 10^17 cm^-2 are closely related to star formation: most of their HI either will become part of the interstellar medium before z=2 or has been ejected from a galaxy at z>3. Cold accretion flows are critical for the success of our simulation in reproducing the observed rate of incidence of damped Lyman-alpha and particularly that of Lyman limit systems. We therefore conclude that cold accretion flows exist and have already been detected in the form of high column density HI absorbers.

Through Thick and Thin - HI Absorption in Cosmological Simulations

We investigate the column density distribution function of neutral hydrogen at redshift z = 3 using a cosmological simulation of galaxy formation from the OverWhelmingly Large Simulations (OWLS) project. The base simulation includes gravity, hydrodynamics, star formation, supernovae feedback, stellar winds, chemodynamics, and element-by-element cooling in the presence of a uniform UV background. Self-shielding and formation of molecular hydrogen are treated in post-processing, without introducing any free parameters, using an accurate reverse ray-tracing algorithm and an empirical relation between gas pressure and molecular mass fraction. The simulation reproduces the observed z = 3 abundance of Ly-A forest, Lyman Limit and Damped Ly-A HI absorption systems probed by quasar sight lines over ten orders of magnitude in column density. Self-shielding flattens the column density distribution for NHI > 10^18 cm-2, while the conversion to fully neutral gas and conversion of HI to H2 steepen it around column densities of NHI = 10^20.3 cm-2 and NHI = 10^21.5 cm-2, respectively.

Hydrogen Self Shielding in a Quasar Proximity Zone

We calculate the distribution of HI within 750 proper kpc/h of a quasar, Lbol = 1.62e13 Lsun, powered by an SMBH, Mbh = 4.47e8 Msun, at z = 3. Our numerical model includes a cosmological hydrodynamic simulation that tracks the self consistent growth and thermal feedback of black holes calculated using GADGET-3 as well as a detailed post-processing ray tracing treatment of the non-uniform ionizing radiation field calculated using SPHRAY, which naturally accounts for the self shielding of optically thick systems. We show that the correct treatment of self shielding introduces a flattening feature into the neutral column density distribution around Log NHI = 20 and that regions with the lowest neutral fractions are not those with the highest density gas. For comparison, we solve a Ricatti equation which determines the equilibrium Hydrogen ionization fractions in the presence of a radiation field that falls off as 1/r^2 with regions above a given gas density threshold completely shielded from ionizing radiation. We demonstrate that these semi analytic models cannot reproduce the HI field calculated using SPHRAY. We conclude by comparing our models of this single proximity zone to observations by Hennawi and Prochaska of the absorption spectra of background quasars which are coincident on the sky with foreground quasars in their Quasars Probing Quasars (QPQ) series of papers. Compared to the QPQ sample, we find a factor of 3 fewer optically thick (Log NHI > 17.2) systems around our quasar, however the dark matter halo that hosts our simulated quasar, Mhalo = 5.25e12 Msun, is less massive than the typical QPQ host halo by a factor of four. Allowing for a linear scaling between halo mass, baryonic overdensity and number of absorbers, we estimate the typical host halo mass in the QPQ sample as 1.92e13 Msun.

SPHRAY: A Smoothed Particle Hydrodynamics Ray Tracer for Radiative Transfer

We introduce SPHRAY, a Smoothed Particle Hydrodynamics (SPH) ray tracer designed to solve the 3D, time dependent, radiative transfer (RT) equations for arbitrary density fields. The SPH nature of SPHRAY makes the incorporation of separate hydrodynamics and gravity solvers very natural. SPHRAY relies on a Monte Carlo (MC) ray tracing scheme that does not interpolate the SPH particles onto a grid but instead integrates directly through the SPH kernels. Given initial conditions and a description of the sources of ionizing radiation, the code will calculate the non-equilibrium ionization state (HI, HII, HeI, HeII, HeIII, e) and temperature (internal energy/entropy) of each SPH particle. The sources of radiation can include point like objects, diffuse recombination radiation, and a background field from outside the computational volume. The MC ray tracing implementation allows for the quick introduction of new physics and is parallelization friendly. A quick Axis Aligned Bounding Box (AABB) test taken from computer graphics applications allows for the acceleration of the raytracing component. We present the algorithms used in SPHRAY and verify the code by performing all the test problems detailed in the recent Radiative Transfer Comparison Project of Iliev et. al. The Fortran 90 source code for SPHRAY and example SPH density fields are made available on a companion website (www.sphray.org).

Radiation-induced large-scale structure during the reionization epoch: the autocorrelation function

The structures produced during the epoch of reionization by the action of radiation on neutral hydrogen are in principle different to those that arise through gravitational growth of initially small perturbations. We explore the difference between the two mechanisms using high resolution cosmological radiative transfer. Our computations use a Monte Carlo code which raytraces directly through SPH kernels without a grid, preserving the high spatial resolution of the underlying hydrodynamic simulation. Because the properties of the first sources of radiation are uncertain, we simulate a range of models with different source properties and recombination physics. We examine the morphology of the neutral hydrogren distribution and the reionization history in these models. We find that at fixed mean neutral fraction, structures are visually most affected by the existence of a lower limit in source luminosity, then by galaxy mass to light ratio, and are minimally affected by changes in the recombination rate and amplitude of mass fluctuations. We concentrate on the autocorrelation function of the neutral hydrogen, xi_HI(r) as a basic quantitive measure of Radiation Induced Structure (RIS.) All the models we test exhibit a characteristic behaviour, with xi_HI(r) becoming initially linearly antibiased with respect to the matter correlation function, reaching a minimum bias factor b~0.5 when the universe is ~10-20% ionized. After this xi_HI(r) increases rapidly in amplitude, overtaking the matter correlation function. It keeps a power law shape, but flattens considerably, reaching an asymptotic logarithmic slope of gamma ~-0.5. The growth rate of HI fluctuations is exponentially more rapid than gravitational growth over a brief interval of redshift Deltaz ~ 2-3.