Gauge preheating with full general relativity

We study gauge preheating following pseudoscalar-driven inflation in full general relativity. We implement the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) scheme to solve the full nonlinear evolution of the metric alongside the dynamics of the pseudoscalar and gauge fields. The dynamics of the background and emission of gravitational waves are broadly consistent with simulations in a Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime. We find large, localized overdensities in the BSSN simulations of order δ = δρ/ρ ∼ 30, and the dimensionless power spectrum of δ peaks above unity. These overdense regions are seeded on length scales only slightly smaller than the horizon, and have a compactness C ∼ 0.1. The scale of peak compactness is shorter than the Jeans length, which implies that pressure of the matter fields plays an important role in the evolution of these objects.


Introduction
The end of cosmic inflation [1][2][3][4] and the subsequent transition to the radiation-dominated hot Big Bang remains one of the most poorly understood epochs in the evolution of the Universe.During this reheating epoch, the accelerated expansion of inflation must end, and the inflaton energy density must be eventually transferred to relativistic degrees of freedom to begin the hot Big Bang.Because the scales involved are so small-smaller than the Hubble radius at the end of inflation-information about the reheating epoch is either erased as the resulting Standard Model plasma achieves thermal equilibrium or is inseparable from the effects of the subsequent nonlinear gravitational evolution of structure formation.
The strongest effects of gauge field production occur when the inflaton rolls the fastest, which typically occurs at the end of inflation.The large couplings required to produce effects observable on scales accessible to the cosmic microwave background (CMB) or gravitationalwave interferometers subsequently generate a high-frequency gravitational wave background during preheating so large that existing bounds on the effective number of relativistic species rules the regime out [86][87][88]. 2 The production of large metric perturbations calls into question whether nonlinear gravitational effects can be safely neglected.Since gravitational wave backgrounds from gauge preheating currently provide the strongest constraint on these models, it is crucial to test the robustness of prior predictions [86][87][88] and whether nonlinear gravity might enhance [96] or suppress [97] the production of gravitational waves.Further, the production of large metric fluctuations indicates the presence of large inhomogeneities in the matter sector which may undergo gravitational collapse under the influence of local gravity.The purpose of this work is to push further the exploration of gravitational signatures of preheating into the regime of local, nonlinear gravity.To that end, building on the results of [98], we initiate a study of preheating into gauge fields using the full machinery of numerical relativity [99][100][101] to follow the nonlinear evolution of the metric.
In this paper, we demonstrate that preheating probes regions where nonlinear gravitational effects are expected to become important, and we explore the degree to which simulations that include nonlinear gravity are required to accurately characterize this epoch.In particular, we perform a careful study of gravitational effects during gauge preheating, focusing on the development of large density inhomogeneities and gravitational wave production.We first demonstrate that linearized gravity is quickly violated.Then, using numerical relativity, we follow the evolution of the metric including the effects of dynamical gravity to show that this breakdown of linearized gravity does not signal the formation of black holes.We demonstrate that, despite producing regions with very large density contrasts δρ/ρ ∼ 30, there is nevertheless no evidence that black holes are formed-no horizons are formed in our simulations.We show that the spatial extent of the regions of large density contrast are smaller than the Jeans length, which indicates that pressure plays an important role in their subsequent evolution.Furthermore, their compactness, which is a measure of whether a region satisfies the so-called 'hoop-conjecture' [102], attains a maximum value of C ∼ 10 −1 .
The paper is structured as follows.In section 2 we define the model and present the equations that govern the evolution and amplification of gauge fields during and after inflation.In section 3 we present results of numerical simulations, focusing on the dynamics of density and gravitational wave fluctuations as a function of the axion-gauge coupling strength.Our conclusions and proposed avenues for future work are presented in section 4. The details of the decomposition of gauge fields in the BSSN formalism are relegated to appendix A, while appendix B details how our initial conditions are set in perturbation theory in the BSSN formalism, and finally appendix C describes robustness checks of our results.
We use natural units, which set ℏ = c = 1, and define the reduced Planck mass M Pl = 1/ √ 8πG.Repeated/contracted Greek spacetime indices are summed via the Einstein summation convention.

