The Effective Field Theory of Cosmological Large Scale Structures

Large scale structure surveys will likely become the next leading cosmological probe. In our universe, matter perturbations are large on short distances and small at long scales, i.e. strongly coupled in the UV and weakly coupled in the IR. To make precise analytical predictions on large scales, we develop an effective field theory formulated in terms of an IR effective fluid characterized by several parameters, such as speed of sound and viscosity. These parameters, determined by the UV physics described by the Boltzmann equation, are measured from N-body simulations. We find that the speed of sound of the effective fluid is c_s^2 10^(-6) and that the viscosity contributions are of the same order. The fluid describes all the relevant physics at long scales k and permits a manifestly convergent perturbative expansion in the size of the matter perturbations \delta(k) for all the observables. As an example, we calculate the correction to the power spectrum at order \delta(k)^4. The predictions of the effective field theory are found to be in much better agreement with observation than standard cosmological perturbation theory, already reaching percent precision at this order up to a relatively short scale k \sim 0.24 h/Mpc.


Introduction
Large Scale Structure Surveys have the potential of becoming the leading cosmological observable in the next decade. They contain a tremendous amount of cosmological information. If we were able to extract information from all the modes that go from the horizon scale ∼ 10 4 Mpc to the non-linear scale ∼ 10 Mpc, we would obtain about 10 4 10 3 ∼ 10 9 (1) independent modes. The Planck satellite in comparison has about (2×10 3 ) 2 ∼ 10 6 modes. Of course, accessing all this information is much harder than for the CMB, due to the short scale non-linearities. There are several aspects to this problem. The first problem is related to our currently limited understanding of the evolution of dark matter on large scales. Non-linear corrections are very important even on scales larger than 10 Mpc, because modes of different wavelengths couple to each other. Understanding these corrections is a problem that affects all large scale structure observables. There are then two additional issues that affect most, but at least not all, observables. One is the fact that most dark matter is clumped in very nonlinear structures (dark matter halos); and the other is the fact that what we often observe are galaxies, and not just dark matter halos, and not even dark matter long wavelength perturbations. The solution to these two last problems requires the correct understanding of the so-called halo-and galaxy-biases. These two problems, while important and deeply interesting in their own right, are very astrophysical in nature, and we do not address them here. Instead here we try to address in a rigorous way the first problem, that is the prediction of the dark matter distribution on scales larger than the non-linear scale. The fact that the universe is characterized by two well separated scales, the Hubble scale, over which perturbations are linear, and the non-linear scale, which indeed characterizes the scale over which gravitational collapse overtakes the expansion of the universe, makes the problem amenable to an Effective Field Theory (EFT) treatment. An effective theory is a description of a system that captures all the relevant degrees of freedom and describes all the relevant physics at a macroscopic scale of interest. The short distance (so called ultraviolet or 'UV') physics is integrated out and affects the effective field theory only through various couplings in a perturbative expansion in the ratio of microphysical UV scale/s to the macroscopic scale being probed. This technique has been systematically used in particle physics and condensed matter physics for many years, but has not been fully used in astrophysics and cosmology. An important early (and recent) application of these techniques in cosmology is the so-called Effective Field Theory of Inflation [1]. In a similar vein, understanding the large scale properties of the universe is very important, and is ready for a careful analysis.
Indeed the situation in the universe is very similar to what happens in the chiral Lagrangian that describes pion interactions in Particle Physics. At very low energies, pions are weakly interacting. These interactions and the size of the fluctuations grow with energy until we hit the Quantum ChromoDynamics (QCD) scale, ∼ 4πF π , at which the pions become strongly coupled. The Chiral Lagrangian [2] offers the correct effective theory allowing arbitrarily precise predictions, up to non-perturbative effects, at energies E 4πF π . In our universe, matter fluctuations are small at large distances and becomes larger and larger as we move up to the non-linear scale. Since the size of the non-linear terms, which are nothing but interactions, grows with the size of the fluctuations, we see that at long distances the universe should be described by some weakly coupled degree of freedom, that becomes more and more interacting as we move closer to the non-linear scale, at which point the fluctuations become strongly coupled. The coupling constant should indeed be represented by the ratio of the considered wavenumber k over the wavenumber at the non-linear scale k N L : k/k N L .
Notice that indeed the size of the density perturbations δρ/ρ on a scale k scales as (k/k N L ) 2 . This scaling suggests the existence of an effective field theory that should allow us to describe with arbitrary precision the universe on scales k k N L , very much as the Chiral Lagrangian represents the right effective field theory to describe pion dynamics to arbitrary precision.
Such an effective theory would have particularly relevant observational implications. Already now, large scale structure surveys such as BOSS or DES are measuring the galaxygalaxy correlation function, so called Baryon Acoustic Oscillations (BAO), at scales of order 100 Mpc. Next generation experiments such as LSST will measure this quantity at about percent precision. These observations contain huge amount of information on Dark Energy and on Inflation, through for example the non-Gaussianity of the primordial perturbations. The BAO scale is about one order of magnitude longer than the non-linear scale, where δρ/ρ ∼ 10 −2 , and therefore physics at this scale must be describable by rigorous perturbative methods. The alternative is to rely on either time consuming numerical simulations, or on analytical approaches that however are limited by some irreducible mistake that is hard to quantify precisely. In an ideal situation, numerical N -body simulations should be quickly done only at small scales, to describe phenomena affected by gravitational collapse, rather than running large simulations to describe weakly coupled physics. This has been indeed recently elucidated in the context of the bias, where it was shown that in order to derive the bias on large scales one needs to run very small simulations in a curved universe [3]. This line of reasoning is indeed very similar to what happens in QCD, where we perform lattice simulation to measure quantities relevant at energies above around one GeV, while we use the chiral Lagrangian for predictions at smaller energies.
The effective field theory (EFT) of the long distance universe was initially developed by some of us in [4]. It was noticed that by concentrating on length scales longer than the non-linear scale, the universe is described by a fluid with small perturbations. The equations of motion of this fluid are organized in a derivative expansion in the ratio of the considered wavenumber over the wavenumber associated to the non-linear scale k N L ∼ 1/10 Mpc −1 . At leading order in derivatives, the fluid has the stress tensor of an ordinary imperfect fluid, characterized by a speed of sound for the fluctuations, a bulk and a shear viscosity, plus a stochastic pressure component. This makes our approach different with respect to the 'standard' approaches both at a quantitative and a qualitative level.
The purpose of this paper is to further develop this effective theory and be able to make observational predictions. The parameters that characterize the fluid, the speed of sounds, bulk viscosity, etc., are determined by the microphysics at the non-linear scale, that we call UV, and cannot be derived from within the effective theory. They have to be either fit to observations, or measured in small N -body simulations. At this point, the EFT becomes predictive. Again, this is very similar to what happens in QCD, where one can measure the pion coupling constant F π in lattice simulations, after which the Chiral Lagrangian becomes predictive.
Our basic method and key results are summarized as the following: • By smoothing the collisionless Boltzmann equation for non-relativistic matter in an expanding FRW background on a length scale Λ −1 , we establish the continuity and Euler equations for an effective fluid. The Euler equation includes an effective stresstensor [τ ij ] Λ that is sourced by the short-modes δ s .
• By taking correlation functions of the stress tensor in the presence of long wavelength fluctuations, we define an effective stress-tensor that is only a function of the long wavelength fluctuations. It takes the form where the various parameters c 2 s , c 2 bv etc. are defined by proper correlation functions of short wavelength and long wavelength fluctuations.
• By directly evaluating the stress tensor from the microphysical theory, i.e., from Nbody simulations, and computing the appropriate correlation functions, we calculate the value of the fluid parameters. For a ΛCDM universe with standard cosmological parameters at redshift z = 0 and smoothing scale Λ = 1/3h Mpc −1 , we find where c 2 comb is the combination of c 2 s , c 2 bv and c 2 sv that is relevant for the leading nonlinear correction to the power spectrum (one-loop in perturbation theory), and c is the speed of light.
• Alternatively, by directly matching the couplings of the effective fluid to the measured power spectrum, we obtain c 2 comb (Λ = 1/3) 0.9×10 −6 c 2 in remarkable agreement with the direct measurement from N -body simulations.
• The fluid parameters carry Λ dependence (as does any 'bare' parameter in an interacting field theory). This cutoff dependence is taken to cancel against the cutoff dependence of the loop integral. As usual in effective field theories, we 'renormalize' the theory by sending the cutoff Λ → ∞ and carefully changing the fluid parameters so that predictions at low wavenumbers are not changed in the process. The finite values of the fluid parameters such as c 2 comb in the Λ → ∞ limit is a direct measure of the irreducible finite error made in standard approaches that approximate the dark matter on large scales as a pressureless ideal fluid. This is an irreducible error that is not recovered even by solving non-linearly the equations for a pressureless ideal fluid, as the various perturbative approaches attempt to do. This occurs simply because the equations they solve are not correct. Our approach, in contrast, should reach arbitrary precision, at least in principle, up to non-perturbative corrections.
• The pressure and viscosity dampen the power spectrum by acting in opposition to gravity, which makes sense intuitively. This is able to help explain the observed shape of the baryon-acoustic-oscillations in the power spectrum relative to standard perturbation theory (SPT).
• More precisely, at one-loop the density-density power spectrum receives a correction δP from the fluid parameters, which we find to be where P 11,l (k) is the linear power spectrum. Since this is negative and grows as a function of k, the power spectrum is reduced compared to SPT at high k's, improving the agreement with the full non-linear spectrum.
• We will find that already at one-loop, the computed power spectrum agrees at percent level with the non-linear one up to k ∼ 0.24h Mpc −1 . This suggest that in large scale structure surveys we should be able to extract primordial information all the way to at least such an high wavenumber, improving greatly with respect to the CMB our knowledge of the origin of the universe.

