Lattice Simulations of Inflation

The scalar field theory of cosmological inflation constitutes nowadays one of the preferred scenarios for the physics of the early universe. In this paper we aim at studying the inflationary universe making use of a numerical lattice simulation. Various lattice codes have been written in the last decades and have been extensively used for understating the reheating phase of the universe, but they have never been used to study the inflationary phase itself far from the end of inflation (i.e. about 50 e-folds before the end of inflation). In this paper we use a lattice simulation to reproduce the well-known results of some simple models of single-field inflation, particularly for the scalar field perturbation. The main model that we consider is the standard slow-roll inflation with an harmonic potential for the inflaton field. We explore the technical aspects that need to be accounted for in order to reproduce with precision the nearly scale invariant power spectrum of inflaton perturbations. We also consider the case of a step potential, and show that the simulation is able to correctly reproduce the oscillatory features in the power spectrum of this model. Even if a lattice simulation is not needed in these cases, that are well within the regime of validity of linear perturbation theory, this sets the basis to future work on using lattice simulations to study more complicated models of inflation.


Introduction
Inflation is the accelerated expansion of the early universe originally introduced as a solution to the horizon and flatness problems of Cosmology [1][2][3][4][5]. In the standard picture this early accelerated expansion is driven by a scalar field, the inflaton, which acts as a source with a negative pressure. The standard theory of inflation predicts scalar and tensor perturbations in the early universe as quantum vacuum fluctuations [6][7][8][9][10][11]. While the tensor ones (i.e. gravitational waves) remain undetected, the scalar perturbations are observed as temperature fluctuations in the CMB radiation [12]. This outstanding prediction makes inflation very appealing as a scenario for the early universe and opens up to many theoretical and observational challenges.
In this paper we aim at studying the inflationary phase of the universe with a lattice simulation. Many lattice simulations have been developed to study the reheating phase at the end of inflation [13][14][15][16][17][18][19][20][21][22], where the inflaton decays transferring its energy to the thermal hot bath of the early universe. In the case of reheating lattice simulations are needed because of the non-linear physics involved.
We will use a lattice simulation to predict the well-known results of some simple single-field models of inflation, particularly for the scalar field perturbation. In these standard models there is no need to use a lattice simulation to study the dynamics of the inflaton, which is well understood with linear perturbation theory. However, this first step is needed to understand the technical aspects involved in simulating the inflationary universe on a discrete lattice. This sets the basis to future work on using lattice simulations to study more complicated inflationary models. Indeed, as we will discuss in section 6, there are many examples of non-standard inflationary models that involve non-linear physics, and where a lattice simulation might be useful.
We will show which aspects need to be taken into account to reproduce the nearly-scale invariant power spectrum of inflaton perturbations for the simplest slow-roll single-field inflationary model. Such technical aspects include, for example, taking into account the modified dispersion relation induced by the finite grid spacing on the modes propagating on the lattice. Our simulation is written in C++ and it is OpenMP parallelized. It is inspired on LATTICEEASY [15], a lattice code for evolving fields in expanding space-times that has been developed to study the reheating phase of the universe.
The paper is organized as follows. In section 2 we start with a quick review of inflation to set the notation and recap the equations needed in the rest of the work. In section 3 we introduce the lattice simulation and we write explicitly the discrete version of the equations of motion. In section 4 we discuss the difference between continuous and discrete dynamics in Fourier space, and how we need to take it into account to reproduce the results from the linear theory. In section section 5 we show the results of the simulation on two simple examples of single-field inflation. Finally, in section 6 we summarize the conclusions and implications of this work.

Single field Inflation
We start with a quick recap of inflation in order to present the notation and the equations that will be used in the rest of this work. A more detailed introduction to inflation can be found, for example, in [23]. In the standard inflationary picture the action of the early universe is the following: where φ is the inflaton and V (φ) its potential. From now on we set M Pl = 1. In the usual perturbation theory approach, both the scalar field and the metric g µν are split into a homogeneous background quantity plus a space dependent perturbation. As we will see in section 3, this is not true for the lattice simulation, where the inflaton is evolved as a whole without any splitting and the metric is left unperturbed. In this section, however, we still split the inflaton in a background quantity plus perturbations φ( x, t) =φ(t) + ϕ( x, t) and discuss them separately. This is only done to illustrate the well established results of linear perturbation theory.
Here and in the rest of the paper, the background metric is assumed to be the spatially-flat Friedman-Lemaitre Robinson-Walker metric, with line element: where we set c = 1.

