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 scalefactor 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.

Gauge fields are present in the Standard Model and its extensions, for example the photon, W and Z bosons. They can have highly nonlinear dynamics during reheating, partly due to their Bosonic nature. Their simulations, however, are challenging -first because of the increase in the number of dynamical degrees of freedom, but more importantly, because of the requirement of preserving certain additional constraints (for example, the Gauss constraint) during the time evolution. In this work, we provide an algorithm and code, GFiRe -a Gauge Field Integrator for Reheating which is capable of simulating a system of charged scalar fields and Abelian Gauge fields on a lattice in an homogeneous and isotropic, expanding universe. GFiRe preserves the Gauss constraint and energy constraint exceptionally well; this remains true even when the expansion is determined self-consistently by the spatially averaged energy momentum tensor at each time step. Our time evolution is symplectic and explicit in time.
Below, we briefly review some relevant literature on nonlinear dynamics of cosmological fields (mostly in the context of reheating), and put our current work in context of the existing literature -highlighting the new and useful aspects of our present work.
Exploring the nonlinear dynamics of scalar fields during reheating is not new. Soon after the importance of non-adiabatic particle production in the oscillating inflaton background was appreciated (a phenomenon known as preheating [2][3][4][5]), it was realized that the subsequent backreaction can lead to a complex, non-linear evolution of the combined inflaton-daughter fields system. For some of the earliest scalar field simulations in the reheating context, see [68,69]. Since then, a number of dedicated lattice codes have been developed for scalar field models. Most of them are based on finite-differencing techniques, e.g., LatticeEasy [70], Defrost [71], HLattice [72], PyCool [73], Gabe [74], but there are also pseudo-spectral ones, e.g., PSpectRe [75] and Stella [76]. Usually, the codes solve the evolution of a system of spatially inhomogeneous, (self-)interacting scalar-field(s) in a spatially homogeneous and isotropic spacetime. In many cases, the expansion rate is determined self-consistently from the acceleration equation, with appropriate preservation of the energy constraint (i.e., the Friedmann equation is satisfied). 1 Beyond scalar fields, lattice simulations of gauge fields have also been carried out. These scenarios can be split into two groups depending on the nature of the scalar coupled to the gauge fields: (i) the scalar is uncharged (χ), and (ii) the scalar is charged (ϕ).
For an uncharged scalar + Abelian gauge fields systems, Lagrangians which include couplings of the form χFF [85][86][87][88][89][90][91], and f (χ)F 2 [92] have been simulated in the context of reheating. For non-Abelian fields, L ⊃ f (χ)trF 2 [93] has been explored. In [88][89][90], gauge invariance at the level of the spatially discretized action is guaranteed by the use of slightly modified Link Variables [94], and the Gauss constraint is preserved during the time evolution. The expansion of the universe is obtained from the spatially averaged energy density. However, their algorithm (which works in the temporal gauge), required that the integration in time is implicit. On the other hand, the authors in [85][86][87] use an explicit time evolution scheme. A Lorenz gauge is used, where all gauge degrees of freedom are made dynamical. Gauge invariance at the level of the spatially discretized action, and the preservation of the Lorenz gauge condition by time evolution is not guaranteed, however, the Lorenz gauge condition is carefully checked a-posteriori after solving the equations of motion on the lattice. The expansion of the universe is once again determined self-consistently.
For the case of charged scalar and gauge fields, lattice simulations with Abelian [27,[95][96][97][98] and non-Abelian [99][100][101][102][103][104][105][106] gauge fields have been carried where the spatially discretized action, and the time evolution respects the Gauss constraint by using standard Link Variables. However, these simulations assume that the scalar and gauge fields are spectators (i.e., energetically subdominant). The energy constraint in an expanding universe is not relevant in these studies since they assume that the Friedmann-Robertson-Walker (FRW) expansion is determined by an unrelated component of the Universe, e.g., a dark matter component or a thermal bath. The evolution of the scalefactor is set to a fixed power-law, a(τ ) ∝ τ p (for example, with p = 1 for radiation domination).
Our present algorithm and code, GFiRe, includes and adds to the desirable features that appear separately in the above mentioned works. In particular, some of its attractive features include: • Our discretized action respects a residual gauge invariance related to time-independent gauge transformations (we work in a gauge where the time component of the gauge field is set to zero). This is achieved by using Link Variables for the spatial discretization. The Gauss constraint follows naturally from residual symmetries of the spatially discretized action.
• Our time evolution scheme is explicit in time, and also symplectic. This scheme treats the time evolution of the scalefactor (whose equations of motion involve spatial averaging of the energy momentum tensor) in the same way it treats the time evolution of the fields. 2 • We show that our symplectic time evolution scheme (of the scalefactor and fields) preserves the Gauss constraint. We show explicitly that the Gauss constraint is respected exactly by non-linear symplectic integrators of arbitrary order. This result allows us to take advantage of the properties of high order symplectic integrators, for example, excellent stability and conservation for relatively large time steps without the need to worry about spurious gauge modes.
We note that when the charged scalar is the inflaton, or more generally, when the fields dominate the energy density of the universe, solving accurately for the FRW expansion and the equation-of-state can be important for the observational consequences of reheating [12,13,29,30]. Such calculations can substantially reduce uncertainties in inflationary observables such as the tensor-to-scalar ratio, r, and the spectral index n s [13]. Self-consistent expansion with good energy-constraint preservation can also be important for studies of non-gaussianity [22,23]. The paper is organized as follows. In Section 2 we set the scene for the symplectic integrators. After reviewing the generic equations of motion of the model, we rewrite the system in a Hamiltonian form and discretize it on a spatial lattice. In Section 3 we introduce the symplectic integrator algorithm and apply it to the Hamiltonian lattice system. We present the results from two test studies with our new algorithm in Section 4. Section 5 is devoted to concluding remarks and a discussion. Throughout the paper we use natural units in which c = = 1 and the reduced Planck mass is m Pl = 1/ √ 8πG. For the metric, we use the + − −− signature.

Scalar Electrodynamics
We consider classical scalar Electrodynamics, also known as the Abelian-Higgs or the Ginzburg-Landau model, with a general Higgs potential. The action for this theory is where ϕ is a complex valued scalar field, g is the determinant of the metric, and R is the Ricci scalar. The covariant derivative and the field tensor are where α(x ν ) is an arbitrary, real valued function of spacetime.
Varying the action with respect to g µν yields the Einstein Equations: where G µν is the Einstein tensor, R µν is the Ricci tensor and T µν is the energy momentum tensor that includes contributions from ϕ and A µ . Varying the action (2.1) with respect to ϕ and A µ yields the Euler Lagrange equations of motion for ϕ and A µ : In particular, the 0 component of the second equation was obtained by varying the action with respect to A 0 . This is the Gauss constraint. For future convenience we define Note that C G = 0 follows from the 0th component of the second line of (2.5). Moreover, even without using C G = 0, the equation of motion for ϕ and spatial part of the second line of (2.5), implies What if we had fixed A 0 = 0 at the level of the action? In this case, the Euler-Lagrange equations do not yield the Gauss constraint (C G = 0). Nevertheless, there is an alternate route to arrive at the correct expression for C G , and also show that ∂ 0 C G = 0. In this route, we can make use of the residual symmetry A i → A i + ∂ i α(x) and ϕ → ϕe −iα(x) of the action (left over even after A 0 = 0 is fixed) to arrive at an infinite number of Noether charges q(x) = C G with ∂ 0 C G = 0 on the equations of motion (see Ch. 10 in [109]). The key thing we miss out on is that we now need to set C G = 0 "by hand" at some initial time.
We now restrict ourselves to a homogeneous and isotropic universe where Here, τ is conformal time, and a(τ ) is the scalefactor. Furthermore, we chose to work in the temporal gauge with A 0 = 0. In this case, the equations of motion (2.5) become where the last equation is the Gauss constraint. The trace, and the 00 component of the Einstein equations (2.4) (assuming FRW spacetimes) yield 3 : where H ≡ a /a, and we have defined the spatially averaged density and pressure as ρ ≡ T 00 /a 2 and p ≡ δ ij T ij /(3a 2 ). In analogy with the Gauss constraint, the energy constraint There is an alternative treatment for arriving at the evolution equation for the scalefactor by first approximating the spacetime to be FRW at the level of the action (rather than at the level of the equations of motion). Explicitly, we first integrate the FRW Einstein-Hilbert term in the action by parts and then extremize the action with respect to the scalefactor, to arrive at the first equation in (2.11). Note that this time, the spatial averaging follows directly from extremizing the action.
In summary, we have shown that even if A 0 = 0 is enforced at the level of the action, we arrive at ∂ 0 C G = 0 via the equations of motion. Separately, assuming an FRW universe at the level of the action yields the same equation of motion for the scalefactor as the one obtained by imposing the FRW symmetries at the level of the equations of motion. Below, we will arrive at the same results in the Hamiltonian formulation. 3 For reference, (2.10)

Hamiltonian formalism in A 0 = 0 gauge and FRW spacetimes
Our algorithm and code is developed from a Hamiltonian formulation of the equations of motion, which we derive in this section. Here, we will work directly in the A 0 = 0 gauge, and assume an FRW universe from the beginning. We will also be laboriously explicit in our expressions since it will ease the route to their corresponding discretized versions.

Action and Hamiltonian
The action (2.1) becomes where the 0 subscript reminds us that we are in the temporal gauge and in an FRW universe. The complex field ϕ is written in terms of two real fields, (2.13) Our gauge choice allows us to apply the Hamiltonian formalism readily, since all variables are dynamical. Their conjugate momenta are The Hamiltonian is found from the Legendre transformation (2.17) 4 The Einstein-Hilbert term in the action, see Eqs. (2.12), in the FRW approximation involves a spatiallyindependent 3m 2 Pl a (τ ) 2 factor integrated over space. The factor can be moved in front of the 3-dimensional integral, yielding a formally divergent 3-volume integral over unity. The latter can be written as d 3 x = (2π) 3 δ 3 (0), after applying the Fourier representation of the Dirac delta function.

Hamilton's equations
The action in Eq. (2.12) is extremised when the Hamilton's equations of motion are satisfied. To write the equations of motion compactly, we use Poisson brackets {·, ·} P , defined in our case as The A and B are functions (not functionals) of the fields, fields momenta and their spatial derivaties, i.e., and we also define etc. We can now express Hamilton's equations as 21) and note that the Hamiltonian from Eq. (2.17) can be rewritten as (2.22) The first and third lines in Eq. gives the dynamical equations of motion for the real and imaginary parts of the Higgs and the spatial components of the gauge field, i.e., the first two lines in Eq. (2.9) (with A 0 = 0).

