We present the methodology of a photon-conserving, spatially-adaptive, ray-tracing radiative transfer algorithm, designed to run on multiple parallel Graphic Processing Units (GPUs). Each GPU has thousands computing cores, making them ideally suited to the task of tracing independent rays. This ray-tracing implementation has speed competitive with approximate momentum methods, even with thousands of ionization sources, without sacrificing accuracy and resolution. Here, we validate our implementation with the selection of tests presented in the "cosmological radiative transfer codes comparison project," to demonstrate the correct behavior of the code. We also present a selection of benchmarks to demonstrate the performance and computational scaling of the code. As expected, our method scales linearly with the number of sources and with the square of the dimension of the 3D computational grid. Our current implementation is scalable to an arbitrary number of nodes possessing GPUs, but is limited to a uniform resolution 3D grid. Cosmological simulations of reionization with tens of thousands of radiation sources and intergalactic volumes sampled with 1024$^3$ grid points take about 30 days on 64 GPUs to reach complete reionization.

We present semi-analytic models of the epoch of reionization focusing on the differences between continuous and bursty star formation (SF). Our model utilizes physically motivated analytic fits to 1D radiative transfer simulations of HII regions around dark matter halos in a representative cosmic volume. Constraining our simulations with observed and extrapolated UV luminosity functions of high redshift galaxies, we find that for a fixed halo mass, stellar populations forming in bursty models produce larger HII regions which leave behind long-lived relic HII regions which are able to maintain partial ionization in the intergalactic medium (IGM) in a manner similar to an early X-ray background. The overall effect is a significant increase in the optical depth of the IGM, $\tau_e$, and a milder increase of the redshift of reionization. To produce $\tau_e =0.066$ observed by Planck and complete reionization by redshift $z_{\rm re}\sim 6$, models with bursty SF require an escape fraction $f_{\rm esc}\sim 2\%-10\%$ that is 2-10 times lower than $f_{\rm esc}\sim 17\%$ found assuming continuous SF and is consistent with upper limits on $f_{\rm esc}$ from observations at $z=0$ and $z\sim 1.3-6$. The ionizing photon budget needed to reproduce the observed $\tau_e$ and $z_{\rm re}$ depends on the period and duty cycle of the bursts of SF and the temperature of the neutral IGM. These results suggest that the tension between observed and predicted ionizing photon budget for reionization can be alleviated if reionization is driven by short bursts of SF, perhaps relating to the formation of Population~III stars and compact star clusters such as proto-globular clusters.