From Dark Matter Particles to Cosmic Fluid
We take dark matter to be fundamentally described by a set of identical collisionless classical non-relativistic particles interacting only gravitationally. This is a very good approximation for all dark matter candidates apart from very light axions. Given that on large scales baryons follow dark matter, we can include them in the overall dark matter description. As we discuss later, we also neglect general relativistic effects and radiation effects. In this approximation, numerical N -body simulations exactly solve our UV theory. The coefficients of our effective fluid can therefore be extracted directly from the N -body simulations, following directly the procedure described in [4]. Here, the UV theory is described by a Boltzmann equation. Therefore, in order to be able to extract the fluid parameters from N -body simulations, we need to derive the fluid equations from the Boltzmann equations and subsequently express the parameters of the effective fluid directly in terms of quantities measurable in an N -body simulation. This is the task of this section.

Boltzmann Equation
Let us start from a one-particle phase space density f n ( x, p) such that f n ( x, p)d 3 xd 3 p represents the probability for the particle n to occupy the infinitesimal phase space volume d 3 xd 3 p. For a point particle, we have The total phase space density f is defined such that f ( x, p)d 3 xd 3 p is the probability that there is a particle in the infinitesimal phase space volume d 3 xd 3 p: We define the mass density ρ, the momentum density π i and the kinetic tensor σ ij as The particle distribution f n evolves accordingly to the Boltzmann equation where φ n is the single-particle Newtonian potential. There are two important points to highlight about the former equation. First, we have taken the Newtonian limit of the full general relativistic Boltzmann equation. This is an approximation we make for simplicity. All our results can be trivially extended to include general relativistic effects. However, it is easy to realize that the Newtonian approximation is particularly well justified. Non-linear corrections to the evolution of the dark matter evolution are concentrated at short scales, with corrections that scale as k 2 /k 2 N L . General relativistic corrections are expected to scale as k 2 /(aH) 2 . This means that we should be able to cover up to wavelength of order 300 Mpc before worrying about per mille General relativity corrections. Furthermore, one of the main goals of this paper is to recover the parameters of the effective fluid of the universe from very short scale simulations valid on distances of order of the non-linear scale. The parameters we will extract in the Newtonian approximation are automatically valid also for the description of an effective fluid coupled to gravity in the full general relativistic setting.
A second important point to highlight in the former Boltzmann equation is about the single-particle Newtonian potential φ n . Following [4], the Newtonian potential φ is defined through the Poisson equation with ρ b being the background density and ∂ 2 = δ ij ∂ i ∂ j . We raise and lower spatial indexes with δ ij . The solution reads Notice that the overall φ( x) is IR divergent in an infinite universe. This is due to a breaking of the Newtonian approximation. We have regulated it with an IR cutoff µ that we will take to zero at the end of the calculation. Our results do not depend on µ, as indeed we are interested in very short distance physics. By summing over n, we obtain the Boltzmann equation for f