Gauss constraint
Where is the Gauss constraint in this Hamiltonian prescription? As discussed earlier, the residual symmetry A i → A i + ∂ i α(x) and ϕ → ϕe −iα(x) of the action (2.12) left over after setting A 0 = 0, leads to an infinite number of Noether charges which essentially leads us to the Gauss constraint [109]. Explicitly, the infinitesimal symmetry transformations are where |α(x)| 1. The conserved Noether current is given by The corresponding constant charge is (2.25) Since Q = const. for arbitrary α(x), then provided lim |x|→∞ α(x)π A (τ, x) = 0, which must hold, since α and π A vanish separately at infinity -the former according to the definition of a gauge transformation and the latter trivially. Hence, the gauge fixed action, (2.12), enforces ∂ τ C G = 0 and so do its Hamiltonian equations.
As a final remark we recall that when the Gauss constraint was derived in the case where the gauge was fixed after deriving the equations of motion, we ended up with C G (τ, x) = 0. In our Hamiltonian formulation, we arrived at ∂ τ C G (τ, x) = 0. We can readily restrict ourselves to C G (τ, x) = 0 at all times in our Hamiltonian formulation by simply imposing the Gauss constraint on the initial time slice, C G (τ in , x) = 0, and then use Hamiltonian equations to get C G (τ, x) = 0.

