The time integration scheme is just one part of achieving second order convergence of the solution. The second part is implementing higher-order integration schemes for the heat source. A future PR will address this, and may be focused on using quadrature methods for the volume integrals of the source term. An example time-step convergence plot is shown with the time-integration for the thermal solve and the quadrature-based heat source integration used together:
Originally posted by @colemanjs in #28 (comment)