Smoothing
Following [4], we construct the equations of motion for the effective fluid by smoothing the Boltzmann equations and by taking moments of the resulting long-distance Boltzmann equation. The smoothing guarantees that the Boltzmann hierarchy can be truncated, leaving us with an effective fluid. indeed, notice that it is not trivial at all that we should end up with an effective fluid. Fluid equations are usually valid over distances longer than the mean free path of the particles. But here for dark matter particles the mean free path is virtually infinite. What saves us is that the dark matter particles have had a finite amount of proper time, of order H −1 , to travel since reheating, and they traveled at a very non-relativistic speed. This defines a length scale vH −1 ∼ 1/k N L which is indeed of order of the non-linear scale. This length scale plays the role of a mean free path, as verified in [4]. The truncation of the Boltzmann hierarchy is regulated by powers k/k N L 1. We define the Gaussian smoothing with Λ 2 representing a k-space, comoving cutoff scale. This will smooth out quantities with wavenumber k Λ, or equivalently with waveleghts smaller than λ 1/Λ. We regularize our observable quantities O( x, t), ρ, π, φ, . . . , by taking convolutions in real space with the filter, defining long-wavelength quantities as Notice that in Fourier space W (k) → 1 as k → 0: our fields are asymptotically untouched at long distances. The smoothed Boltzmann equation becomes Fluid equations are obtained by taking successive moments creating in this way a set of coupled differential equations known as Boltzmann hierarchy. As we will explain in more detail later, it will be sufficient for the purposes this paper to stop at the first two moments (one-loop approximation). The first two moments will give the continuity and momentum equations in the approximation in which the fluid is described by the Navier-Stokes approximation, with the addition of a stochastic term. We obtaiṅ Let us define the various quantities that enter in these equations. We define the long wavelength velocity field as the ratio of the momentum and the density The right hand side of the momentum equation (19) contains the divergence of an effective stress tensor which is induced by the short wavelength fluctuations. This is given by where κ and Φ correspond to 'kinetically-induced' and 'gravitationally-induced' parts: Note that we have subtracted out the self term from w ij l , as necessary when passing from the continuous to the discrete description in the Newtonian approximation, and used that ∂ 2 φ = 4πGa 2 (ρ − ρ b ) and ∂ 2 φ l = 4πGa 2 (ρ l − ρ b ) to express Φ l in terms of φ and φ l . In the limit in which there are no short wavelength fluctuations, and Λ → ∞, κ l and Φ l vanish. In App. A we provide the above expression written just in terms of the short wavelength fluctuations.

Integrating out UV Physics
The effective stress tensor that we have identified is explicitly dependent on the short wavelength fluctuations. These are very large, strongly coupled, and therefore impossible to treat within the effective theory. When we compute correlation functions of long wavelength fluctuations, we are taking expectation values. Since short wavelength fluctuations are not observed directly, we can take the expectation value over their values. This is the classical field theory analog of the operation of 'integrating out' the UV degrees of freedom in quantum field theory, now applied to classical field theory. The long wavelength perturbations will affect the result of the expectation value of the short modes, through, e.g., tidal like effects. This means that the expectation value will depend on the long modes. In practice, we take the expectation value on a long wavelength background. The resulting function depends only on long wavelength fluctuations as degrees of freedom. In this way, we have defined an effective theory that contains only long wavelength fluctuations. Since long wavelength fluctuations are perturbatively small, we can Taylor expand in the size of the long wavelength fluctuations. Schematically we have For the precision we pursue in the rest of the paper, we will stop at linear level in the long wavelength fluctuations, though nothing stops us from going to higher order. By the symmetries of the problem, the resulting stress tensor must take the following form This is the stress tensor of an imperfect fluid. p b is the background pressure that is induced by short distance inhomogeneities even in the absence of long wavelength fluctuations. c 2 s is the speed of sounds of the fluctuations: δp = c 2 s δρ. The parameters c bv and c sv are the coefficients for the bulk ζ and the shear η viscosity respectively, with units of velocity. They are related to η and ζ by the relation η = 3ρ b c 2 sv /(4H), ζ = ρ b c 2 bv /H . ∆τ ij represents a stochastic term, that takes into account the difference between the actual value of τ ij in a given realization and its expectation value 1 . We will come back to this term shortly, but it is worth noting that neglecting this term in the above equations reproduces the familiar Navier-Stokes equations.
Finally, the ellipses (. . .) represent terms that are either higher order in δ l , or higher order on derivatives of δ l . Indeed, higher derivative terms will be in general suppressed by k/k N L 1, and, as typical in effective field theories, we take a derivative expansion in those. Astrophysically, these terms would corresponds to the effects induced by a sort of higherderivative tidal tensor. Once we expand in derivatives of the long wavelength fluctuations, we take the parameters in (25) to be spatially independent, but time dependent.
The coefficient δp b , c s , c sb , c sv are determined by the UV physics and by our smoothing cutoff Λ, and are not predictable within the effective theory. They must be measured from either N -body simulations, or fit directly to observations. This is akin to what happens in the Chiral Lagrangian for parameters that can be measured in experiments or in lattice simulations, such as F π . We first define the correlation functions that will allow us to extract these parameters from small N -body simulations.

Matching Correlation Functions
It useful to define the following quantities from the stress tensor Then, according to (25), we have In order to extract the parameters of the effective fluid, we multiply each of these functions with long wavelength fields, and take expectation values. By forming suitable combinations of these, the parameters of the effective fluid can be extracted. We will need the following set of correlation functions where δ l = δρ l /ρ b . From these we obtain the following expressions for the parameters of the effective theory where c 2 v = c 2 sv + c 2 bv is the sum of the viscosity coefficients. By extracting the correlation functions in (29) from N -body simulations, and performing the ratios in (30), we should be able to extract the parameters of the effective theory. Notice that the ratios are supposed to be spatially independent. Such behavior is expected to hold at large distances x Λ −1 where the higher derivative terms are negligible.
In simulations, we should in principle also measure the stochastic components of the stress tensor. In the two point function at one-loop, at leading order in derivatives, it enters just the correlation function of the trace. This amounts to measuring We will see that the effect of this stochastic term is accidentally higher order in δ l and so does not enter at leading order. After all these parameters have been measured in N -body simulations, the EFT is prone for perturbation theory. It is alternatively possible to perform directly perturbation theory and fit the results to observations. We will be able to perform both approach and check that we obtain the same result. We will describe in detail how to measure these quantities in simulations in App. D and we will give the results of these measurements in Sec. 4. For the moment we will instead directly move to apply perturbation theory with our EFT.
These parameters can either be measured from N -body simulations directly or kept generic and then extracted by fitting the results to observables. We are able to do both and verify that we obtain the same result. We will describe in detail how to measure these quantities in simulations in App. D and we will give the results of these measurements in Sec. 4. First we develop perturbation theory within the EFT that will allow us to make predictions and to extract these parameters from observations.

Perturbation Theory with the EFT
We now proceed to perform perturbation theory within our EFT. The non-linear equations of motion that we need to solve are , Ω m is the present day matter fraction, a 0 is the present day scale factor, usually taken to be equal to 1, and . . . represent higher order terms (in sense that we will explain shortly) that come from the expression of the short wavelength stress tensor τ ij in terms of long wavelength fluctuations. Our theory is defined on scales longer than the non-linear scale. For this reason we have δ l 1. The longwavelength velocity v l and the long-wavelength φ l are small even inside the the non-linear scale and they are even smaller at larger distances. This means that we can reliably solve the above non-linear equations iteratively around the linear solution. Such an iterative solution is very similar to what is done in quantum field theory, where the solution to the quantum non-linear equations is organized in Feynman diagrams. Indeed we can organize the various perturbative terms around Feynman diagrams even in this case. The result is very similar to what is computed in the in-in formalism, for example when computing quantum corrections to inflationary correlation functions [26]. Indeed, the calculation we are going to do shares many of the features that are present in normal quantum field theory computations: cutoff, renormalization, running, and so on are all concepts that will appear and prove useful as we proceed. They have nothing to do with the word 'quantum' in 'quantum field theory', rather they have to do with the 'field theory'. Our calculation is for a classical field theory and shares all these features.