Hamiltonian formalism in
The next step towards solving the Abelian-Higgs model numerically is the discretization of the fields on three dimensional spatial lattice. We assume a comoving cubic grid with periodic boundary conditions and N 3 points. The comoving length of the edge of the unit cell is equal The scalar fields, {ϕ 1,x , ϕ 2,x }, are defined on the lattice points with discrete spatial coordinates x and the gauge fields, A j,x , are defined on the lattice links, connecting the adjacent points x and x + n j . The unit vectors are The standard links, U j,x , and plaquettes, P ij,x , are [94] Fig. 1 for a visual representation. We can now use these variables to define a lattice action in the A 0 = 0 gauge, which still has a residual time-independent gauge symmetry. Similarly to the continuous case from the previous section, we will show that in the Hamiltonian formalism, the lattice version of the Gauss constraint is recovered from the conserved charges corresponding to the residual gauge invariance of the lattice action.

Lattice action and Hamiltonian
The links and the plaquettes from Eq. (2.29) can be used to define spatial derivates on the lattice (no summation is assumed over repeated indices) Hence, the lattice version of the continuous action from Eq. (2.12) becomes where "c.c" refers to complex conjugate. Notice that the three-dimensional integral over space has become a summation over the lattice points d 3 x → N 3 x b 3 , and that the action has been rescaled by a factor of b 3 -a rescaling of the total action has no effect on the Euler-Lagrange or Hamilton's equations of motions. The continuous system of five fields capturing infinitely many degrees of freedom and the scalefactor has now become one having 5 × N 3 + 1 generalized coordinates, with conjugate momenta (c.f. (2.16)) The Legendre transformation yields the lattice Hamiltonian (c.f. (2.17)) (2.33)

Hamilton's equations
To extremise the lattice action in Eq. (2.31) we need to use the corresponding Hamilton's equations. We can again write them concisely in terms of the lattice Poisson brackets (c.f.
x ] T and all quantities are evaluated at equal times. A and B are functions of the lattice fields and lattice momenta, i.e., The lattice Hamilton's equations can be now expressed as where again the first and third lines are the definitions of the lattice conjugate momenta from Eq. (2.32), whilst the second and fourth lines are the evolution equations for the scalefactor and the lattice fields, respectively. We are again missing the (lattice version of the) Gauss constraint. Just like before, it can be recovered from a residual gauge symmetry of the lattice action, Eq. (2.31), as we show in the following section.