Background evolution
The expansion of the universe is governed by the Friedman equations: where ρ and p are the background energy density and pressure of the scalar field, and the dots represent derivatives in cosmic time t. At the background level they are written as ρ =φ 2 /2 + V and p =φ 2 /2 − V . From the second Friedman equation we can see that ifφ 2 /2 V the universe undergoes an accelerated expansion. This is usually called slow roll condition. The rate of acceleration is described by the following parameter: The universe expands in an accelerated way if < 1. Typically, this parameter is much smaller than 1 during slow-roll inflation 1, and approaches 1 at the end of inflation. Varying the action (2.1) with respect to the field yields to the Klein-Gordon equation: where H = a −1 ∂ τ a. This equation will determine the motion of the background value of the inflaton on the inflationary trajectory.

Quantum perturbations
As mentioned in the introduction, the quantum fluctuations of the inflaton are responsible for the fluctuations of the CMB. In this framework the inflaton is promoted to a quantum operator: where a and a † are the creation and annihilation operators satisfying [a k , a † k ] = δ( k − k ). Varying the action (2.1) at second order in perturbation leads to the Mukhanov-Sasaki equation for field fluctuations: where v( k, τ ) = aϕ( k, τ ) is the Mukhanov variable 1 and k = | k|. This equation is an harmonic oscillator with a time dependent mass term: The time dependence of the mass introduces a formal ambiguity in choosing the physical vacuum state of the theory. This ambiguity is solved by choosing the so called Bunch-Davies vacuum, which is the ground state of the inflaton in Minkowski space corresponding to the asymptotic past of the theory. This vacuum is associated to the mode functions: where m is the mass of the inflaton. These mode functions serve as the initial conditions for eq. (2.7), which can be solved numerically together with the background equations of section 2.1 to determine the evolution of field fluctuations at linear order. In the case of a free massive scalar field of mass m, the effective mass appearing in the Mukhanov-Sasaki equation simplifies to: In this case we can write an analytical solution to the Mukhanov-Sasaki equation with the initial conditions given by eq. (2.9) in the asymptotic past τ = −∞. The solution can be written as: where H (1) ν is the modified Hankel function of the first kind. We can use this solution to write the power spectrum of inflaton perturbations as P φ (k) = |ϕ( k, τ )| 2 . This gives a theoretical prediction for the dimensionless power spectrum of inflaton perturbation: (2.12)

Lattice simulation
The idea behind a lattice simulation is simple and it consists of simulating the dynamics of continuum fields on a finite cubic lattice. The lattice is defined as a collection of N 3 points separated by comoving lattice spacing dx = L/N , where L is the comoving physical size of the box. To any given field in continuous space, we associate N 3 values to each point of the cubic lattice. For the inflaton φ this corresponds to: To enlighten the notation, from now on we will write i ≡ i 1 , i 2 , i 3 for discrete vector indices. Here and for the rest of this work, we use periodic boundary conditions on the lattice. Contrarily to what is done in perturbation theory, in the simulation we do not split in background and perturbation quantities. Indeed, the inflaton on the lattice is evolved altogether using the Euler-Lagrange equations in real space: where L[φ] is a discretized version of the Laplacian operator. We will mainly consider the following second order stencil for the Laplacian operator: In appendix C we also consider other stencils for the Laplacian, and show how the choice of the Laplacian affects the dynamics of the simulation. The space-time is evolved neglecting the perturbations of the metric, which is assumed to be the unperturbed FLRW metric (2.2). Indeed, the scale factor appearing in eq. (3.2) is evolved using the second Friedman equation: In principle, metric perturbations can be included in a full numerical GR computation, as we see for example in [24,25]. We choose to neglect them because, as it is well known, their coupling with φ is slow-roll suppressed 2 . The equations of motion of the system (3.2)(3.4) determine the evolution of the inflationary universe and of the inflaton field on the lattice. We solve them numerically with a second order staggered-leapfrog integrator inherited from LATTICEEASY [15]. As wee can see, the equations of motion for the evolution of the system are purely classical, even if the inflaton is a quantum field. The quantum behavior of the inflaton will be captured by the initial conditions that, as explained in section 4.3, are randomly generated over the lattice. In this picture the uncertainty associated to the quantum nature of the field will be replaced by a statistical uncertainty over the different random realizations of the lattice simulation. This semi-classical approximation is common and it constitutes the working assumption of most of the lattice simulations in the context of inflation.