Organization of the perturbation theory
The simplest way to organize our perturbation theory is to use the fact that, in any order of magnitude approximation, φ is constant at all scales, of order 10 −5 , and that well inside the horizon the Newtonian approximation holds. For length scales longer than the equality scale, at the linear level we therefore have We see that as L → 0, δ l grows and indeed becomes of order one at L ∼ λ N L 2 . At distances larger than the non-linear scale, we therefore expand in powers of δ l , keeping in mind that the additional fluctuations scale as in (33). Let us estimate the relative size of the terms.
Loop corrections: it is easy to estimate from the non linear structure of the equations that ∂ i v i ∼ Hδ l . Notice that in the power spectrum we need to take two non-linear corrections, or alternatively look at the cubic corrections. The non-linear terms scale as Loop corrections therefore scale as δ 2 l , peaked at the highest possible scale within the theory Λ. Pressure and Viscosity terms: these terms result from integrating out the modes higher than the Λ scale. So, naively they should scale as the δ at the extreme UV scales beyond the effective theory. Since the theory in the UV is strongly coupled, very large corrections are expected, and the result cannot be extrapolated from the linear regime, even at the order of magnitude level. This is why we will measure the parameters such as c 2 s from N -body simulations. What we will find is that these parameters are of order 10 −5 . This happens because the combination of short modes that generates the parameters like c s are such that short modes that have virialized do not contribute 3 . Since these terms scale as φδ, the contribution is peaked at those modes that have just become non-linear δ ∼ 1, but not yet virialized. We therefore expect the parameters c 2 s 's to be of order φ ∼ 10 −5 . We will see that the fact that short modes entered the horizon in the radiation era makes this a bit of an overestimate. At this point we are ready to estimate the size of these corrections: Notice that thanks to the strongly coupled UV theory (or thanks to the fact that non-linear structures virialize), we have that for c 2 s ∼ 10 −5 , the contribution from these terms is larger than the one loop contribution in the low energy theory. This is so because the theory is strongly coupled in the UV. We conclude that one insertion of these terms counts at least as a one-loop term.
Stochastic terms: Let us now continue on to evaluate the effect of the stochastic terms. This is a bit more complicated. Let us evaluate the relative effect on the power spectrum. The structure of the equations leads to the following approximate non-linear solution In the power spectrum we therefore have as the stochastic part must be correlated with itself. Due to virialization, we expect that the correlation function of τ should be Poisson like on independent pixels of order the non-linear scale k −1 N L . We therefore estimate (38) This is indeed confirmed by calculations in perturbation theory [4]. We therefore have This tells us that on scales longer than the non-linear scale, the contribution of the stochastic pressure is parametrically smaller than the pressure effects. Since in this paper we will stop at one-loop order, we can therefore neglect this correction. It should be noted that these contributions scale parametrically differently than the loop contributions, and so, depending on the scale considered, they might be more relevant that a 2-loop contribution. We notice that something similar happens also at the level of dissipative fluids, where we generically include dissipative terms through the Navier-Stokes equations, but we neglect stochastic terms. Higher Derivative Terms: When we take the expectation value of the short-distance stress tensor in a background of a long mode, we have Taylor expanded both in the long wavelength fluctuations and in their derivatives. Higher power corrections scale with powers of δ l , and so corresponds to higher-loop terms. Higher derivative terms instead scale as powers of k 2 /k 2 N L ∼ δ l , where we have taken the squared because of rotational invariance. This shows that higher derivative terms scale nicely as loop terms. If we allow the cutoff to remain finite, we should also include higher derivative terms that scale as k 2 /Λ 2 .
General Relativistic and Radiation Corrections: In this paper we will neglect general relativistic corrections and all non-linear contributions coming from the fact that the universe was radiation dominated at early times. General Relativistic corrections scale as For the high scales where non-linear corrections are relevant, for example at the BAO scale k BAO ∼ 10 −2 , these corrections are of order 10 −4 , and so uninteresting from this point of view. Radiation is the dominant component of the universe at early times. Neglecting it amounts to neglect corrections that scale as a/a eq ∼ 10 −3 , where a is the scale factor and a eq is the scale factor at matter radiation equality. Inclusion of these corrections in perturbation theory has been studied in [20], and it gives a small correction to power spectrum, and small, but potentially measurable, corrections to the three-point function, corresponding to f N L ∼few. Both General Relativistic effects and radiation effects do not represent an intrinsic limitation of our EFT. They can be straightforwardly included in our formalism, by simply improving the equations of motion we use in this paper.
In summary, we see that apart from the stochastic terms and some higher derivative terms, all the remaining terms: loops, pressure, higher derivatives and higher powers of δ l from τ , scale as powers of δ l , which is our main ordering parameter. Stochastic terms instead contribute at leading order as δ 3/2 l , so they count as one loop and a half. Cutoff-dependent higher derivative terms scale as k 2 /Λ 2 .
Cutoff dependence and effective expansion parameter: So far, it looks like that our expansion parameter is the highest δ l we have in our theory, which is δ l (k ∼ Λ) ∼ Λ 2 /k 2 N L . However, the situation is even better than this. So far, we have defined our theory with a regulating cutoff at k ∼ Λ. Because of this, all our intermediate results depend explicitly on Λ: c s (Λ), c sv (Λ), etc. and loops need to be cutoff at Λ. This induces an explicit Λ dependence plus higher derivative terms of order k 2 /Λ 2 . However, the sum of all the diagrams will be independent of Λ. Indeed c s (Λ), etc. should be really thought of as oneloop counterterms. Upon carefully choosing the counterterms c s (Λ), etc., any Λ dependence cancels apart from terms in k/Λ that should be removed by higher derivative corrections in the stress tensor that we neglected. In order to resolve such error with the least effort, we will choose the counterterms at a fixed cutoff in such a way as to have the theory agree with observations at a certain renormalization scale k ren. , and then we will extrapolate our results to Λ → ∞, effectively letting the residual terms in k/Λ vanish. In this Λ → ∞ regime, the loop term gets dominated by the regime in which one of the modes has wavenumber of order the non-linear scale, while the other has a wavenumber of order of the external wavenumber. In this way, loop terms scale as one-power of δ l as the counterterm 4 . At this point, the expansion parameter of the EFT will be δ l ∼ k 2 /k 2 N L evaluated at the scale of the external modes, with no residual Λ dependence even in the expansion parameters.
Again, this is very similar to what happens when one computes loop corrections in the Chiral Lagrangian. After regulating the chiral theory with a cutoff Λ, there are naively two expansion parameters. If E is the energy scale of the process, we have E/F π and E/Λ. After renormalization and by sending Λ → ∞, we are left only we E/F π as an expansion parameters.
In summary, the expansion parameter of the EFT is δ