Gauss constraint
Thanks to the links and plaquettes the lattice action given in Eq. (2.31) is invariant under the time-independent symmetry transformations where In their infinitesimal form (α x 1), The Noether theorem then implies that for arbitrary α x , there is a conserved charge (2.39) It reduces to 5 (2.40) Note that J (τ ) = 0 holds for arbitrary α x . Hence, at each lattice point (2.41) This is the lattice counterpart to the Gauss constraint given in Eq. (2.26). Thus, the gaugefixed lattice action, Eq. (2.31), again enforces ∂ τ C latt G,x = 0 and so do its Hamilton's equations. The desired Abelian-Higgs theory with C latt G,x (τ ) = 0 is recovered by imposing the Gauss constraint on the initial conditions and then using the lattice Hamilton's equations to evolve the system.
The remaining question is whether there is a way to discretize in time the lattice fields, so that their evolution (according to the corresponding discretized in time lattice Hamilton's equations) still respects the Gauss constraint in Eq. (2.41). We show in the next section, that there is a scheme based on symplectic discretization which does that.

GFiRe implementation
As we argue below, we can numerically time-evolve the Abelian-Higgs system prescribed by the lattice action in (2.31) and meet the Gauss constraint (2.41) everywhere on the lattice using symplectic integrators.
In GFiRe, we use a symplectic prescription to integrate the Hamilton's equations on the lattice (Eqns. (2.36)), which we write compactly as , the integer j index of z runs over the 2 × 1 + 5 × N 3 phase space variables. In the rest of the section we give the details of the numerical integrator.

Symplectic integrators
We first present the structure of the symplectic integrator used by GFiRe in an abstract form. It might be useful to refer to Fig. 2 for a schematic picture of the integrator.
We need two assumptions about the form of the Hamiltonian. First, we assume that the Hamiltonian splits into three mutually non-commuting terms Second, we also assume that each Hamiltonian term is a function of only one variable from each pair of generalized coordinate and corresponding conjugate momentum. In other words, for a given j we have either ∂H i /∂q j = 0 or ∂H i /∂p j = 0 or both (here p j is the conjugate momentum to the generalised coordinate q j ). This assumption allows us to use standard symplectic integrator techniques [107].
The formal solution to the Hamilton's equations (3.1), can be written as 6 for an arbitrary time interval ∆τ . The numerical approximation for any finite time interval ∆τ is (see [107]): where ∆τ is the time step and K (k) is an operator involving Poisson brackets with the Hamiltonian terms from Eq. (3.2). The choice of the integrator order, k, determines the exact form of the operator. Thanks to our two assumptions, we can use a symplectic method of integration based on operator splitting techniques [107]. The integrator of lowest order, k = 2, is

5)
6 Note that in this formal solution, H is constant in time, since it is evaluated on the solution.
whereas higher-order operators are obtained recursively 7 The numerical solution can be made arbitrarily close to the true one as one increases k and/or decreases ∆τ . Hence, an integrator of arbitrary accuracy has been reduced to the operation of exp[K 1 ∆τ ], exp[K 2 ∆τ ] and exp[K 3 ∆τ ] in some particular combination and taking the appropriate values for the time interval ∆τ for each operation. Let's see how this can be implemented for the Abelian-Higgs system.

Evolution on the lattice
The lattice Hamiltonian, H latt 0 , from Eq. (2.33), splits into three non-commuting terms, just like in Eq. (3.2), Aj,x 2   , Note that each of the terms is a function either of a generalised coordinate or its conjugate momentum, but never of both consistent with our second assumption in the previous, more formal subsection.
We can now use Eqs. (3.4), (3.5) and (3.6) to evolve the complex scalar, 3-vector gauge field, scalefactor and their conjugate momenta. To this end we must find the action of each exp[K i ∆τ ] on the lattice phase space variables. In deriving the expressions below, we make an extensive use of the operator identity for our Hamiltonian (2.33): when applied to any z j . Note that this is not an approximation in the context of our Hamiltonian.
After acting with exp[K 1 ∆τ /2] on all z j we get (see Eq.

9)
7 There is a more optimal way of choosing the "weights" for higher order composite integrators. Such weights prevent the number of compositions from growing too rapidly as we go to higher and higher integrators (see [107]).
with the rest of z j unchanged. Similarly, the action of exp[K 2 ∆τ /2] yields leaving the rest of the generalized coordinates and momenta unchanged. Finally, the action of exp[K 3 ∆τ ] gives 11) with no changes in the generalised coordinates. Below, we state one of the main benefits of using the above scheme for time evolution.

Preservation of the Gauss constraint
We verify explicitly that after each individual step, Eq. (3.9) and/or Eq. (3.10) and/or Eq. (3.11), the lattice Gauss constraint, Eq. (2.41), is respected exactly, C latt G,x (τ ) = const. 8 To see this, first note that the action of exp[K 1 ∆τ /2] (see Eq. (3.9)) trivially yields C latt G,x → C latt G,x . The exp[K 2 ∆τ /2] step (see Eq. (3.10)) gives Finally, the exp[K 3 ∆τ ] (see Eq. (3.11)) operation yields where the ∂V terms cancel due to the U (1) symmetry of the potential, the links, U j , terms in the first two big brackets cancel with the links terms in the other two big brackets, whereas the sums of the Plaquette terms vanish by virtue of the identity 1,2,3 j,l iP jl,x + c.c. = 0.
(also see Eq. (2.29)). Another, perhaps more impressive way of stating our result is that when we evaluate C latt G at τ + ∆τ by using the above expressions for time evolved z j , all terms proportional to ∆τ n for arbitrary n vanish. Thus, regardless of the order of the time integrator, k, the Gauss constraint is always satisfied to machine precision at each lattice site. 9 As we showed above, for physical consistency, initially we must set the constant on the right hand side in the Gauss constraint equation, Eq. (2.41), to zero everywhere on the lattice. The symplectic evolution now guarantees that it remains zero at later times.

