Interesting! Actually it's not the dt, TimeUnits bit, since that happens at the end of the linear solve anyway. So if I skip the linear solve, I still need to scale the units back to their normalized values.
Instead, it looks like in this test problem, the first few time steps meet the criteria for skipping the call to the HYPRE solver. Eventually, the species get enough out of equilibrium so that the solves kick in.
So in the 'gold standard', on those first few steps I end up doing a trivial linear solve and adding roundoff-level values to the existing radiation field. In the updated solver I skip those initial solves, meaning that I do not add roundoff-level values anywhere until the physics warrants changing the solution. It seems that this slight difference results in "Output times" that differ in the 14th digit.
The end result of this test is still an overall error (in comparison with the analytical solution) that is orders of magnitude more accurate than the tolerances. My recommendation is that we declare the new solver "successful", and just update the gold standard to reflect this changed digit.