Gauge preheating and the BSSN formalism
We consider a pseudoscalar inflaton, φ, minimally coupled to Einstein gravity, and coupled to a U(1) gauge field, A µ , described by the action where R is the Ricci scalar and F µν = ∇ µ A ν − ∇ ν A µ is the field strength tensor whose dual is The Levi-Civita tensor is which is written in terms of the permutation symbol (the Levi-Civita symbol), We work with the convention that ϵ µ              Though inflation driven by a monomial potential is strongly disfavored [104,105], we choose a quadratic potential merely to provide a simple model for the inflaton's dynamics during preheating.For further discussion of the effects of potential choice on gauge preheating, see [87].We consider the standard shift-symmetric, dimension-5 axial coupling between the pseudoscalar and the gauge field Numerical investigations of preheating typically employ a homogeneously expanding, Friedmann-Lemaître-Robertson-Walker (FLRW) background spacetime with metric Here τ is conformal time, and the scale factor a(τ ) evolves according to the Friedmann equations and Primes denote a derivative with respect to conformal time, τ , and overbars generally indicate averaged quantities on constant-τ hypersurfaces.The (averaged) energy density and pressure are ρ ≡ − T 0 0 and P ≡ δ j i T i j /3, respectively.In an FLRW spacetime, the equations of motion for the scalar field and gauge fields are In these equations, repeated Latin (i.e., spatial) indices are implicitly contracted with the Kronecker delta function.Fixing the (flat-space) Lorenz gauge η µν ∂ µ A ν = 0, eq.(2.9b) and eq.(2.9c) may both be recast into the second-order differential equations taking the form (2.10)

Numerical Relativity
To implement nonlinear gravity, we use the BSSN decomposition [99,100] where the metric is decomposed as The lapse α and the shift β i parameterize (nondynamical) gauge degrees of freedom.Threedimensional hypersurfaces are measured by the spatial metric3 which is further decomposed into a conformal factor ϕ and a unit-determinant spatial metric γij .We denote spatial covariant derivatives (those compatible with the spatial metric γ ij ) with D i , and indicate trace-removed quantities as Three-dimensional hypersurfaces are defined relative to the temporal coordinate t with normal vector The evolution of the metric components is specified by the set of first-order differential equations, where D i is the 3-dimensional covariant derivative, K is the trace of the extrinsic curvature tensor, Ãij is the traceless part of the extrinsic curvature, see eqs.(A.7) and (A.8), and R ij is the spatial projection of the Ricci tensor.The sources, ρ, S, and S ij , are calculated from the stress-energy tensor, see appendix A.1 for details.The BSSN system introduces new degrees of freedom so that Einstein's equations are computationally stable.Numerical solutions must satisfy the Hamiltonian and momentum constraints, and to sufficient precision throughout the simulations. 4ince the lapse α and shift β i are purely gauge degrees of freedom, we are free to specify them via convenient evolution equations.To allow for the formation of compact structures, such as black holes, while also trying to stay near the FLRW background [106] we take a Bona-Massó slicing condition for the lapse [101,107], as well as a hyperbolic gamma driver condition for the shift, where we set η = 100 as a choice that minimizes the violation of the constraints eqs.(2.16) and (2.17) at late times (see appendix C).In the homogeneous limit, this slicing reduces to which is the normal conformal-time evolution equation for the scale factor when α → a.
We track the scalar field φ and its conjugate momentum, Π, which is the component of its (covariant) four-gradient normal to spatial hypersurfaces, We additionally promote the spatial derivatives of the scalar field to dynamical quantities, and split the vector potential A µ into its components along and orthogonal to spatial hypersurfaces, and respectively, such that standard four-potential is simply reconstructed as The electric and magnetic fields are then given by In practice, we can fully evolve the system by evolving A, A m , and the purely spatial vector E m .In terms of these variables, the scalar field's Euler-Lagrange equation, eq.(A.10), reduces to the first-order system In the BSSN system, the metric in eq.(2.11) is not conformally related to η µν -unlike the FLRW metric eq.(2.6)-and it is more convenient to chose the covariant Lorenz gauge (2.29) The auxiliary field Z is a dynamical constraint-damping degree of freedom that can increase the computational stability of the system [108][109][110][111].In practice, we do not find including Z meaningfully improves the stability of our simulations, but we retain it below for completeness.With this choice, the equations of motion for the gauge field sector are (2.30d) Here the "source vector" has components (2.31b) Gauss's law requires the divergence of the electric field satisfy which is an additional constraint on the system that must be satisfied throughout the simulation, see appendix C. Finally the source-terms in eq. ( 2.15) evaluate to ) ) which close the dynamical system.