Approximate post-inflationary initial conditions
Within the inflationary paradigm, the (almost) homogeneous inflaton field dominates the Universe during and shortly after the end of inflation. One can assume that the inflaton is the real part of the complex salar (i.e., slow-roll inflation is realized along the real axis in the complex field space, towards the minimum of the scalar field potential, see Fig. 3). At the beginning of our simulations of preheating, the initial variables on the lattice are ϕ 1,x =φ 1 + δϕ 1,x , ϕ 2,x = δϕ 2,x , π 1,x =π 1 + δπ 1,x , π 2,x = δπ 2,x , where the choice for the initial value of a is conventional and the expression for π a is simply the Friedmann equation. At later times, we use this expression to check the 'energy conservation' of our integrators. We observe violations that scale correctly with the time step -O(∆τ k ) for k-th order integrators. Theφ 1 andπ 1 variables are determined from the homogeneous inflationary dynamics, whereas the δ-field perturbations have a power spectrum set by the quantum vacuum fluctuations. The initial fluctuations for each lattice field can be expanded in terms of Fourier modes. 10 These can be then written in terms of products of mode functions and 'stochastic' complex numbers. The stochasticity in the initial field perturbations is what allows the classical simulations to approximately capture aspects of the quantum uncertainty in the actual vacuum fluctuations. We also must satisfy the Gauss constraint, Eq. (2.41), with the constant term (3.16) The random ('stochastic') complex numbers provide the classical counterpart of the quantum creation and annihilation operators, and take the values Here X I k is a uniform deviate on (0, 1) and Y I k is a uniform deviate on [0, 1). The mode 10 The lattice Fourier conventions we use are Figure 3. The complex scalar potential profile, see Eq. (4.1), for two separate cases v = 0 (left panel) and v = 0 (right panel). After inflation ends, the field is primarily rolling down along ϕ 1 -axis, followed by oscillations about the minimum, which can lead to resonant particle production in the gauge fields as well as the ϕ field. Note that in both cases, we begin withφ 1 (τ in ) v.
functions are given by the flat spacetime expressions which approximate well the modes of interest, i.e., subhorizon modes, k H, at the end of inflation, regardless of the coupling strength, parametrised by the comoving Compton wavenumber, k C ≡ eφ 1 , for the gauge fields. Note that the normalization factor of (2π) 3/2 /L 3/2 lat in Eq. (3.16) is needed to make the (initial) two-point functions of field perturbations, averaged over the lattice volume, independent of the comoving box size, L lat = N b, and equal to the continuous ones, see for example, [70].

Numerical studies of the scalar electrodynamics
To test our algorithm, we study the non-linear dynamics in models with the scalar field potential (see Fig. 3) Note that this setup corresponds to the well-known quartic inflationary scenario [110,111]. As such, this particular shape of the potential is in conflict with CMB observations [112]. Nevertheless, since the main focus of this section is on the post-inflationary dynamics, we content ourselves with interpreting Eq. (4.1) as the effective form of the inflaton potential after inflation, remaining agnostic about the details of V during inflation. For simplicity, we also setπ 1 = 0 initially. Shortly after the end of inflationφ 1 (τ ) begins to oscillate. Its oscillations can amplify exponentially fast the initial field fluctuations, δϕ 1 , δϕ 2 and δA. The phenomenon can be understood in terms of parametric resonance [113] and is more generally known in the context of inflation as preheating [5]. The field fluctuations can grow very quickly, until mode-mode couplings become non-negligible. This marks the end of the linear stage of preheating [9] and the onset of backreaction. The ensuing non-linear dynamics can be quite rich and is strongly dependent on v. We now consider two different scenarios -v = 0 and v = 0, separately and study the linear and non-linear regimes with our numerical prescription.
For the results below, we typically rely on a N 3 = 512 3 lattice. We use a fourth order symplectic integrator (k = 2 in the notation of Eq. (3.6)) with a time step ∆τ = 0.028( √ λφ 0 ) −1 and comoving edge length L latt = 60( √ λφ 0 ) −1 . The comoving lattice spacing ∆x = L latt /N . Note that we always have ∆τ < ∆x. By the end of the simulations a(τ f )∆x < 1/( √ λv) where √ λv is the mass scale in the valley of the potential in the v = 0 case. In the v = 0 case, the effective inverse mass scale case grows with the scalefactor, hence we always resolve the relevant spatial and temporal scales at late times, if we resolve them initially.

