  Primordial magnetic fields (PMFs) may explain observations of magnetic fields
on extragalactic scales. They are most cleanly constrained by measurements of
cosmic microwave background radiation (CMB) anisotropies. Their effects on
cosmic recombination may even be at the heart of the resolution of the Hubble
tension. We present the most detailed analysis of the effects of PMFs on cosmic
recombination to date. To this end we extend the public magneto-hydrodynamic
code {\sl ENZO} with a new cosmic recombination routine, Monte-Carlo
simulations of Lyman-$\alpha$ photon transport, and a Compton drag term in the
baryon momentum equation. The resulting code allows us, for the first time, to
realistically predict the impact of PMFs on the cosmic ionization history and
the clumping of baryons during cosmic recombination. Our results identify the
importance of mixing of Lyman-$\alpha$ photons between overdense- and
underdense- regions for small PMF strength. This mixing speeds up recombination
beyond the speed-up due to clumping. For non-helical PMFs with a Batchelor
spectrum we find a surprising dependency of results on ultra-violet magnetic
modes. The present study shall serve as a theoretical foundation for a future
precise comparison of recombination with PMFs to CMB data.