Discrete dynamics in Fourier space
In this paper we are interested in reproducing the results from the well known linear theory of inflaton perturbations. For this reason, even if the lattice simulation evolves fields in real space, we need to explicitly define how we compute quantities in Fourier space. This is crucial in order to implement the initial conditions, which are given in Fourier space, and to correctly interpret the outputs of the code. As we will see in section 4.2, this will be non-trivial due to the modified dispersion relation induced by the discretization of space on the modes propagating on the lattice.

Discrete Fourier transform and reciprocal lattice
We start with the definition of lattice modes. For a given discrete field f i living on the lattice we define its Discrete Fourier Transform (DFT) in the following way 3 : In our notationf is the Fourier transform of f for any given field. The Fourier fields live on the reciprocal lattice where we associate to each point the following comoving momentum: With this definition, we can write the inverse DFT (iDFT) as follows:

The modified dispersion relation
In this section we describe the effect of the discretization on the propagation of modes on the lattice. In continuous space, the Fourier transform (FT) of the Laplacian operator for differential equations is quite simple and reads: As it is well known, this relation gets modified on the lattice [26], where we transform the field with the Discrete Fourier Transform (DFT) (as pointed out for example in [17]). It can be easily derived from eq. (3.3) and eq. (4.3) that: where we introduced the effective modes k eff as: Contrarily to what happen in the continuous case, this relation differs significantly from the value of the k−modes of the reciprocal lattice of eq. (4.2) k eff,l1,l2,l3 = k lat,l1,l2,l3 . Indeed, k lat and k eff are only equal in the limit l 1 , l 2 , l 3 N . Note that the expression for k eff of eq. (4.6) depends on the definition of the lattice Laplacian of eq. (3.3). A different choice of the numerical stencil for the Laplacian would lead to a different expression for k eff .
This effect is quite general and we can interpret it as a modified dispersion relation induced by the discrete spacing on the modes propagating on the lattice. Indeed, if we look at the equation for a free, massless scalar field on the lattice we can see that modes will propagate with energy ω(k lat ) = k eff = k lat , which is different from the usual ω(k) = k of continuous space.
In the case of inflation, we can see the effect of the modified dispersion relation by looking at the lattice version of the Mukhanov-Sasaki equation, which is obtained from eq. (2.7) by replacing k with k eff : where v l = aφ l in analogy to the continuous case. We will use this equation to predict the evolution of perturbations in the simulation. In a similar way to section 2, we can write an exact solution to this equation in the case of a free massive scalar field: From this solution we can write the lattice dimension-less power spectrum as which differs from eq. (2.12) by the presence of k eff instead of k inside the Hankel function. We will discuss in detail the consequences of this in section 5. Notice that the modified dispersion relation will also influence the ground state in Minkowski space, i.e. the Bunch-Davies initial conditions, by changing the expression for the ω k appearing in eq. (2.9). We will take this into account when discussing the initial conditions in section 4.3. From the discussion above, and in particular from eq. (4.8) and eq. (4.9), it follows that the dynamics of the classical Fourier modes on the lattice is equivalent to the one of the continuous (quantum) modes ϕ( k, τ ) if we interpret the effective modes k eff as the physical modes actually probed by the lattice simulation. It follows that, when computing the dimensionless power spectrum P φ from the simulation, we should multiply the dimension-full power spectrum |φ l | 2 by k 3 eff instead of k 3 lat if we want to match the continuous result 4 . In other words, the non-scaled dimension-full power spectrum of the simulation |φ l | 2 will not scale as ∼ k −3 lat , as one would naively expect, but it will scale as ∼ k −3 eff . We will discuss this in more detail when looking at the results of the simulation in section 5.
Moreover, this will have consequences on the effective spacial resolution of the simulation. Instead of probing modes up to 5 k lat,max = 2π √ 3N Nyquist /L, where N Nyquist = N/2, it will probe physical modes up to: This means that the effective range of physical modes evolved by the simulation will be reduced by a factor of 2/π 0.64. As already mentioned, a different definition of the Laplacian would lead to a different expression for k eff and to a different value of k eff,max . In appendix C we show the comparison between the k eff associated to different stencils and we discuss the consequences on the dynamics of the simulation.