Results
In this section we present the results of our simulations.We begin in section 3.1 by describing our software and the computational choices we make in our simulations.We then compare the evolution of the background, including the full effects of nonlinear gravity, with our previous FLRW simulations.We then look for signs that gravity might lead to collapse in section 3.2, and finally study the resulting gravitational wave spectra in section 3.3.

GABERel, initial conditions and background evolution
We extend GABERel [98,112] to treat gauge fields in addition to scalar fields (as presented in Ref. [98]) in full numerical relativity using the BSSN formalism.Such an analysis is the only way to assess whether nonlinear gravitational physics is important and whether it leads to the collapse of gravitationally bound objects.
The BSSN simulations we present here use grids with N 3 = 384 3 points and a comoving box length L = 7.5 m −1 , with a = 1 at the end of inflation; therefore the initial, physical box size is L 0 = e −2 L at the beginning of the simulation.We use the standard fourthorder Runge-Kutta method for time integration with steps of size ∆t = ∆x/10 = L 0 /(10N ).The spatial discretization uses fourth-order finite-difference stencils (with upwind variants for advective derivatives and centered differences otherwise).We begin simulations two e-folds before the end of inflation, which ensures that even the largest-scale modes in the simulation begin with nearly Bunch-Davies initial conditions.Each mode has a uniform-random phase and an amplitude sampled from the Rayleigh distribution with variance set by the Bunch-Davies vacuum.The fields' time derivatives are set in the Wentzel-Kramers-Brillouin (WKB) approximation following the standard prescription [98,113].We filter out high-wavenumber modes from the initial conditions as in Ref. [98], setting the cutoff scale to We can also compare the fully nonlinear system to the FLRW system from previous work.To simulate the FLRW system, we use the same software described in Refs.[87,88], which also uses a fourth-order spatial discretization and time evolution scheme.Because FLRW simulations are less computationally expensive, we use grids with N 3 = 512 3 points and a larger box length of L = 15 m −1 .The initial conditions are cut off at a wavenumber The FLRW simulations are therefore rather different from the BSSN ones, not just in the physical content and its representation but also in numerical implementation.The only quantitative comparisons that can be made between the two are statistical, for which reason the differing choices of simulation volumes and grids are relatively inconsequential.Indeed, agreement between the two methods (in regimes where it is expected) provides a robust test of the results' independence of the numerical procedure.
The axion-gauge field coupling in eq.(2.9a) leads to a tachyonic instability in the gauge fields whenever φ0 ̸ = 0 for momenta that satisfy [46,48,57] During inflation, this tachyonic instability leads to the exponential enhancement of one (helical) polarization of the gauge field relative to the other.Since the width of the instability in eq.(3.1) is proportional to the axion velocity, the largest effects typically occur near the end of inflation where the inflaton velocity is the largest.These dynamics complicate the setting of initial conditions for preheating simulations.As detailed in Refs.[87,88], the initial conditions are set by solving for the dynamics of the linearized equations of motion of the background and field fluctuations during inflation, again until two e-folds before inflation ends.At this time, the initial conditions for gauge field fluctuations hardly depart from the Bunch-Davies vacuum on the scales present in the simulation (as noted above).The dynamics of the homogeneous mode of the inflaton, however, are impacted by gauge-field backreaction; in both FLRW and BSSN simulations we set φ0 and φ′ 0 using the full numerical results.We consider values of α g between 8, the smallest for which preheating is efficient, and 14, which is roughly the largest value currently allowed by ∆N eff constraints on gravitational wave production derived in Refs.[86][87][88].
The larger the axial coupling, the wider the instability-thereby increasing the efficiency by which energy is transferred from the homogeneous mode to the gauge fields.At low couplings, α g ≲ 9, the gauge fields never fully dominate the energy budget of the Universe and we consider preheating to be incomplete.At couplings just above, with 9 ≲ α g ≲ 10, the majority of the energy contained in the homogeneous mode of the inflation is transferred to the gauge fields, but the process takes several oscillations.If 10 ≲ α g ≲ 12, the homogeneous mode of the field breaks down during the first oscillation.In this regime, there is a substantial amount of backreaction onto the modes of the inflaton, as well.For the highest couplings, α g ≳ 13, resonance is so strong that the homogeneous mode of the field never crosses zero before it decays.In these cases, backreaction onto the inflaton field is suppressed.At the highest coupling, the backreaction stalls the evolution of the inflaton on its potential, and inflation briefly restarts.These four regimes can be seen in fig. 1 for the FLRW case (left panels) which show both the energy contained in the gauge field as a function of time and the mean of the inflation.More details on this evolution can be found in [87,88].The right panel of fig. 1 shows the same quantities when implemented in the BSSN scheme.The inclusion of nonlinear gravity has little to no effect on the qualitative structure of preheating at the level of the background, for the wide range of parameters we have simulated. .Note that at early times (N ≲ 0) the gauge field energy densities differ between the two methods only due to the differing choice of cutoffs in initial conditions, as the dominant contribution is from vacuum modes (and unphysical).Nonetheless, the system exhibits large density contrasts, especially for larger couplings.As a first attempt, we can look for local gravitational effects using a linearized scheme.In conformal Newtonian gauge, where the scalar part of the metric is we can solve for the Newtonian potential Φ with the 00 component and the (divergence of) the 0i components of the Einstein equations.These respectively are and which we solve using the same method as described in [98].Figure 2 shows the statistics of the Newtonian potential calculated from our FLRW simulations in fig. 1.In this work, we calculate the Newtonian potential passively-with no feedback onto the evolution of the fields-to show that the Newtonian potential becomes too large, Φ > 0.25, to be able to treat gravity to linear order for the specific situations we consider.At this value, −g Newt = a 4 (1 − 4Φ) locally changes sign and linearized gravity, in conformal Newtonian gauge, necessarily breaks down.The left panel shows the absolute value of the minimum value of the Newtonian potential on each slice, the right panel shows the square root of the variance of the Newtonian potential across the grid.The solid horizontal lines denote the value Φ = 0.25, the point at which linearized gravity breaks down.Curves are not smooth only because these quantities are calculated at relatively infrequent intervals.