Nonlinear dynamics in the v = 0 case
We begin with v = 0, for which the scalar field potential profile takes the form of a quartic bowl, see Fig. 3. At the start of the simulations, the background scalar field oscillates along the real axis in the complex ϕ plane, as shown in the left panel of the figure. Depending on the ratio e 2 /λ, these oscillations can lead to the exponential amplification of the complex scalar and/or the gauge field fluctuations. We set e = √ λ . This choice is rather interesting, since it entails significant gauge field, δA, resonant particle production, but virtually no Higgs particle production 11,12 [113]. Our lattice simulations capture this subtle effect.
Field Power Spectra: In Fig. 4, we plot the evolution of the gauge field power spectrum. Soon after the simulations commence, there is a rapid excitation of a broad range of δA comoving wavenumbers, as depicted by the dark red curves in the figure. The power spectrum eventually stops growing, when backreaction kicks in. Later on, non-linear effects such as rescattering excite even a wider range of comoving modes and drive the gauge field power spectrum into a long-lived single-broad-peak configuration.
Unlike the gauge field power spectrum, the power spectra of the two scalar field components, shown in Fig. 5, do not feature the early rapid growth, as expected for our parameter choice, Eq. (4.3). A broad range of comoving ϕ modes starts getting excited only around the time of backreaction of δA onφ 1 , as depicted by the light red and orange curves in Fig. 5. At late times, mode-mode couplings again drive the ϕ power spectra into a long-lived broad single-peaked configuration.
Energy fraction and equation of state evolution: The two-stage evolution (an initialφ 1 oscillatory stage with δA growing, followed by a long-lived steady-state nonlinear stage) can be also observed in the evolution of the fractional energies given in the left panel in Fig. 6. For τ between 0 and ∼ 70( √ λφ 0 ) −1 , most of the total energy, E tot is stored in the oscillatingφ 1 , in the form of kinetic and potential energy. Afterwards, as non-linear effects become important the energy gets quickly redistributed across all components. At around τ ∼ 200( √ λφ 0 ) −1 , the system enters the long-lived non-linear stage, characterized by a steady energy equipartition [96]. Thereafter, the ϕ self-interaction potential energy vanishes, f pot ≈ 0, whereas f kin ≈ f grad and f elec ≈ f magn . Our numerical algorithm also allowed us, for the first time, to compute the self-consistent expansion of the FRW background sourced by an inhomogeneous scalar electrodynamics system. The evolution of the mean equation of state, w, is shown in the right panel in Fig. 6. We find that w = 1/3 throughout, which is to be expected. During the initial oscillatory stage, we have a single homogeneous oscillating scalar field with a quartic self-interaction, which implies the well-known result of 1/3 for the mean equation of state [115] (in our notation, the only non-zero f i s are f K = 2f V = 2/3). Later on, since the ϕ self-interaction potential energy vanishes, the real and imaginary parts of ϕ, as well as the components of the gauge fields behave as massless radiation, again implying a radiation-like equation of state (in our notation in the radiation limit, f K + f elec ≈ 1/2, since the magnetic and electromagnetic components are approximately equal, as well as the Higgs kinetic and gradient energies).
Lattice snapshots: Individual snapshots of the field configurations and their energy densities on the lattice at any given time reveal a rich spatial structure in the fields at both the linear and nonlinear stages. In Fig. 7, we provide an example of snapshots of the fractional electric and magnetic field densities at three different times. The initial resonance instability leads to a growth of large length-scale modes with a somewhat larger fraction in electric fields. The third panel reveals a more scrambled configuration at late times (after backreaction). While we do not do so here, plotting the vector field configurations (rather than scalar energy densities), or pseudoscalar quantities such as (E · B) also provides useful insight into the complex underlying dynamics.
Energy and Gauss constraint preservation: To keep track of the violation of the energy conservation in our simulations, we consider the quantity where C E was defined in section Eq. (2.11). For the simulation whose results we have been discussing so far, the evolution of E is shown in the left panel in Fig. 8. Note that it is easy to achieve a very small degree of energy violation, < 10 −5 , with a fairly large time step, due to the high order of the symplectic time integrator. Furthermore, the energy violation is quite stable and grows very slowly, due to the time-reversability of symplectic integrators.  .7) for the case where v = 0. The energy violation is very stable which is a generic property of symplectic integrators. The Gauss constraint violation starts growing from being close to machine precision initially due to finite differencing noise during the initial homogeneous oscillatory stage. As the Higgs fragments, the violation settles down to a constant small value. Note that we have defined E and G to be positive definite.
We characterize the violation of the Gauss constraint with the quantity (see Eq. (2.26)) (4.7) We show its evolution at an arbitrary lattice point, x, in the right panel in Fig. 8. The algorithm performance is excellent, with the violation never exceeding 10 −6 and for most of the time remaining close to machine precision. The brief increase in the violation is observed only during the oscillatory phase, while the ϕ field is still homogeneous. During this period, the calculation of spatial derivatives and/or differences of products of fields might lead to numerical errors due to attempting to compute small differences between large numbers, a phenomenon known as differencing noise. The observed growth in the violation during the oscillatory stage could be explained by such numerical errors.