One-loop Perturbation Theory
We are now ready to implement perturbation theory for the power spectrum at quartic order in δ l , that is at one-loop. At this order, the equations we are going to solve are the ones in (32) with ∆J and the . . . terms neglected.
Let us write the equation for the vorticity w i l = ijk ∂ j v k . Neglecting the stochastic terms that we argued are small, we have In linear perturbation theory the vorticity is driven to zero, and this occurs even the more so at this order in perturbation theory, as the source is proportional to w l . While at higher order one could expect vorticity to be generated, at this order, and therefore for the purposes of this paper, we can take it to be zero. This means that we can work directly with the divergence of the velocity Using a as our time variable, the equations (32) reduce to aHθ l + Hθ l + 3 2 where H = a −1 ∂a/∂τ , subscript 0 for a quantity means that the quantity is evaluated at present time, we have set a 0 = 1, represents ∂/∂a and As we discussed, the parameters c s , c bv and c sv are time dependent and must be measured in the simulations as a function of time. For the purposes of this paper, we will make the simplifying assumption that their time dependence can be inferred in perturbation theory. In other words, we will measure them at one time and deduce their values at different times by perturbation theory 5 .

Perturbative Solutions
Since the correlation function of matter overdensities is small at large distances, we can solve the above set of equations (43) perturbatively in the amplitude of the fluctuations. For the computation of the power spectrum at one loop, it is enough to solve these equations iteratively up to cubic order. Order by order, the solution is given by convolving the retarded Green's function associated to the linear differential operator with the non-linear source term evaluated on lower order solutions. At second order we obtain Let us explain some of the relevant expressions that appear here. G(a,ã) is the retarded Green's function for the second order linear differential operator associated with δ that is obtained after substituting θ in the second equation of (43) with the value obtained from the first, and linearizing. In doing this, it is important to neglect all the terms of order c 2 s because, in our power counting, they count as non-linear terms. This is given by G(a,ã) = 0 for a <ã .
5 As we will see, these parameters need to cancel the Λ dependence associated to the regularized loops. The part of these parameters that depends on Λ can be therefore reliably inferred in perturbation theory. However, the part of these parameters that is Λ-independent and that represents the finite contributions should be measured in simulation or in observation. We will assume that the time dependence for these two components is the same. We will check that this is an accurate approximation in an upcoming paper [27]. For this approximation, we stress that since these are 1-loop terms, it is important to know them up to a relative factor of order δ l 1.
For a ΛCDM cosmology the result can be expressed 6 as a hypergeometric function, although its form is not particularly illuminating. For all calculations presented here it is sufficient to numerically solve the above differential equation. This can be easily accomplished by replacing the δ(a −ã) on the RHS of the first equation with zero, but starting with the boundary conditions being G(a,ã)| a=ã = 0, and ∂ ∂a G(a,ã)| a=ã = 1/(ãH(ã)) 2 . In principle, it is possible to include in the linear equations that determine the Green's function and the growth functions also the higher-order linear terms proportional to c 2 s and c 2 v . Doing this amounts to resumming the effect of these pressure and viscous terms. The resulting linear equation can be easely solved numerically, finding for example that the growth factor becomes k-dependent, being the more suppressed the higher is the wavenumber [21]. However, it is not fully consistent to resum these terms without including the relevant loop corrections.
D(a) represents the growth factor at scale-factor-time a. In particular, we have written the linear solution as with a 0 being the present time, and δs 1 representing a classical stochastic variable with variance equal to the present smoothed power spectrum with P 11,l (k) being the smoothing of the linearly computed power spectrum at present time A very useful simplification is due to the fact the growth factor and the Green's function are k-independent. This is due to the fact that at linear level we can neglect the pressure and viscosity terms that would otherwise induce a k-dependence. Because of this, the convolution integrals that would couple time integration and momentum integration nicely split into separate time integrals and momentum integrals that can be simply performed separately. We have tried to underline this in (45) by adding suitable parenthesis. Iterating, we obtain the solution for δ at cubic order δ (3) . For brevity, we report it in App. B. Notice that in the terms in δ (2) (and δ (3) ) we have neglected the contribution from the pressure and the viscosity, which count as third order terms. They give where c 2 comb is given by and it is the combination that is relevant at one-loop order. For the terms multiplying c 2 s and c 2 v in the second equation of (43), we can substitute the linear relation Notice that δ (2) and δ (3) are the same as the standard ones used in perturbation theory, with just two differences. The first is the important smoothing of the sourcing power spectrum. This makes the convolution integral that we are going to perform next rapidly converging, but also Λ dependent. The second difference is of a more technical nature, and relies on the fact that in standard perturbation theory the time dependence of the non-linear solution is approximated by the linear growth factor D elevated to the power 2 and 3 for δ (2) and δ (3) , while the momentum dependence is approximated to be the same momentum dependence as in standard EdS universe. This procedure is exact in EdS, but not so in other space times. Some studies [18] (see also [29,30,31]) have checked that this is correct up to percent level on the full power spectrum. Since however percent accuracy is the target of next generation experiments, we decide to perform the correct computation, which is not so very complicated to set up in any case. For the purpose of comparing with the literature and to gain familiarity with the EFT setup with simpler formulas, we provide results obtained with this approximate treatment of the perturbed solutions in App. C.

Diagrams
By contracting the non-linear expression we obtain the non-linear corrections. There are three diagrams at order δ 4 l . After including the linear contribution, we have l ( k, a 0 )δ (2) l ( q, a 0 ) , P 13 (k, a 0 ) = 2 δ where the . . . means that we have removed a factor of (2π) 3 δ (3) ( k + q) from the expectation value. P 11 represents the unsmoothed linear power spectrum, as the linear theory does not need to be regularized. The term P 13, c 2 comb is supposed to remove the Λ dependence that comes from P 13 . It is a counterterm diagram. Strictly speaking, we would need a counterterm diagram also from P 22 , which is provided by the two-point function of the stochastic source ∆J i in (32). As we discussed, this term is supposed to count as a δ 5 l term, and therefore we neglect it. This means that the Λ dependence associated with P 22 is very weak at this order in the calculation. The full stochastic term will be included in a following paper [27].
The expressions for P 22 , P 13 , P 13, c 2 comb are given by where cos(θ) = k · q/(k q); The convolution integrals in P 22 and P 13 are the sign that these are one loop diagrams. P 13, c 2 comb does not have a convolution integral as it is a one-loop counterterm. The diagrams are pictorially represented in Fig. 1. Since α( k, − k) = β( k, − k) = 0, there are no non-1PI diagrams. Notice that in P 13 and P 22 we have already carried out at least a part of the angular integration. The sign of the P 13, c 2 comb in (57) is also quite intuitive. For positive c s or c v , the gravitational collapse is slowed down, and so this contribution tends to decrease the gravitational collapse.