Density contrast and gravitational collapse
While the main features of the preheating story remain unchanged in the presence of nonlinear gravity, we now look to see if the additional nonlinear interactions provided by the gravitational sector affect the modes and scales that participate.Specifically, we search for hints that these interactions may lead to gravitational collapse.
The first place we look is at the power spectrum of density fluctuations.In fig. 3 we plot the dimensionless power spectrum of δ ≡ δρ/ρ, where the power spectrum is defined from the two point correlation function, When calculating ∆ 2 δ we directly Fourier Transform the ρ as calculated in eq.(2.33a) on the hypersurfaces of the simulation and ignore spatial dependence of the conformal factor.This measure is often cited as the litmus test for the need to incorporate nonlinear effects (see, for example, Ref. [114]) as well as a useful measure to determine whether primordial black holes are produced.This statistic was first used as a way to diagnose potential primordial black hole formation in Ref. [115] and continues to be used, see, for example, Refs.[31,116].
Figure 3 shows that, in many cases, dimensionless power spectra are of order ∆ 2 δ (k) ∼ 1 for large α g .This is reinforced by the fact that perturbation theory breaks down at these times; more precisely, in cases where δ(k) approaches unity, the linearized Einstein's equations are violated as we saw in fig. 2 [106].varying by panel, comparing results from simulations implementing full general relativity via the BSSN scheme (solid blue lines) and FLRW simulations (dashed red).All results are evaluated two e-folds after the end of inflation.The smaller, higher-frequency peak in the FLRW simulations are a numerical artifact that corresponds to the scale of the cutoff of the initial conditions.
As one would expect, the simulations with nonlinear gravity show more numerical effects at high wavenumber; nonetheless, the physical parts of the power spectra are very consistent with the FLRW counterparts.The only exception to this is, perhaps, the very highest values of the self-coupling where there seems to be somewhat higher power at lower scales and somewhat lower power at intermediate scales.We speculate that this is due to the fact that local gravity can lead to some clumping in the nonlinear simulations.
In order for clumps to actually collapse into black holes, overdensities must overcome their own internal pressure.In linear theory, the Jeans scale is used to determine whether the scales of interest can collapse.Modes with wavelength longer than the Jeans length undergo gravitational collapse, while those with shorter wavelength are pressure supported.This scale defines to the (physical) Jeans scale, where c 2 s = δP/δρ is the sound speed.The comoving Jeans scale is k J = ak phys J .Since the gauge fields are the main component of the universe, we estimate c 2 s = 1/3 to approximate the sound speed for a radiation fluid.At the end of our simulations-where the power spectra in fig. 3 are evaluated-the Jeans scale k J is much smaller than the peak frequency of the dimensionless power spectrum of the density fluctuations.This indicates that the radiation pressure of the gauge fields is playing an important role in the evolution of these structures.To use α g = 13 and α g = 14 as examples, we can track the location peak of the dimensionless power spectrum throughout simulation which we can compare to the Jeans scale and the Hubble scale as can be seen in fig. 4. The differences in the evolution of the comoving scales is due to the fact that, for the largest coupling, the gauge field backreaction briefly restarts inflation.Fig. 4 shows us that Jeans-scale modes are excited near the end of inflation; however, these modes are still small.As the density contrast grows, the peak of the power spectrum moves to slightly larger comoving wavenumber, while the comoving Jeans scale shrinks.By the time we reach two e-folds after the end of inflation, the peak mode of the power spectrum is approximately a factor of five larger than the Jeans scale.
A number of different tools are used to identify the existence of black holes when studying critical collapse in numerical relativity [101].In 1 + log slicing, where (∂ t + β i ∂ i )α = −2αK, the vanishing of the lapse α can sometimes be used as an indicator that black holes have formed [117][118][119][120][121][122][123].The same is true for our slicing condition, eq.(2.18), where the only difference is the background clock which is constructed to mimic a conformal-time FLRW solution in the homogeneous limit.In the weak gravity limit, we can compute the inhomogeneity of α to measure how much our slices vary from the FLRW limit, with smaller values of α corresponding to deeper gravitational wells.Figure 5 shows 2-dimensional slices of three different simulations evaluated two e-foldings after the end of inflation.In these figures we can see spatial variation in the density contrast-which is often large, O (10).The lapse, however, does not deviate much more than O 10 −1 .We can extend this to a statistical test by calculating the variance of the density contrast δ and the lapse α throughout the simulation.As we can see in fig.6, the variance of the lapse increases with the size of the axial coupling α g .However, the deviation from the average is never larger than O(0.3).
As a final metric, we look to determine whether regions in our simulations have passed thresholds that one might consider sufficient to produce black holes.We begin by examining just one of the more dramatic runs and calculating the compactness [124][125][126][127] generalized to an expanding spacetime [128] to calculate whether the overdensities on the final slice indicate that black holes might form.In eq.(3.8) δM is the instantaneous over-mass enclosed in some (proper) radius, R.That is, we subtract off the background density when computing the mass.Note that we use a different definition of compactness than that in Ref. [129], where the authors measure the  The comoving Jeans scale compared to the dimensionless power spectrum.The left panels show the range of modes in the peak-most power-bin of the dimensionless power spectrum (blue, shaded region) compared to the Jeans scale (red), the Hubble scale (green), and the longest resolvable mode in the box (black).We only calculate the location of the peak of the dimensionless power spectrum and the Jeans length after the end of inflation when the power spectrum raises above its initialized shape and we can approximate the Universe as radiation dominated.The right panels show the dimensionless power spectrum at the end of inflation (blue, dashed lines) and two e-folds after the end of inflation (blue, solid lines) with vertical lines showing the Jeans corresponding jeans scale at the end of inflation (red, dashed lines) and two e-folds after the end of inflation (red, solid lines).The top panels correspond to α g = 13 and the bottom panels correspond to α g = 14.
total integrated ρ within a (not necessarily spherical) region where ρ/ρ > 5%, although both definitions reduce to a statement of the hoop conjecture [102] about a spherically symmetric overdensity and when ρ/ρ ≫ 1.
To test this, we look at the final slice for the largest coupling we test, α g = 14, and compute the compactness of the largest over-density in the box.Centering coordinates about the location of maximum δ(x), we can calculate the enclosed over-mass and, following [129], we approximate the radius of this overdensity from In both eq.(3.9) and eq.(3.10) the final sum is over all points inside a given radius r. Figure 7 shows the compactness as a function of distance away from the center-a measure that is maximized around C ∼ 10 −1 , yet is still smaller than the critical value of 1/2.In fig.7 we also plot the FLRW approximation to the compactness where δM FLRW = a 3 ρ∆x 3 i (δ i − 1) and R FLRW = a∆x where the scale factor is calculated in the FLRW limit, a 2 = ē 4ϕ , see appendix B.

