GFiRe: a Gauge Field integrator for Reheating
  We present a new numerical algorithm and code, ${\sf GFiRe}$, for solving the
non-linear evolution of Abelian gauge fields coupled to complex scalar fields
in homogeneous and isotropic spacetimes. We adopt a hybrid approach to solving
the system: the spatial derivatives are discretized using standard Lattice
Gauge Field Theory techniques, whereas the time evolution of the fields and
scale factor is implemented with explicit, composite, symplectic integrators.
An important property of our compound algorithm is that the discretized Gauss
constraint is respected exactly, regardless of the order of the symplectic
integrator. This remains true even when the background expansion is computed
"self-consistently"; that is, when the expansion history is computed using
spatial averaged components of the energy momentum tensor in the simulation
volume. Hence, our code can also be used in cases where the fields dominate the
energy density of the universe, for example, during reheating after inflation.
  We test the algorithm in scenarios of reheating where the inflaton is a
complex scalar field with a potential $\propto(2|\varphi|^2-v^2)^2$ and is
coupled to an Abelian gauge field. Tracing the evolution of the system through
complex dynamics (including resonant excitation of fields, backreaction,
formation of solitons, and changes in the equation of state) in a
self-consistently expanding universe, we find the energy conservation violation
($<10^{-4}$) to be very stable and the Gauss constraint violation ($<10^{-6}$)
to be dominated by differencing noise.