Cutoff-(in)dependence
Each diagram is dependent on the cutoff Λ: P 22 and P 13 through the smoothed linear power spectrum, while P 13, c 2 comb depends on Λ through c 2 comb , which is Λ dependent because it arises from integrating out the short distance fluctuations. The Λ dependence of P 13 is to be cancelled by P 13, c 2 comb , while the one of P 22 from the stochastic fluctuations. Consider the sum of the P 13 and P 13, c 2 comb terms. In order for this sum to be Λ independent, both must have the same k-dependence in the relevant regime. By inspection of (57), we see that P 13, c 2 comb goes as k 2 , which implies that P 13 should behave in the same way in this regime. This is in fact the case, as can be readily verified by taking the k → 0 limit of (56). In particular, we can define a Λ-independent renormalized parameter c comb;ren. defined at a renormalization scale k ren. and a Λ-dependent counterterm parameter c 2 comb;ctr. (Λ) such that (a, Λ). However, for the purposes of this paper, we can approximate them to be equal, and we will check this approximation in a forthcoming paper [27]. We can therefore extract the time dependence of c 2 comb (a, Λ) from the k → 0 limit of P 13 . We obtain c 2 comb;ctr. (a, Λ) = c 2 comb;ctr. (a 0 , Λ) where and A plot of the time dependence of the speed of sound is given in Fig. 2.  comb as inferred using the correct time dependence from P 13 and instead using the approximate time dependence derived from the growth functions (see App. C). Starting from very early times, we see that c 2 comb grows as a functions of time, peaks at about a 0.7, and then decreases near the present epoch, probably as due to the onset of dark energy. c 2 comb is positive, implying that this term tends to slow down the collapse of structures.
We determine the value of c 2 comb;ctr. (a 0 , Λ) in two different independent ways: one involving fitting to observation of the power-spectrum derived from simulations, and the other involving direct measurement from simulation. In the first, simplest, way, we determine c 2 comb;ctr. (a 0 , Λ = ∞) by matching the one-loop EFT power spectrum at k = k ren. to the power-spectrum extracted from simulations (or directly from precise observations in the future!). We can do this at various values of Λ, but we take the Λ → ∞ limit in order to drive to zero any effect from higher derivative terms down by powers of k/Λ. At this point, we can derive c 2 comb (a 0 , Λ = ∞) by running the value at Λ = ∞ down to a finite Λ. The formula is given by The limit k ext → 0 is necessary in order to suppress higher derivative terms down by powers of k ext /Λ. The running of c 2 comb is plotted in Fig. 3. c 2 comb (a 0 , Λ = ∞) 6.2 × 10 −7 and is the value obtained by fitting to data using k ren. = 0.16 h Mpc −1 7 . We see that as Λ decreases, we integrate out more and more modes, and c 2 comb grows. We see that the result matches with the one obtained by the second method for measuring c 2 comb , that is by using N -body simulation to extract directly the correlations in (30). We perform measurements at Λ = 1/3 h Mpc −1 and Λ = 1/6 h Mpc −1 , and we see that the match is extremely good. We describe more precisely how these measurements are derived in Sec. 4. We take this as an extremely promising indication of the strength of our approach.
In order to elucidate the effect of the higher derivative terms, we plot in Fig. 4 the value of c 2 comb (a 0 , Λ = 1/3), for various values of the external k ext . It is only for vary low k ext 's that c 2 comb becomes k ext. independent. Finally, there is a third method in which we could have derived c 2 comb (a 0 , Λ). By keeping Λ finite, we could have fit the analytical results to N -body simulations by including higher derivative terms proportional to powers of k/Λ. Indeed, unless Λ → ∞, the largest of these terms are not negligible and need to be included to get the correct c 2 comb (a 0 , Λ). A description of this approach in detail is given in App. E, and leads to the same results for c 2 comb (a 0 , Λ).

Fluid parameters from N -body simulations
If this language of effective field theory is to be born out, we must be able to take the fundamental theory, integrate away UV effects, and find agreement in terms of the parameters described above. Fortunately we have, in the form of simulation, exactly those calculations in the fundamental theory. By smearing the positions of simulated particles with a normalized Gaussian function of width Λ, we are able to introduce a soft UV cutoff of order 1/Λ. For correlations on scales longer then the cutoff we can directly measure c 2 comb corroborating the perturbative analysis presented above.
Specifically we consider random downsamples of positions and velocities of the Consuleo simulation 8 [28] from 2.7 × 10 9 particles downsampled to 1, 000, 053 particles distributed over (420 Mpc) 3 . Even this incredibly coarse resolution allows us to measure the following fields The dependence on the renormalization scale is a measure of the importance of higher loops. We see that as Λ → ∞, c 2 comb decreases as more and more modes are included within the regime of validity of the EFT. However, the fact that as Λ = ∞, c 2 comb = 0 is an indication of the fact that the fundamental theory is not described by a pressureless ideal fluid, but by indeed freely streaming dark matter particles. Data points with 1 σ error bars represent the value obtained from N -body numerical simulations using the methods described in Sec. 4 using two different smoothing lengths Λ −1 . Given the error bars from numerical evaluation, the measured values are in remarkable agreement with what inferred from renormalizing using the power spectrum. δ l , Θ l , ∂ 2 δ l , ∂ 2 Θ l , and A s to measure c 2 comb to within standard errors of 10 percent. While the complexity of the first four fields go linearly in the number of particles, A s is more expensive. The fact that we can achieve such consistency with such a small number of particles is not only remarkable, but numerically quite convenient. We describe the details of the analysis in App. D.1, and here simply provide a summary of the results.
The spatial dependence of measured c 2 comb are plotted for Λ = 1/3 h Mpc −1 and Λ = 1/6 h Mpc −1 in Fig. 5. These regions were chosen to maximize numerical stability as described in App. D. Fitting a constant to the value of c 2 comb for the displayed values after the UV cutoff gives in units of c 2 , we obtain c 2 comb (Λ = 1/6) = 1.26 ± .1 × 10 −6 .
The RG flow between the two measured values is consistent with the prediction from perturbation theory as seen in Fig. 3. Furthermore the measured value from Consuelo agrees nicely with that predicted from matching to the nonlinear CAMB power spectrum.  (62). We see that only as k ext → 0, c 2 comb becomes k ext independent. This is so because at high k any higher derivative terms suppressed by powers of k/k N L are important.