Gravitational waves
One of the most striking features of gauge preheating is the strength of the production of gravitational radiation.During and after gauge preheating, gravitational wave production can be so strong that the resulting radiation density stored in gravitational waves leads to a shift in the effective number of relativistic species N eff that is large enough to be measured or constrained by current and future CMB experiments [86][87][88].In some models, the shift is large enough that the allowed values of the axial coupling are already constrained by existing CMB measurements.These results were obtained using numerical simulations that ignore Figure 6: Dynamics of energy transfer and the metric lapse during preheating.Plotted is the fraction of energy in the gauge fields, ρA /ρ (solid blue, with scale given by the left of each panel) and the root-mean-squared and maximum deviation of the lapse from its average over space (solid green and transparent, dashed green, respectively, with scale given by the right of each panel).Each panel corresponds to a simulation with axial coupling strength indicated on the plot.the backreaction of these large metric fluctuations.In this section, we explore the robustness of predictions of gravitational wave production during gauge preheating in the regime where gravity is nonlinear.
Gravitational waves are the traceless part of the spatial metric, and are contained in γij .In FLRW simulations, the gravitational wave spectrum is computed passively-see, e.g.[24], by calculating the transverse-traceless parts of the matter stress tensor, where P is the projection operator, This sources the transverse-traceless part of the metric via Einstein's equations We can then calculate the power in gravitational waves from the gravitational wave stressenergy tensor [102] T where The energy density in gravitational waves at the time of emission is then given by Since the BSSN system evolves the entire metric by solving the full nonlinear Einstein equations, the gravitational wave spectrum can be directly extracted by projecting out the transverse-traceless part of the extrinsic curvature, h ′ ij ≈ −2 Ãij [96].The BSSN analog to eq. (3.16) is then which is analogous to the quantities used in Refs.[97,130].
The spectral energy density of gravitational waves today can be calculated using the standard transfer functions [21,22], assuming that the signal is evaluated during a radiationdominated period and remains radiation dominated until matter-radiation equality.The spectral energy density today is given by and the frequency by g ⋆S (a r ) g ⋆S (a 0 ) In the preceding expressions g ⋆ is the number of ultra-relativistic degrees of freedom evaluated at reheating (a r ) or today (a 0 ).For simplicity, we take g ⋆ (a r )/g ⋆ (a 0 ) ≈ 100.
In fig.8 we show the gravitational wave spectra from simulations of gauge preheating in both full nonlinear gravity, as well as in rigid FLRW spacetime.Despite the presence of strong nonlinearities, we observe remarkably consistent outputs across the range of couplings we considered.These results suggest that neglecting nonlinear gravity has at most an O(1) effect on the resulting gravitational wave spectrum, even in regions where the metric is highly nonlinear.

