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 multilevel 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 (I) 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 calculated directly at each redshift step, and (3) correct treatment of the He I atom, including the triplet and singlet states. We find that x(e)(=n(e)/n(H)) is approximately 10% smaller at redshifts less than or similar to 800 than in previous calculations, as a result of the nonequilibrium of the excited states of H that 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 the ionization fraction. 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.