Initial field fluctuations on the lattice
We now explain how we generate the initial conditions on the lattice taking into account the modified dispersion relation discussed in the last section. The first step is defining the discrete version of eq. (2.6): lφ j e +i 2π N i· l , (4.12) 4 As we will see, this prescription is useful to compare the power spectrum of the simulation with the continuous theory. However, the reader should keep in mind that only the quantity P (lat) φ (k lat ) of eq. (4.10) is related to the actual variance of the field over the N 3 points of the lattice. 5 More details about k lat,max can be found in appendix B.
where in the last equality we show the comparison with our definition of Fourier modesφ j . We have also introduced the discrete creation and annihilation operators: Here u l are the discrete mode functions, whose exact expression is: (4.14) At the beginning of the simulation the comoving size of the lattice will be smaller than the horizon L aH. For this reason, eq. (4.14) will practically reduce to the Bunch-Davies vacuum of eq. (2.9) for most of the modes. The extra normalization factor L 3/2 is introduced to correct for the finite volume of space 6 . Note that in LATTICEEASY, and in most of the lattice simulations in the context of reheating, the initial fluctuations are usually generated using k lat instead of k eff in eq. (4.14). Once the discrete mode functions u l are determined, field fluctuations are generated in Fourier space using eq. (4.12), initiating every given modeφ l as a Gaussian random number with variance |u l | 2 . We review this procedure in appendix A.

Results of the simulation
We now proceed showing the results of the simulation. The main model that we consider is a standard slow-roll inflationary potential for the inflaton V (φ) = 1 2 m 2 φ 2 . We will focus on this model to show the differences between continuous and discrete dynamics described in section 4. As a further example, we consider also a similar potential but with a step added on top of it. All numerical values in this section are given in Planck units M Pl = 1. More details on how the outputs are computed in the numerical code can be found in appendix B, where we also discuss energy conservation in B.1. In appendix A we discuss in detail how we generate the initial conditions on the lattice.