Conclusions
In this paper we have extended studies of gauge preheating after pseudoscalar-driven inflation to include the effects of nonlinear gravitation.To facilitate this study, we implemented numerical relativity using the BSSN formalism in our simulation software.
The evolution of the energy density and fields was nearly indistinguishable from simulations done in a FLRW spacetime.Including the effects of nonlinear gravity made no qualitative difference in the value or evolution of the averaged background quantities, such as the expansion rate, the average value of the scalar field, or the energy densities of the scalar and gauge fields.
In our FLRW simulations we observed the emergence of regions with very large fractional overdensities δρ/ρ, which routinely exceeded unity.The existence of these regions leads to the failure of linearized gravity, requiring software that properly treats evolution of the metric into the nonlinear regime.In our BSSN simulations, we observe power being shifted from large to small scales due to the gravitational interactions.
In our highest-coupling BSSN runs, we found regions with fractional overdensities as large as δρ/ρ ∼ 30.However, despite the development of these large density contrasts we found no evidence for the formation of black holes such as the presence of horizons or the vanishing of the lapse function.Closer examination of these simulations revealed that the scales where the density power spectrum peaks are within the Jeans length.This suggests that pressure of the radiation plays an important role in the evolution and stability or instability of these very overdense regions.We also computed the compactness of these regions, and investigate the dynamics of the large overdense regions to study their fate.In particular, it would be interesting to investigate whether the large overdensities subsequently undergo gravitational collapse, or whether they simply decay.Our simulations have so far focused on the large density fluctuations that are generated on sub-horizon scales during reheating.These are necessarily inside the Jeans length at production.It would be interesting to study the horizon reentry of large curvature perturbations that are produced near the end of inflation in these models to see if their collapse to black holes can be verified in full general relativity.