Results
The O(δ 4 l ) result of the computation of the power spectrum with our EFT is presented in Fig. 6. On top, we plot the ratio of the one loop power spectrum compared with the non-linear fit provided by the CAMB software with high precision settings, evaluated with the following cosmological parameters: Ω Λ = 0.75 , Ω m = 0.25, Ω b = 0.04, h = 0.7, n s = 1. The linear power spectrum is also obtained from CAMB with high precision settings. We take these data for all perturbative calculations done in the paper. Often in the literature the results of perturbation theory are plotted as ratio of the perturbation theory result versus a no-wiggle power spectrum. In this way the oscillatory features are still present in the plot, though they come mostly from the linear theory. We give this in the bottom part of the plot. l prediction from our EFT is compared with the CAMB non-linear output in the top, and to the no-wiggle power spectrum in the bottom, as well with the linear theory and Standard Perturbation Theory (SPT). The results from the EFT agree at percent level with the nonlinear theory up to k 0.24h Mpc −1 , when some high scale power seems to be missing. Results should improve already by going to δ 5 l order. The results are remarkably better than using SPT. The no-wiggle power spectrum we use is given by P δδ,No−Wiggle = 5.1 · 10 6 q log 2 (13q + 2e)/(54 q 2 (14 + 731/(457q + 1)) + log(13q + 2e)) 2 .
This plot is obtained after renormalizing the EFT prediction to match the power spectrum at k ren. = 0.16h Mpc −1 . In detail, we perform the calculation at several increasing values of Λ, we choose the value of c 2 comb (a 0 , Λ) to match the simulations' non-linear output, and then we extrapolate to Λ → ∞. This gives us c 2 comb (a 0 , Λ = ∞) 6.2 × 10 −7 . Notice that naive considerations of virialization gave an estimated value for c 2 comb of order 10 −5 [4]. The obtained smaller numerical value fits well with the decrease in the transfer functions for wavenumbers that are higher than the equality scale.
The result for the power spectrum agrees at percent level with the CAMB non-linear fit up to k 0.24h Mpc −1 , where the EFT prediction begins to be smaller than the N -body simulation result. Results obtained with the approximate time dependence described in App. C are close to these ones, at percent level. As discussed in App. C, it is potentially dangerous to trust this approximation at high k's for percent level precision, and luckily it is not very hard at all to perform the correct perturbation theory. We expect that the inclusion of the stochastic pressure and of higher order diagrams should improve the fit in the UV, possibly allowing us to fit the simulations to even higher k's. Notice how the counterterm P 13, c 2 comb decreases the power spectrum, compensating for the overshooting of SPT. It is difficult to interpret the percent disagreement that we have at moderate slow scales such as at k 0.12h Mpc −1 . At face value, it looks like that the computed power spectrum presents oscillations that are too large. These disagreements might be improved with the inclusion of higher order terms, or it might even be that at this level of precision, the results from the CAMB non-linear fit or from the N -body simulations might not be precise enough. These same improvements should reduce also the dependence on the renormalization scale: by changing the renormalization scale our results change by about 2%, at k ∼ 0.24h Mpc −1 . This dependence can be taken as a measure of the contribution from higher order terms.

Conclusions
Large scale structure surveys have the potential of becoming the next leading observational window on the physics of the early universe, potentially greatly improving what we are already learning from the CMB. Large Scale Structure physics is however much more complicated than the CMB due to the presence of large matter clustering at small scales. Since at non-linear level different scales are coupled, these non-linearities affect even large scale perturbations that are mildly non linear and so potentially treatable in a perturbative matter. In this paper we have developed the effective field theory of cosmological large scale structures in order to achieve a reliable predictability. Calculations in the effective theory are performed in k/k N L . The effective field theory is a cosmological fluid description for cold dark matter, and by extension all matter including baryons which trace the dark matter. The microphysical description is in terms of a classical gas of point particles, which we have smoothed at the level of the Boltzmann equation. We have exhibited and computed the various couplings that appear in the effective field theory, namely pressure and viscosity by matching to N -body simulations, finding c 2 s ∼ 10 −6 c 2 , etc. . We have developed the perturbative expansion for the power spectrum, which we have carried out at the O(δ 4 l ). The fluid parameters arise from UV modes and alter standard perturbation theory. We have found that the corrections lead to a power spectrum in percent agreement with the full nonlinear spectrum as obtained by CAMB up to k 0.24h Mpc −1 .
It is a peculiar coincidence that k 0.24h Mpc −1 is also the maximum k at which the popular technique of Renormalized Perturbation Theory (RPT) [9] works. While RPT is a very nice technique to compute non-linear corrections to the power spectrum, we stress that our approach is different at a qualitative and a quantitative level. At a qualitative level, RPT tries to solve, as exactly as possible, non-linear equations for a pressureless ideal fluid. In our approach, instead, we try to solve non-linear equations for a different fluid. This has quantitative effects, as it is shown from the fact that as Λ → ∞ our effective parameters like c 2 comb do not vanish. What is more important about our EFT is that in can be improved. By performing higher order computations and by adding suitable counterterms, in principle arbitrary precision for reconstructing the power spectrum, or indeed any dark matter observables, can be achieved by going to a sufficiently high order in perturbation theory, on scales k k N L . Techniques such as RPT or the renormalization group approach [15] for example are still very nice techniques to perturbatively solve some non-linear equations, resumming many diagrams. It would be interesting to apply those techniques to solve the equations of our EFT. We leave this to future work.
The effective field theory approach to large scale structure formation is complimentary to N -body simulations by providing an elegant fluid description. This provides intuition for various nonlinear effects, as well as providing computational efficiency, since the numerics required to measure the fluid parameters are expected to be computationally less expensive than a full scale simulation. Indeed we are able to achieve excellent agreement with the small down-sample of the full simulation we examined here. Of course, since the couplings are UV sensitive, it still requires the use of some form of N -body simulation to fix the physical parameters, either by matching to the stress-tensor directly or to observables. But this matching is only for a small number of physical parameters at some scale and then the constructed field theory is predictive at other scales.
There are several possible extensions of this work. A first extension is to go beyond the one-loop order to two-loop, or higher. This will require the measurement of several new parameters that will enter the effective stress-tensor at higher order, including its stochastic terms. Another extension is to compute the velocity fields and to include the small but finite contributions from vorticity, or to compute higher order N -point functions, which can probe non-Gaussianity. Finally, another extension is to consider different cosmologies; in this work we have presented results on dark energy in the form of a cosmological constant. But one could equally consider other models for dark energy. This would presumably alter the value of the fluid parameters in a way that could be measured either from new simulations or determined by observation. Hopefully, our effective field theory for large scale structures will help us use large scale structure surveys to uncover the physics of the beginning of the universe.
Similarly, one can prove that Φ ij l satisfies So altogether we obtain the effective stress-tensor where We see that [τ ij ] Λ is sourced by short wavelength fluctuations plus higher derivative corrections. Note that by taking the derivative ∂ j this leading piece becomes where the first term in (76) is given by and the second term in (76) can be expanded as and this term should be included since it involves the background piece ρ b , and so it includes a first order contribution.