Nonlinear dynamics in the v = 0 case
We move on to the v = 0 case, for which the profile of the Higgs potential, V (ϕ), resembles a sombrero hat, see right panel in Fig. 3. We set v = 1.32 × 10 −2 m Pl , amplitude ofφ 1 oscillations, see Eq. (4.2), the initial parametric resonance phase is unaffected by v. We still have significant δA resonant particle production. Again parametric resonance does not develop in the Higgs due to our choice of e, as explained in Section 4.1. Only once δA begins to backreact, there is significant amplification of a broad range of comoving Higgs modes. After backreaction, the power spectra of the Higgs and the gauge fields again settle into stable broad single-peaked configurations. Since the power spectra plot are qualitatively similar to the v = 0 case, we have relegated them to an appendix.
Cosmic strings: Plotting the evolution of the fields in real space, reveals a phenomenon that cannot be picked out from the evolution of the power spectra. Note that the v = 0 Higgs potential (right panel in Fig. 3), can support the non-trivial field configurations known as topological strings [116]. They can be generated during thermal phase transitions via the Kibble mechanism in the form of cosmic string networks (for reviews see, e.g., [14,15,117]). Strings can be also produced after backreaction due to parametric resonance [27,[118][119][120], just like in our case. Since strings are characterized by a non-zero integer topological number, known as the winding number, n, we plot the lattice points with n = 0 at four different times in Fig. 9. The first panel in Fig. 9 is at the start of the simulation. All lattice points have n = 0, consistent with the inflationary initial conditions, see Eqs. (3.14) and (3.16). Towards the end of the resonant particle production and the onset of backreaction we observe copious formation of strings and string loops with a sub-Hubble correlation length, as shown in the second panel in Fig. 9. The strings then interact, 13 reconnect into loops and gradually evaporate via classical radiation. We see features developing on loops, which split from the larger loop to form smaller loops, which then decay away. The last large loop in our simulation is seen (third panel in Fig. 9) around the time the fields enter the long-lived steady-state non-linear stage. At this stage, the string core is resolved by roughly 10 points per linear dimension. A more dedicated recent study of decay of single loops can be found in [121]. At late times, there are no strings in our box, which is consistent with the conservation of the net winding number and our initial conditions having zero n. 14 We note that the string configurations never dominate the energy density in our box.
Energy fraction and equation of state evolution: The evolution of the mean fractional energies, f i , is shown in the left panel in Fig. 10. Just like in the case with zero v from Section 4.1, we again have two distinct regimes with a brief transitionary period inbetween. During the initial oscillatory phase, most of the energy is stored in the oscillatingφ 1 , again with f kin ≈ 2f pot . This approximate equality (which is exact for a quartic Higgs potential) becomes less accurate with time, since the amplitude ofφ 1 ∝ a −1 and the v 2 term in V becomes 14 Another late-time configuration which is consistent with our initial conditions and the conservation of the topological charge is a pair of parallel strings with winding numbers of unlike signs (and therefore zero net winding number). We observed such final field configurations with both strings stretched across the same pair of opposite faces of the cubic box. The two strings were stationary with respect to the comoving lattice and not reconnecting into loops and evaporating for the duration of the simulations. The probability for such scenarios was quite low, < 10%, for an ensemble of initial field realizations, see Eq. (3.17), and non-vanishing only for very small box sizes, much smaller than the Hubble radius around the time of backreaction and string formation. increasingly important. Backreaction occurs around τ ∼ 80( √ λφ 0 ) −1 and is followed by a swift redistribution of energy. By τ ∼ 200( √ λφ 0 ) −1 , the last string loops start evaporating and the fields settle on a long and steady approach to equipartition.
After this point, at almost every location on the lattice, the ϕ field has settled into the valley V (|ϕ| = v/ √ 2). The radial component of the ϕ field oscillates above the vev with amplitude v in a spatially inhomogeneous manner. Note that the radial Higgs component behaves as a massive scalar, since the bottom of the valley at |ϕ| = v/ √ 2 is quadratic. The azimuthal Higgs component (the massless Nambu-Goldstone boson) is 'eaten up' by the gauge field like in the conventional Higgs mechanism, rendering the gauge field a massive vector field. Hence, both the radial Higgs degree of freedom and the gauge field behave as massive fields. They have massive non-relativistic modes, whose energy density redshifts as ∝ a −3 , slower than the energy density of the relativistic modes (which scales as ∝ a −4 ). This implies that after a few e-folds of expansion, both the radial Higgs component and the gauge field should behave as pressureless matter, with potential and kinetic energies equal to each other and much greater than the gradient and curl energies. In our notation, this means that for the radial component of the Higgs, f kin → f pot , whereas for the gauge field f elec → f grad (recall that the 'potential' mass term for the gauge field is contained in f grad , see Eq. (4.4)) and f magn → 0. Indeed, this is the late-time behaviour shown in the left panel in Fig. 10.
The self-consistent evolution of the mean equation of state, w (see Eq. (4.5) for its definition) is given for the first time in the right panel in Fig. 10. Initially, it has a radiationlike value, due to the oscillations of the backgroundφ 1 . The growing deviation from 1/3 is due to the v 2 in V becoming increasingly important as the amplitude of the oscillatingφ 1 is redshifted. After backreaction, w steadily approaches 0, since both the gauge field and the radial component of the complex Higgs field are massive. Their relativistic pressures are redshifted away, leaving behind pressureless scalar and vector dust with w → 0. Energy and Gauss Constraint: The energy conservation is shown in the left panel in Fig.  11. We have used the same lattice parameters as for the simulation from Section 4.1. The energy conservation is still excellent, ≤ 10 −4 . It is almost identical to the one for the v = 0 case given in the left panel in Fig. 8, worsening only slightly at late times. The reason for this slightly worse performance for v = 0 can be traced back to the fact that we work with a fixed conformal-time step, ∆τ . For v = 0, i.e., a quartic Higgs potential, the typical frequency scales always decrease with time as ∝ a −1 , implying that their product with the cosmic-time step, a(τ )∆τ , is constant.
On the other hand, for the massive case, v = 0, the typical frequency scales are constant, implying that their product with the cosmic-time step grows like ∝ a(τ ), thereby increasing the time-integration error. Even though it was not necessary for this study, this small degradation in energy conservation can be easily alleviated by decreasing the conformal-time step only slightly. This takes advantage of the fact that the order of the time integrator, k, is high and the energy conservation is quite sensitive to the time step. The local truncation error in the time integration is O(∆τ k+1 ), see Eq. (3.6), and the total accumulated error is O(∆τ k ). For k = 4, the energy conservation can be improved by one or two orders of magnitude, when we decrease the conformal-time step only by a factor of 10 1/4 or 10 1/2 , respectively, as shown in Fig. 12.
The Gauss constraint violation is shown in the right panel in Fig. 11 for a random lattice point. We find that it is qualitatively identical to the one for the v = 0 case given in Fig. 8. It is also insensitive to the conformal-time step, which is expected for a quantity dominated by differencing noise.

Discussion
We have presented a novel prescription for numerically evolving Scalar Electrodynamics in FRW spacetime (sometimes also referred to Abelian-Higgs system, or explicitly, a charged scalar minimally coupled to Abelian gauge fields) . Our prescription combines two different well-known techniques in numerical simulations of related systems, one for handling the spatial derivatives, and another for temporal derivatives. The spatial discretization is carried out according to the Lattice Gauge Field theory prescription (using the lattice links formalism). The time evolution is performed with symplectic integrators. The algorithm allows for the selfconsistent evolution of the FRW scalefactor and fields, while respecting the Gauss constraint exactly on the discretized lattice. To the best of our knowledge, this is the first explicit-intime algorithm which guarantees the preservation of of the Gauss constraint, while solving for the expansion of the universe self-consistently. In addition, the time integrator can be made of arbitrary high order, without violating any of the algorithm's properties. 15 The use of symplectic integrators was inspired by already existing scalar field lattice codes such as HLattice [72], Defrost [71] and PyCool [73]. They are known for their excellent energy conservation and stability, when solving for the expansion of the universe self-consistently. Perhaps, one of the reasons why such integrators have not been employed in gauge field studies is the uncertainty of whether they will respect the Gauss constraint equation -a non-dynamical equation which always appears in gauge field theories, but never in pure-scalar models.
The issue with the Gauss constraint equation has long been solved in flat spacetime, by discretizing the Higgs field on a 4d rigid lattice, and defining the gauge fields on the lattice links. This method has been successfully extended to FRW spacetimes, with the scalefactor evolving according to some fixed power-law, e.g., a(τ ) ∝ τ , not determined by the evolution of the Higgs and gauge fields [27,96]. 16 Given these successes, we employed a combination of the two approaches to achieve both good self-consistent expansion and preservation of the Gauss constraint. We showed how to combine them in a very straightforward way. To use symplectic integrators, we needed a simple Hamiltonian. The Hamiltonian for Scalar Electrodynamics was simplified substantially by working in the A 0 = 0 gauge. This gauge choice made all kinetic terms canonical, which in turn gave use a symplectic integrator for the time evolution. The remaining spatial derivatives were discretized by defining the complex scalar and the 3-components of the gauge field on a 3d spatial lattice, similarly to Lattice Gauge Field theory (with the use of Link variables). This automatically allowed the system to respect the residual gauge freedom in A 0 = 0 gauge. More importantly, this spatial discretization yielded a well-defined expression for the lattice version of the Gauss constraint, consistent with the residual symmetries. This version of the Gauss constraint turned out to be respected exactly by the symplectic integrators.
Numerical Studies: To test our algorithm and code, we investigated two different reheating scenarios. In both, the role of the oscillating inflaton was played by the real component of the complex scalar field. In the first scenario, the scalar potential was a simple quartic minimum (i.e., 'unbroken'), whereas in the second one, the scalar field potential allowed for spontaneous symmetry breaking (a Higgs-like potential). The lattice simulations captured the initial preheating phase in which the gauge field was excited non-perturbatively due to parametric resonance in the oscillating homogeneous scalar field background. They also revealed the subsequent non-linear stage after backreaction and fragmentation of the condensate. All qualitative predictions of linear analyses for the resonant particle production were reproduced. The subsequent non-linear stage included interesting non-linear phenomena such as the formation and disappearance of Nielsen-Olesen strings (in the case with the symmetry breaking potential).
The self-consistent FRW expansion was computed throughout -the energy conservation violation was stable, being always ≤ 10 −4 and scaling appropriately with the size of the time step, ∝ ∆τ k for a kth-order integrator. As expected, the case with unbroken scalar field potential lead to a radiation-like equation of state during the oscillatory as well as the non-linear stages. Similarly, a late-time matter-like state of expansion predicted for the spontaneously broken case (since along with the radial part of the scalar, the gauge fields are now massive), was reproduced during the non-linear stage. In both cases the Gauss constraint was preserved, with violations < 10 −6 ; these violations were dominated by numerical errors due to subtraction of large quantities with very small differences (differencing noise).
Limitations: We note that the combination of a symplectic time integrator and a Lattice Gauge Field type of discretization does not seem to work for non-Abelian gauge theories, e.g., an SU (2) Yang-Mills theory with a Higgs doublet. In this case, the Gauss constraint is not respected by the symplectic integrator. A reason for that could be the non-linear nature of the Gauss constraint in non-Abelian theories. Another class of gauge field models which is not well-suited for our prescription is the one featuring an axion, χ. Both Abelian and non-Abelian theories with a Chern-Simons type of interaction, ∝ χFF , cannot be integrated symplectically due to the non-canonical structure of the gauge field kinetic term. However, an interaction of the form ∝ χF 2 could work at least in Abelian theories, since the kinetic term of the gauge field is canonical up to an axion-dependent rescaling, meeting the conditions for usage of symplectic integrators.
We hope that GFiRe will be used for many more cosmological studies of non-linear gauge field dynamics. We plan to make GFiRe public in the near future.
A Power Spectra for v = 0 case We provide the time evolution of the power spectra for the gauge fields as well as the real and imaginary components of the ϕ field below in the v = 0 case. See corresponding discussion in the main text.  Figure 13. The evolution of the power spectra is qualitatively the same as in the v = 0 case.