B.1 Metric perturbations
We define the perturbed metric where h µν is a small perturbation with scalar-vector-tensor decomposition Here C i and G i are transverse vectors-namely, We expand the BSSN variables to linear order about a homogeneous and isotropic background spacetime as α = α 0 + δα (B.3a) Expanding the BSSN metric in these variables, we identify e 4ϕ 0 = α 0 (τ ) (B.9) We can linearize the equation of motion for γij to obtain remembering that ∂ t in the BSSN formalism corresponds to conformal time for our particular gauge choices.

B.2 Initial conditions
The linearized tensor parts of a ij and δγ ij , eq. (B.3), comprise genuine, propagating degrees of freedom, so their initial condition is specified independent of constraints.We therefore choose D ij = 0 and D ′ ij = 0 from eq. (B.2) initially.The remaining degrees of freedom are either gauge choices or set by constraints.We take β i = δβ i = 0 as an initial condition.Correspondingly, we set F = 0 and G i = 0 from eq. (B.2), amounting to the fixing of two vector and one scalar gauge modes.The initial conditions are simplest to solve for when taking B = 0 as the final gauge choice.
We define the short-hand Finally, Gauss's law expands to and, since we take ∂ i ∂ t A i = 0 initially, we satisfy this constraint by setting

C Robustness checks
For the six BSSN simulations presented here, we calculate the violation of the Hamiltonian and momentum constraints, eqs.(2.16) and (2.17), and of Gauss' Law, eq.(2.32), to ensure that we stay on the solution surface of the problem.In fig. 9 we plot these constraints vs N ≡ ln a to show that they are bounded throughout the simulation.For the results presented in this work, we take a conservative cutoff of N = 2 as a final time.This allows us to consider all of the BSSN simulations in regime where the constraints are small.In addition, the results in section 3 are well reproduced by simulations with 256 3 rather than 384 3 grid points; the latter is the largest grid that is computationally tractable with our resources.We choose this largest grid for our main results in order to maximize the range in scales over which the BSSN results may be reliably compared to the FLRW results, e.g., in fig. 3.
For FLRW simulations, the first Friedmann equation, H 2 = ρ/3M 2 Pl , serves as a constraint that measures the analog of energy conservation in flat space; the simulations satisfy it to one part in 10 4 or better.The gauge constraints are likewise satisfied to one part in 10 3 ÷ 10 2 (see Refs. [40,41,86,87] for further discussion of the robustness of these methods).The FLRW results are all consistent with those of lower-resolution simulations presented in Refs.[86][87][88].