B Expression for δ (3)
The iterative solution for δ (3) is given by where D i 's represent the result of the integration of the Green's functions and the other time-dependent coefficients. They are given by Notice again the great simplification that occurs due to the fact that the growth factors and the Green's function do not depend on k, so that the time integrals and the momentum integrals decouple.
C O(δ 4 l ) Power Spectrum with Approximate Treatment In the main part of the paper, we performed perturbation theory with our EFT in a rigorous and exact way. However, it is possible to perform an approximate treatment that simplifies quite a bit the formulas. In this way it is simpler to follow the derivation and we therefore present it here. If the universe were to be EdS, then the solution for δ (n) would be δ (n) ∝ a n D(a) n . Thanks to this, all formulas simplify remarkably. Our universe is of course not of the EdS form, because of the cosmological constant. But it is tempting to extend the results obtained in EdS to the ones in ΛCDM universe by replacing in the EdS formulas the EdS growth factor with the growth factor in ΛCDM. In this we obtain HereP 22 andP 13 are the time-independent one loop contributions given bỹ (r 2 − 2rx + 1) 2 P 11,l (kr, Λ)P 11,l (k √ r 2 − 2rx + 1, Λ) , P 13 (k) = k 3 1008π 2 P 11,l (k, Λ) (82) ∞ 0 dr 3 r 3 r 2 − 1 3 7r 2 + 2 log r + 1 1 − r − 42r 4 + 100r 2 + 12 r 2 − 158 P 11,l (kr, Λ) .
The counterterm contribution is given bỹ where the time dependence of c 2 comb can be inferred in perturbation theory from (30) to be Notice that, within this approximation, the k dependence and the time dependence of c 2 comb δ (1) l is the same as the one of the source of δ (3) l in the high k limit. This approximation is quite a good numerical approximation, and here below in Fig. 7 we present results of comparisons for P 22 and P 13 , where we see that the disagreement is at percent level. This ratio is plotted in Fig. 8. Notice that at k ∼ 0.24 Mpc −1 where the non-linear corrections are of order of a few ten percents, a percent error in the loop calculation leads to an error that is dangerously close to one percent. Therefore, for percent precision, using this approximate treatment at the high k's that we can reach with the EFT corresponds to pushing the boundaries of safety. It is further important to stress that this is not a parametrically good approximation in δ l , as is our loop expansion. As noted in [29,30,31], this is an expansion in the smallness of the ratio Ω m (a)/(∂ log D/∂ log a) 2 . A good fit is (∂ log D/∂ log a) 2 ∼ Ω(a) 1.2 , explaining the possibility to make this approximation. Luckily, we find it not so hard to implement directly the correct perturbation theory, which is the one we present in this paper.  This being equal to one would justify the approximations done in this appendix. We see that at late times the difference is quite large. Since this is a correction to the one-loop term, this error can be acceptable for a one-loop calculation. On the contrary the approximation can become more harmful if one goes to higher loops. The actual result on the power spectrum of the approximation is even better than what shown in the plot, as more of the clustering happens before dark energy domination.

D Measuring in N -body simulations D.1 Efficient calculation
We break the calculation into three parts: the calculation of primary fields, the calculation of secondary fields, and the calculation of correlations. Primary fields are dependent solely upon the positions or velocities of the simulation particles, and secondary fields are dependent upon the primary fields, and we care about correlations between particular secondary fields.
The immediate goal is to calculate the following expressions: based upon the two-point correlation functions: P IJ (r). These two point correlation functions can be seen as the expectation value of the product of field I with field J, but with I evaluated at all points a distance r from all points r : In the quantities being calculated the overall normalization cancels. The cost of evaluating each of these fields at all positions is prohibitive, so we approximate these correlations by measuring some large N number of pairs of points (pseudo)-randomly chosen but at fixed r: where N set (r) is a set of N pairs of points randomly selected from the space to be separated by a distance r.
The following are the specific correlation functions used: P δΘ (r) = δ l ( r + r )Θ l ( r ) Ω(r ,r),r P ΘΘ (r) = Θ l ( r + r )Θ l ( r ) Ω(r ,r),r (93) These are functions of the gravitational short-mode field A s , the over density δ, and the velocity divergence Θ, and relevant spatial derivatives. These fields are to be calculated from the observed (or simulated) positions and velocities of point-sources at a fixed moment in time (redshift=0 initially). We smear the observation of the position of these sources with a gaussian, introducing a (soft) ultraviolet cutoff.
To calculate c 2 s and c 2 v in terms of A s , δ, Θ, ∂ 2 δ, ∂ 2 Θ we use a number of secondary fields. While it would be possible to numerically estimate the relative necessary spatial derivatives, it is simple enough to explicitly carry the operations out analytically, and treat them as independent secondary fields, avoiding the introduction of specious numerical error.
have presented these fields, in detail, to emphasize that all such secondary fields rely on a relatively small number of primary fields, and it is only the primary fields which need concern themselves with the explicit number of particles down sampled from the simulation.

D.2 Stability of Measurement and Region Selection
The fluid parameter of interest c 2 comb is calculated as a ratio of polynomial functions of correlations. In fig. 9 we plot the numerator and denominator of these ratios for Λ = 1/3 (h/M pc) and Λ = 1/6 (h/M pc). It is worth noting strong confirmation of the effective field theory description is how well the numerator and denominator of the c 2 comb ratio tracks each other. With these sorts of statistical measurements, one should be careful of small fluctuations causing misleading signal near zero over zero regions. For the measurements described in the paper we select a region where the denominator is 3σ above zero.

E Renormalizing at finite Λ
The procedure as outlined in the main text involves renormalization for Λ = ∞ at some chosen k ren by fitting P 1−loop δδ for c 2 comb (Λ) against observation at that particular k ren . The power spectrum at all other k become predictions of the EFT. We could imagine to perform the same procedure at finite Λ. In this case, however, higher derivative terms suppressed by powers of k/Λ should be included. These terms do indeed vanish as Λ → ∞, but at finite Λ and k they are not negligible. In fact, as shown in Fig. 10, without the addition of higher derivative terms, the power spectrum deviates from the Λ = ∞ as k → Λ. One can indeed check that the values of c 2 comb that is obtained by fitting in this way is off with respect to the correct value as obtained from running down to finite Λ the value of c 2 comb at Λ = ∞ by about 15%, depending on the cutoff used, see Fig. 11. Instead, by allowing for higher derivative terms, the correct value of c 2 comb is derived. Figure 10: Prediction of the non-linear power spectrum without the addition of higher derivative terms, as we send Λ → ∞, normalized to the non-linear power spectrum. We see that if we keep Λ finite, non-included higher derivative terms that scale as powers of k/Λ are important. Indeed the results improves as Λ = ∞, which is the correct procedure. Figure 11: Ratio of the value of c 2 comb (Λ, k ren ) as obtained from running from c 2 comb (Λ = ∞, k ren ) versus the one obtained by fitting directly the result of the EFT at finite Λ without the inclusion of higher derivative terms in k/Λ. At low k ren , the error is particularly pronounced for Λ = 1/6, as in that case k ren /Λ is not very small. Inclusion of higher derivative terms reduces the mismatch to few percent.