How exactly did the Universe become neutral?

We present a refined treatment of H, He I, and He II recombination in the early Universe. The difference from previous calculations is that we use multi-level atoms and evolve the population of each level with redshift by including all bound-bound and bound-free transitions. In this framework we follow several hundred atomic energy levels for H, He I, and He II combined. The main improvements of this method over previous recombination calculations are: (1) allowing excited atomic level populations to depart from an equilibrium distribution; (2) replacing the total recombination coefficient with recombination to and photoionization from each level directly at each redshift step; and (3) correct treatment of the He I atom, including the triplet and singlet states. We find that the ionization fraction x_e = n_e/n_H is approximately 10% smaller at redshifts <~800 than in previous calculations, due to the non-equilibrium of the excited states of H, which is caused by the strong but cool radiation field at those redshifts. In addition we find that He I recombination is delayed compared with previous calculations, and occurs only just before H recombination. These changes in turn can affect the predicted power spectrum of microwave anisotropies at the few percent level. Other improvements such as including molecular and ionic species of H, including complete heating and cooling terms for the evolution of the matter temperature, including collisional rates, and including feedback of the secondary spectral distortions on the radiation field, produce negligible change to x_e. The lower x_e at low z found in this work affects the abundances of H molecular and ionic species by 10-25%. However this difference is probably not larger than other uncertainties in the reaction rates.