Figure 1 :
Figure1: Energy fraction in the gauge fields (top panels) and homogeneous component of the inflaton (bottom panels) for FLRW simulations (left panels) and BSSN simulations (right panels).Note that at early times (N ≲ 0) the gauge field energy densities differ between the two methods only due to the differing choice of cutoffs in initial conditions, as the dominant contribution is from vacuum modes (and unphysical).

2 Figure 2 :
Figure 2: Statistics of the Newtonian potential calculated passively from FLRW simulations.The left panel shows the absolute value of the minimum value of the Newtonian potential on each slice, the right panel shows the square root of the variance of the Newtonian potential across the grid.The solid horizontal lines denote the value Φ = 0.25, the point at which linearized gravity breaks down.Curves are not smooth only because these quantities are calculated at relatively infrequent intervals.

Figure 3 :
Figure3: Spectra of density fluctuations δ = δρ/ρ in simulations with axial coupling α varying by panel, comparing results from simulations implementing full general relativity via the BSSN scheme (solid blue lines) and FLRW simulations (dashed red).All results are evaluated two e-folds after the end of inflation.The smaller, higher-frequency peak in the FLRW simulations are a numerical artifact that corresponds to the scale of the cutoff of the initial conditions.

Figure 4 :
Figure 4:The comoving Jeans scale compared to the dimensionless power spectrum.The left panels show the range of modes in the peak-most power-bin of the dimensionless power spectrum (blue, shaded region) compared to the Jeans scale (red), the Hubble scale (green), and the longest resolvable mode in the box (black).We only calculate the location of the peak of the dimensionless power spectrum and the Jeans length after the end of inflation when the power spectrum raises above its initialized shape and we can approximate the Universe as radiation dominated.The right panels show the dimensionless power spectrum at the end of inflation (blue, dashed lines) and two e-folds after the end of inflation (blue, solid lines) with vertical lines showing the Jeans corresponding jeans scale at the end of inflation (red, dashed lines) and two e-folds after the end of inflation (red, solid lines).The top panels correspond to α g = 13 and the bottom panels correspond to α g = 14.

Figure 5 :
Figure 5: Two-dimensional spatial slices of the density contrast, δ, (top panels) and lapse (normalized to its average over all space; bottom panels), α/ ⟨α⟩.Columns display results for simulations with axial coupling α g = 10, 12, and 14 from left to right.All results are evaluated two e-folds after the end of inflation.

Figure 7 :
Figure7: The left panel shows the compactness C as a function of distance from the center of the largest overdensity present in the final slice of the simulation with α g = 14 (black, solid).We also plot the FLRW approximation to the compactness (red, dashed).The right panel shows the z = constant slice through the center of the overdensity.

3 e 2 Pl 1 ∇ 2
ik•x −k 2 d 3 y e −ik•y f (y) (B.16)to denote the solution to Poisson's equation, ∇ 2 y(x) = f (x).The constraint on the vectors may be solved straightforwardly via1 2α 0 ∂ t C m = − 1 M S m .(B.17)This sets the vector contribution a mn as −∂ (n ∂ t C m) /2α 0 .For the scalar degrees of freedom, with B = 0 one may directly solve the scalar part of the momentum constraint as 2 2 2