Standard slow roll potential
We first discuss the result for a simple harmonic potential for the inflaton where m = 0.51 · 10 −5 . This value is chosen to roughly match the COBE normalization of the power spectrum of curvature perturbation P ζ 25 · 10 −10 . The initial average value of the inflaton is chosen to beφ in = 14.5. Its velocity is determined solving numerically the background Klein-Gordon equation (2.5) and it is given byφ in = −0.8152m. With these values the universe is in the middle of the inflationary phase, and there are N e 53 e-folds 7 left before the end of inflation. The system is evolved until a = 10 3 (N e 6.9) which means that at the end of the simulation we will still be in the inflationary phase.
In the lattice simulation the inflaton is evolved as a whole, meaning that there is no splitting between background and perturbation quantities. However, we still perform the splitting at the output level in order to compare our results with linear perturbation theory. In fig. 1 we show the evolution of the background value of the inflatonφ and its velocityφ as functions of the number of e-folds N e . These quantities are computed from the simulation as averages over the N 3 points of the lattice. In the same plot, we also show the evolution of H and ε . The initial bump in the field velocityφ is a numerical effect due to the initial stabilization on the inflationary trajectory. From these plots we clearly see that we are in the middle of the inflationary phase, being ε 1 andφ constant We now come to the dynamics of field fluctuations and to the importance of the modified dispersion relation discussed in section 4. We show results from a run of the code with L = 1.4/m and N 3 = 128 3 . This translates to: where H in = H in is the initial value of the Hubble rate. The modes are almost all sub-horizon at the beginning of the simulation. We evolve the system until a = 10 3 , which means that the modes will be all super-horizon at the end of the simulation.
In the left panel of fig. 2 we show the dimensionless power spectrum of the inflaton P φ at the end of the simulation, plotted against lattice modes k lat of eq. (4.2). We compare the power spectrum computed from the simulation with the theoretical prediction for discrete dynamics of eq. (4.10), which is shown as a green line. From this plot we can see that the discrete power spectrum is quite different from the almost scale-invariant power spectrum of the continuous theory, given by eq. (2.12) and depicted as a blue line in the plot. This is a manifestation of the different dynamics in Fourier space between the discrete and the continuous case and it is a manifestation of the modified dispersion relation. However, as we have discussed in section 4.2, continuous and discrete dynamics are equivalent if we interpret k eff of eq. (4.6) instead of k lat as the physical modes probed by the lattice simulation. In fig. 3 we can see a comparison between k eff and lattice modes k lat . The departure from the diagonal is a manifestation of the modified dispersion relation discussed in section 4.2. The one dimensional quantities k eff and k lat are obtained averaging eqs. (4.2) and (4.6) over spherical bins on the lattice 8 .
In the right panel of fig. 2 we show the power spectrum from the simulation computed interpreting k eff as physical modes and we compare it to the theoretical prediction of eq. (2.12). From these plots P φ k lat /aH k eff /aH  we see that only in this case we are able to reproduce with precision the nearly scale-invariant spectrum of single-field inflation. Note that interpreting k eff as the physical modes will also reduce the resolution in Fourier space, that is computed from eq. (4.11) and it is given by k eff,max 53.50H in . In appendix C we also consider results of simulations with different stencils for the discrete Laplacian and compare the corresponding effective momenta.
In fig. 4 we show the evolution of the power spectrum during the simulation, plotted at different times as a function of physical modes and going from the early-time Bunch-Davies state to the final scale-invariant state.

Potential with a step
As a further example, in this section we show the results of the code for a model with potential: This potential is analogous to the harmonic potential of the last section but with a step localized at φ step . This model has been studied in [27], where the authors show that the presence of the step causes oscillations in the power spectrum of scalar perturbations. Here we show results for the same parameters of the last section. The only difference here is that we use L = 0.6/m as comoving size of the box, that corresponds to k eff,max 124.84H in . Moreover, we have three extra parameters s, d and φ step . We choose φ step = 14.35, and we run the simulation with different values of s and d.
In fig. 5 we show the background quantities in the case s = 0.01, d = 0.005. We can see here that the step of the potential causes a bump in all the background quantities, but without changing significantly the slow-roll dynamics of the inflaton. Indeed, ε is still smaller than 1 during the simulation and the departure of the inflaton from the slow-roll trajectory is small.
In fig. 6 we show the evolution of the power spectrum during the simulation for s = 0.01 and d = 0.005. Here we can clearly see that the presence of the step introduces oscillations in the power spectrum.
In fig. 7 we show the final power spectrum of a simulation run with s = 0.001, d = 0.005 and we compare it with the result obtained solving the Mukhanov-Sasaki equation (2.7) with a numerical integrator, which serves as a theoretical prediction. For this simulation we increased the number of lattice points to N 3 = 256 3 and the box size to L = 1.2/m in order to improve the spatial resolution. In the right panel of this figure we show the result obtained by interpreting k eff the physical modes, while on the left panel we show the lattice power spectrum. From the right plot we can see that the matching between the theoretical prediction and the lattice simulation is not perfect, in particular for the largest modes of the simulation. However, the lattice code is able to correctly reproduce the oscillations, that have the same frequency and a similar amplitude compared to the theoretical prediction. In this example we can again see that interpreting k eff as the physical modes is important in order to get a sensible result.

Conclusions and outlook
In this work we studied some well known models of single-field inflation with a lattice simulation. After introducing the numerical code, we discussed the effects induced by the discrete lattice that need to be understood in order to reproduce the well known results of single-field inflation for the scalar field perturbation. We have seen that an important role in reproducing the results of the continuous theory is played by the identification of the physical modes that are probed by the lattice simulation. In order to do so, we studied how the finite grid spacing influences the dispersion relation of modes propagating on the lattice. Taking into account the modified dispersion relation was necessary in order to correctly interpret the power spectra computed from the simulation and to properly initialize the scalar field on the lattice. This allowed us to match the results from the linear theory in two cases. The first was a simple harmonic potential for the inflaton, while the second was a similar potential but with a step added on top of it.
As we already stated, a non-linear lattice simulation is not needed to understand the physics of these simple models of inflationary potential, that lies well within the regime of validity of linear perturbation theory. However, there are many non-standard models of inflation that involve non-linear processes and that have been studied in the literature because of their interesting phenomenology. Examples are the geometrical destabilization of inflation [28], or inflationary models with abelian [29][30][31][32][33] and non-abelian [34][35][36][37][38] gauge fields. These models potentially manifest strong backreaction effects, that invalidate the use of perturbation theory [28,[39][40][41][42]. For this reason, lattice simulations might turn out to be an important tool in order to compute observables from these kind of models, such as the 2-and 3-point functions of scalar (and tensor) perturbations. This work constitutes a first step in this direction. Indeed, even if some examples of lattice study during inflation already exist in the literature (see [43] for example), this is the first time in which a lattice simulation is used to recover the 2-point function of the inflaton in single-field models with such precision. This is a necessary result if one wants to use a lattice simulation to compute observables from more complicated models. Moreover, a precise computation of the 2-point function is a necessary step in developing a code that is able to compute 3-point statistics during inflation, which might be important in many applications.

A.1 Background quantities
The initial background values of the inflatonφ in and its velocityφ in are set by requiring that the universe is in the middle of the inflationary phase. This will depend on the inflaton potential, that we will write explicitly in section 5. The scale factor a is simply set to 1 at the beginning of the simulation, while its derivative in program time a is computed from the first Friedman equation (2.3) using only the background energy density and pressure of the field and neglecting gradient contributions: Then, after the field fluctuations are generated, the value of a in is updated to include the gradient term, computed as a lattice average of (∂ i φ) 2 /2. Note that quantum vacuum sub-horizon fluctuations should not contribute to the Friedmann equations. However, we still include them in generating the initial value of H, and this is a consequence of our semi-classical approximation. This important in order to evolve the discrete system in a consistent way. Indeed, neglecting gradient contributions will result in an effective residual curvature in the second Friedmann equation, that we use to evolve the scale factor during the simulation. Moreover, including the gradient term in a consistent way allows us to check that energy is conserved in the discrete system, i.e. to ensure that there are no numerical errors propagating on the lattice during the simulation. More about the energy conservation check can be found in appendix B.1.

A.2 Initial fluctuations of the field
The field fluctuations are generated starting from the expression of the mode functions of eq. (4.14). The extra factor of L 3/2 is a common normalization [15,16] introduced to correct for the finite volume of space. In order to understand this, take the two point function of the field: We can clearly see that this scales as L −3 due to the presence of the finite-volume delta function δ L . If we want the quantity φ 2 , and the two-point functions in general, to be independent of the physical size of the lattice we have to normalize the mode functions by a factor of L 3/2 . In our classical simulation we take a statistical point of view, interpreting the quantum creation and annihilation operator as stochastic variables that take different values at each realization. In this picture the creation and annihilation operators of eq. (4.12) are initiated as: where X l and Y l are random variables uniformly distributed between 0 and 1 for each l. This is equivalent to generating the Fourier modes of the field as Gaussian random numbers with variance |u l | 2 as follows:φ where u l is given by eq. (4.14) with a = 1 and τ = 0. From here, we first apply the iDFT eq. (4.3) and then add the background value of the inflaton to obtain the initial field configuration on the lattice. The fluctuations of the time derivative of the scalar field ∂ τφ l are generated in the same way using the time derivative of the mode functions ∂ τ u l and using the same realizations of X l and Y l . Note that we do not adopt the same procedure of LATTICEEASY for generating the initial field configuration, which is known to have a bug, as first noticed in [16]. However the Fourier transforms in our code are computed through the same routines of LATTICEEASY, contained in the script FFTEASY that is publicly available online [44].

B Outputs of the code
We now summarize the various outputs of the code. The background quantities are computed as averages over the box. We output quantities such as the average of the fieldφ, its derivativeφ , and the energy density and pressure of the field ρ lat and p lat . The code will also output the power spectrum of field fluctuations P φ ( k), which in all relevant applications will depend only on the absolute value of the momentum P φ (k). In order to do so, we first take the DFT to obtain |φ l1,l2,l3 | 2 . Then, after normalizing by a factor L −3 to get the physical power spectrum of the mode functions (see the discussion in appendix A.2), we average over spherical bins to obtain the one dimensional isotropic power spectrum P . This is done averaging |φ l1,l2,l3 | 2 over all lattice points such that int(| l|) = , where is the bin number. Then, a comoving momentum k lat, is associated to each bin by averaging the absolute value of eq. (4.2) over the bin. Note that the procedure for associating the momentum to each bin is different from the one of LATTICEEASY, where the momenta associated to the bins are simply 2π /L. This leads to a distortion in the output momenta of LATTICEEASY, which is independent of N and can lead to a difference of up to 20% in the IR 9 . We output the power spectrum for modes only up to the Nyquist frequency k Nyquist = 2π √ 3N Nyquist /L, where N Nyquist = N/2, because they contain all the physical information.
The dimensionless power spectrum P (lat) φ (k) is obtained multiplying P by k 2 lat, /(2π 2 ) (left plots of figs. 2 and 7). However, as we discuss in section 5.1, we can successfully reproduce the results of the continuous theory at all scales only if we multiply the dimension-full power spectrum P (k l ) by k 2 eff, /(2π 2 ) instead of k 2 lat, /(2π 2 ), where k eff is obtained averaging eq. (4.6) over the same spherical bins (right plots of figs. 2 and 7).

B.1 Energy conservation
In order to check energy conservation during the evolution of the system, we define the following quantity: We use this quantity, which should be close to 1, to quantify energy conservation in our code. In fig. 8 we show the plots of energy conservation in the two examples of section 5. From these plots we see that energy is conserved at 10 −7 level in both cases, and the same holds for all the examples discussed in this paper.

C Different stencils for the Laplacian operator
In this section we consider different stencils for the Laplacian. We refer to the Laplacian considered in the main text of eq. (3.3) and its corresponding effective momentum as L (2) [φ] i1,i2,i3 and k (2) eff , where the 2 refers to the second order of the stencil. The first one we consider is the following 4th order stencil, which has a similar structure of L (2) [φ] but involves more points: c a1,a2,a3 φ i1+a1,i2+a2,i3+a3 , (C.1) 9 Note that this distortion can also be relevant in generating the initial conditions. Next, we consider the anisotropic second order stencils defined in [45]. We display the coefficients associated to these stencil as: For each stencil L i we refer to its corresponding effective momentum as k i eff . We avoid writing the lengthy expressions for all the effective momenta, but we plot them in fig. 9 for a lattice with N = 128 and L = 1.4/m. All the k iso,i eff are real, with the exception of k iso,3 eff which becomes purely imaginary around k lat 435m (we show the absolute value of k iso,3 eff in the plot). From this plot we can see that only L (4) performs better than L (2) in terms of k eff,max and in terms of overall deviation from k lat , while the other isotropic stencils are significantly worse in this sense. The isotropic stencils, however, might perform better under other points of view. For example, these stencils do not have directional dependence in the second order truncation term in real space [45], contrarily to L (2) and L (4) .
In fig. 10 we show the final power spectrum computed from simulations with different stencils for the Laplacian. We run these simulations with the 1 2 m 2 φ 2 model and with the same parameters of section 5.1. We compare results from L (2) , L (4) and the isotropic stencil L iso,1 . This figure is analogous to fig. 2, and the dashed lines in the left plot are the analytical predictions for discrete k eff /m k lat /m dynamics computed from eq. (4.10). In all these cases, we can see that the identification k eff ↔ k allows us to recover the continuous result (right plot of fig. 10). The same result can be obtained with the other isotropic stencils L iso,i , that we do not show in order to make the plots more readable. For L iso,2 and L iso,3 , however, this is true only up to a certain momentum cutoff after which k eff (k lat ) starts decreasing, making the modes unphysical (see fig. 9).