Systematic Renormalization of the Effective Theory of Large Scale Structure

A perturbative description of Large Scale Structure is a cornerstone of our understanding of the observed distribution of matter in the universe. Renormalization is an essential and defining step to make this description physical and predictive. Here we introduce a systematic renormalization procedure, which neatly associates counterterms to the UV-sensitive diagrams order by order, as it is commonly done in quantum field theory. As a concrete example, we renormalize the one-loop power spectrum and bispectrum of both density and velocity. In addition, we present a series of results that are valid to all orders in perturbation theory. First, we show that while systematic renormalization requires temporally non-local counterterms, in practice one can use an equivalent basis made of local operators. We give an explicit prescription to generate all counterterms allowed by the symmetries. Second, we present a formal proof of the well-known general argument that the contribution of short distance perturbations to large scale density contrast $\delta$ and momentum density $\mathbf\pi(\mathbf k)$ scale as $k^2$ and $k$, respectively. Third, we demonstrate that the common practice of introducing counterterms only in the Euler equation when one is interested in correlators of $ \delta$ is indeed valid to all orders.


Introduction
A robust and accurate understanding of the gravitational clustering of Dark Matter is one of the main goals of cosmology. The smallness of the smoothed density and velocity on large scales suggests that some appropriate perturbation theory should converge to the right answer in this regime. It has recently become clear [1,2,3] that a consistent treatment of this problem requires the addition of an infinite hierarchy of effective interactions to the fluid equations of Standard Perturbation Theory (SPT) [4]. 1 Although the lowest order effective interactions have been discussed in details (in the density [2,3,17,18,19] and velocity [20] power spectrum, and in the density bispectrum [21,22], including non-Gaussian initial conditions [23]), a complete prescription to generate all the relevant operators to all orders in perturbation theory is still missing. In addition, the current renormalization procedure -a procedure that makes the perturbative expansion physically meaningful -has been carried out only at the lowest order for a limited set of interesting observables, and not in a fully systematic way.
In this paper, we introduce a renormalization procedure that works systematically to all orders in perturbation theory and give an explicit prescription to construct all effective counterterms allowed by the symmetries. Here we summarize our main findings: • Systematic renormalization In SPT one solves the continuity and Euler equations order by order:δ The nonlinear terms relate short and long wavelength modes to each other and lead to "loop" diagrams in the perturbative calculation of correlation functions. To be meaningful, loops need to be regularized and new interactions or "counterterms" must be included. With the exception of [20], all treatments so far have considered counterterms only in the Euler equation. Moreover, the SPT kernels (and the associated diagrams) combine several qualitatively different contributions together, making it practically impossible to separate embeddings of lower order loop diagrams into a higher order diagram. Because of this, the connection among various terms in the loop expansion, e.g. 1-loop 2-point correlator and 1-loop 3-point correlator, is less clear than in the common treatment of for example quantum field theories. Here, we introduce a systematic renormalization that avoid these drawbacks and is analogous to the one commonly adopted in quantum filed theories (including the Feynman diagramatics). For additional clarity, the procedure is illustrated using the explicit examples of 1-loop power spectrum and 1-loop bispectrum.
Our systematic renormalization makes explicit the implications of the long memory of the short-scale modes, of the equivalence principle and neatly organizes the generation of new counterterms. Nevertheless, it has its own drawbacks. First, it requires the introduction of counterterms in the continuity equation. These are indeed necessary to discuss velocity correlators, but they are degenerate with those in the Euler equation if one is only interested in correlators of δ. Second, systematic renormalization requires that all counterterms are non-local in time. A practical remedy to these shortcomings is presented in section 6.
• Mass-weighted velocity/momentum description To renormalize the one and two-loop power spectrum and one-loop bispectrum, it is sufficient to include counterterms to the Euler equation [2,20,18,21,22]. It is well-known that this approach works to all orders and for all statistics of δ, and is equivalent to a particular definition of velocity in terms of momentum current. 2 In section 6.1 we summarize and sharpen those arguments.
• Double softness It is well-known (see e.g. [24]) that the conservation of momentum and locality implies that interactions among modes of characteristic wave-number q can generate a longer wavelength mode k ≪ q suppressed at least as δ(k) ∝ k 2 , namely the density is "double soft". Analogously, the momentum is single soft, |π| ∝ k. After reviewing a general argument for this property due to Peebles, we prove it directly from the equations of motion in section 5 (details are in appendix B). Double softness can be seen and implemented most straghtforwardly using mass-weighted velocity or momentum. On the other hand, in terms of the more general velocity used in [20] and in the systematic renormalization, double softness implies a highly non-trivial relation among counterterms to all orders.
In this paper, we have extensively used the same nomenclature as in quantum and statistical field theory to stress the similarity with the problem at hand and to make the discussion more intuitive. When needed, we have adapted and refined some of the definitions. For the convenience of the reader, a glossary of the technical terms is provided in appendix E.
The rest of the paper is organized as follows. In section 2, we review standard perturbation theory in a compact "vector" form and introduce diagrammatic rules, which we elucidate for the case of the tree-level bispectrum. In section 3, we discuss how interactions give rise to loops and the separation between 1-particle irreducible (1PI) and 1-particle reducible (1PR) diagrams and stochastic and non-stochastic counterterms. The one-loop power spectrum is used as an example. In section 4, we present our systematic renormalization procedure and discuss the implication of the long memory of the short-scale modes and of the equivalence principle. The one-loop bispectrum is used as a concrete example of our general remarks. Section 5 is devoted to the discussion of double softness. Finally, in section 6, we give the explicit form of all allowed counterterms and put the density-only description on a firmer footing. Some technical details are left to the appendices.

Preliminaries
At very large scales dark matter can be treated effectively as a pressure-less fluid of density contrast δ = (ρ −ρ)/ρ, and velocity v i , whose dynamics (in a matter dominated model of expansion rate H) are governed by the continuity and Euler equations: where the gravitational potential satisfies ∇ 2 φ = 3 2 H 2 δ. Starting from zero vorticity, i.e. ∇×v = 0, it remains zero under this evolution. Hence, we can take the divergence of the latter equation and use the velocity divergence as the independent variable θ = ∂ i v i : (We will comment on the generation of vorticity in section 4.3.) At shorter scales the above description is incomplete and higher derivative corrections must be included. In practice, the need for such corrections (counterterms) can be seen from the dependence of the nonlinear solutions on the choice of the ultraviolet cutoff. The coefficients of the counterterms depend on the cutoff in such a way that the final result is cutoff-independent, though it does depend on the unknown details of the ultraviolet physics. That is, the counterterms capture the contribution of short distance physics by introducing new effective interactions among long-wavelength modes. Their coefficients are to be fixed either by a first-principle calculation or, more realistically, using simulations and observations. In this paper, we will derive the general form of these counterterms. For this purpose, it is useful to develop a systematic approach to renormalization-the procedure of removing cutoff dependence-as in conventional relativistic field theories. The main difference is that here we are dealing with an initial value problem, whose implications becomes clear in the sequel. In particular, it helps to have a more symmetric treatment of the primary variables δ and θ, and a more explicit diagrammatic representation of perturbation theory. Throughout the paper we work in an Einstein-de Sitter (EdS) universe but to a good approximation the conclusions remain unchagned in ΛCDM cosmology. Following [25,26] we introduce a doublet field 3) The equations of motion can then be written (in momentum space) as and non-vanishing components of the symmetrized vertex functions γ (s) are given by The vertices can be read from (2.1): (2.9)

Perturbation theory and diagrammatic rules
Given some initial condition Ψ a (k, η 0 ) = (δ 1 (k), δ 1 (k)), the above system can be solved perturbatively to find Ψ a (k, τ ). At first order, we have where repeated indices are summed over and D + ab is the growing solution of the linearized system. Note that it is k-independent. D + ab (η) is related to the conventional growth factor through The higher order solutions can be expressed, recursively, by combining two lower order ones via the interactions vertices γ abc : where the retarded Green's function, or propagator, is given by for η > η ′ , and 0 otherwise. The growth function D + ab (η, η 0 ) coincides with the growing part of g ab (η, η 0 ) and, since η 0 is practically −∞, it can be approximated by g ab (η, η 0 ). This perturbative expansion is diagrammatically illustrated by Feynman diagrams as (2.14) Every propagator is represented by a line, and the cubic interactions by vertices. The arrows indicate the flow of time, and clearly the first index of γ abc , which corresponds to the out-going line, has no symmetry property with the last two which refer to the in-going lines. A typical diagram contributing to Ψ (4) is where here and in what follows the indices are often omitted to avoid clutter. Note also that the momentum integrals in (2.12) must be understood as running up to some ultraviolet cutoff Λ beyond which the effective theory is not applicable anymore. Renormalization ensures that the final results do not depend on the choice of Λ.
Using the above Feynman rules, the following recipe can be used to calculate correlation functions perturbatively. An external leg Ψ (n) a is given by the sum over all tree diagrams constructed by successive bifurcations through trilinear interaction vertices to get to n initial fields δ 1 (k), as in (2.15). Any trilinear vertex merges a Ψ (n) and a Ψ (m) to result in an (n + m)-th order Ψ (n+m) field. These perturbative expansions for the external fields are then correlated using the initial statistics of δ 1 (k). We assume Gaussian initial condition (2.16) which reduces the correlation functions of several Ψ fields to a sum over all possible pairings of the initial fields that appear in their perturbative expansion. The correlation of initial fields is called contraction and is reperesented by a dot. If all initial momenta in a given diagram are fixed in terms of the external ones, the diagram is called tree-level. One example is Otherwise, if there are undetermined initial momenta it is a loop diagram. There is one undetermined momentum for each loop and it must be integrated over. The tree-level contribution to an N -point correlation function is of order 2(N − 1) in δ 1 , and an n-loop contribution is of order 2(N + n − 1). It is common to label various diagrams by the orders of their external legs in perturbation theory. For instance, the 1-loop contribution to the power spectrum where Ψ (1) is correlated with Ψ (3) is often called P 13 .

Example: tree-level bispectrum
Let us conclude by illustrating the above formalism in the example of tree-level 3-point correlation function of density contrast In the conventional SPT formalism this is given in terms of the F n kernels: where ). (2.20) In the doublet formalism, the relevant diagram is where we demonstrated the SPT diagram on the right for comparison. Following the Feynman rules, one obtains Using the above expressions for g ab and γ abc (and, of course, including other permutations), it is easy to show that this agrees with (2.19).

Loops
Loops are encountered in the perturbative calculation of correlators when the initial momenta are not fixed in terms of external momenta, and therefore, are integrated over. The ultraviolet part of these integrals include short wavelength modes which are not expected to be well described by the effective field theory. The integrals are often divergent. One should therefore regulate the integrals (say by sharply cutting them off at momentum Λ), and use the cutoff dependence to introduce new effective interactions (counterterms) that capture the contribution of the short scale physics to the low momentum observables as an expansion in k/Λ. This procedure provides a systematic way of inferring the form of higher derivative corrections to be added to the long-wavelength effective theory, and in most cases, it covers all possibilities that are compatible with the symmetries of the system under consideration. 3 Hence, understanding better the properties of loop diagrams is the key to identify the general structure of higher derivative corrections. This is the subject of this section.

Reducible and irreducible diagrams
Consider the P 13 contribution to the 2-point correlation function Suppressing time-integrations and indices the diagram would have the following momentum structure To renormalize this loop integral a linear counterterm must be added to the equations of motion. Its leading contribution to the 2-point correlation function combines with the original diagram to make the full result cutoff-independent. Now suppose the same loop diagram of (3.1) is inserted into a higher order graph, e.g. should renormalize (3.4). We will see that such a choice of the counterterm is indeed possible, and in particular, the diagram (3.5) does renormalize diagram (3.4). Hence, embedding lower order loop diagrams in higher order ones does not lead to genuinely new corrections (counterterms). We call these loop diagrams reducible, or in analogy with particle physics diagrammatics, 1-Particle Reducible (1PR). Their characteristic property is that there exist an internal line l that is not immediately connected to any external line through a contraction such that cutting l divides the diagram into two disconnected pieces. All other diagrams are called 1-Particle Irreducible (1PI) diagrams, (3.1) is an example.
Since the counterterms are to be added to the equations of motion as new vertices, it makes sense to work not with the full diagrams representing the correlation functions but directly with the time-evolution diagrams. These show the evolution of a set of initial fields into one higher order final field at a later time. The loops correspond to integrating over some of the initial momenta, and hence, to allowing them to become "hard", i.e. of the order of the ultraviolet cutoff Λ. We represent these hard momenta by thick lines. Moreover, to exclude 1PR diagrams one would simply amputate all other ("soft") lines. That is, the relevant diagrams are made of (i) several hard lines extending all the way to the initial condition δ 1 , (ii) one soft out-going line, and (iii) any number of soft in-going lines which are immediately connected to the hard lines. An example is (3.6) The presence of high-momentum hard lines makes these diagrams cutoff dependent. They are renormalized by adding new vertices to the effective theory that have the same structure of the soft lines, and the same cutoff dependence in their coefficients: Note that the contribution of the new vertex to the final field always comes with the propagator of the out-going soft field which is common between (3.6) and (3.7). Since we are interested in the structure of the corrections to the equations of motion and not the resulting final field, we will also amputate the final propagator.

Stochastic vs. non-stochastic corrections
There is another important classification of counterterms that is unique to the case of initial value problem. The high-momentum initial fields were seen to result from integration over unfixed initial momenta. If in a given time-evolution diagram all of these can be paired and contracted with each other, the associated counterterm is called a non-stochastic one, e.g. non-stochastic diagram: The counterterms of these diagram capture the fact that short scales are affected by large scales.
Otherwise, if at least one pair of the initial hard fields do not have exactly opposite momenta, it corresponds to the probability that short scale fluctuations average into a long-wavelength mode, e.g.
The associated counterterms are called stochastic. These are by nature indeterministic and hence contain a stochastic field (usually denoted by J(x, η)) about which we can only make statistical statements.
One can contract and integrate over the momenta of all initial hard fields that are paired in a time-evolution diagram (which include all of the hard lines in a non-stochastic diagram). This is because the loop is independent of how the ingoing soft lines are embedded in a more complicated diagram. For instance, the hard modes δ 1 (q) and δ 1 (−q) in fig. 3.8 can be contracted regardless of whether we consider diagram of (3.1) or (3.4). After doing so these diagrams (and the associated counter-terms) incorporate the average response of the short modes to the initial long modes. 4 On the other hand, the unpaired initial hard modes will be contracted by unpaired modes coming from other external lines as in Our primary focus will be on non-stochastic counterterms since in our universe the stochastic ones are expected to be sub-dominant at long wavelengths (see below).

Example: One-loop power spectrum
As an explicit example let us consider the 1-loop non-stochastic contribution to the power spectrum (3.1). This corresponds to whose amplitude is given by (recall that amputation removes the final propagator) Here the time variable is defined as a = e η , g a 2 1 (η 1 , η 0 ) = g a 2 2 (η 1 , η 0 ) is the same as the growth factor, and the summation of repeated indices is assumed. Using the explicit expressions for the propagator and vertices, and restricting to q ≫ k, we obtain the following expression for the ultraviolet part of the loop diagram (3.14) The velocity dispersion, which is defined as is the cutoff dependent coefficient that we aim to renormalize by the addition of a counterterm proportional to k 2 . This will be done in the next section.

Renormalization and counterterms
In this section, we first explain how to read off the counterterm associated to each amputated loop graph. Then, we discuss various properties of these counterterms, give their general structure, and conclude with an explicit application of these ideas to the 1-loop bispectrum.
The amplitude of a generic amputated diagram (such as (3.6)) with the final vertex inserted at time η, outgoing momentum and index (k, a), and a set of soft ingoing lines characterized by {k n , b n , η ′ n }, and hard ingoing lines characterized by {q i } (recall that hard lines are all initial fields In the non-stochastic case the hard lines {q i } can be paired and contracted to give P lin (q i ) and integration over their momenta results in The C ab of last section was an example of this with only one ingoing soft line. The C kernel is a cutoff dependent quantity. This diagram is renormalized by the addition of a counterterm of the following form to the equation of motion (c.f. To cancel the cutoff dependence of the loop diagrams, the dependence of c a... on {k n } has to match that of C U V a··· , namely the cutoff dependent part of C a··· when it is expanded in powers of k/q. Let us analyze this more closely. In the ultraviolet regime the momentum dependence of C a··· can always be brought via Taylor expansion into a sum of terms, each of which factorizes into a product of a function of {q i } times a function of {k n }. We are interested in the functional dependence on {k n }. Consider a vertex at η n where the n th soft line is attached to the diagram. From the equations of motion (2.1) and (2.2), the leading order vertices of the theory are Thus, the soft line can realize one of the following six possibilities: These are the only ways in which c a··· may depend on the n th field, up to possible multiplication by positive factors of k i n , arising from vertices appearing at later times in the diagram. The contribution from these later vertices is guaranteed to be analytic in k i n because (i) by definition of an amputated graph, a soft line n attaches immediately to a hard line, say with momentum q n , thus (ii) the later vertices depend on k n only through q n + k n , and (iii) since k n ≪ q n one can Taylor expand in k i n . Moreover, c a··· has to have a cutoff dependent coefficient with exactly the same time-dependence and overall amplitude as the one arising from the integration over {q i } in C U V a··· . However, there can be a cutoff independent piece with different time-dependence. This counterterm by construction renormalizes all possible embeddings of the original amputated diagram into any higher order diagram. Next, we will discuss the manifestation of characteristic properties of the cosmic fluid in the form of these new, and yet rather abstract, corrections.

Long memory effect and non-locality in time
Counterterms capture the effect of the short wavelength modes on the long wavelength physics. They result from replacing hard loop diagrams with new effective interactions among soft lines (i.e. integrating out the short-distance physics). The result of this procedure is guaranteed to be spatially local since the short modes can probe the long modes only in a region of size Λ −1 over which the long modes can be accurately Taylor expanded in powers of k/Λ. However, before reaching the virialization scale, the short modes evolve with approximately the same time scale H −1 as the long modes. Hence, we should expect the counterterms to be generically non-local in time as a short mode can be influenced by long wavelength modes at some earlier time and respond at a much later time. The hard lines in a loop diagram can therefore be collapsed to a point spatially but not temporally [15,18].
This expectation is manifest in our general formula (4.3) and in the explicit 1-loop example. In order to systematically renormalize the amputated diagram of (3.12) one needs to add a temporally nonlocal quadratic vertex, corresponding to a linear counterterm in the equations of motion: where we omitted the standard cubic vertices. Here k 2 c ab has a cutoff dependent piece that (up to a minus sign) exactly matches C U V ab in (3.14) so that the sum of the two cancels, up to a regularization independent piece. As will be seen in section 4.5.1 this counter-term automatically renormalizes the 1PR contribution (3.4) to the bispectrum. On the other hand, if we were only interested in renormalizing the P 13 diagram (3.1), this vertex could be replaced by a local one (the conventional speed of sound and viscosity counterterms) As we will see in section 6 this reduction is always possible in perturbation theory. However, it is incompatible with systematic renormalization. For example, using the local-in-time counterterm in (4.7), on finds that the 1PR diagram (3.4) requires new counterterms [21].

Large scale flow and the Equivalence Principle
The origianl system (2.1) is invariant under the addition of a uniform bulk flow [27,28] x where n is an arbitrary time-dependent vector. This ensures that a very long-wavelength perturbation, whose main effect on much shorter scales is to produce a uniform acceleration, is locally unobservable, in accord with the Equivalence Principle. 5 Local observables can depend on the tidal field ∂ i ∂ j φ, the shear ∂ i v j , their spatial derivatives, and their convective time-derivatives The counterterms which characterize the response of the short scale physics to long-wavelength modes should also follow the same rules. The short modes can only respond to the locally measurable quantities. However, when the counterterms are non-local in time as in (4.6), they depend on these local observables along the large scale flow. The same considerations imply that the time integral must be taken along the fluid trajectory, which is implicitly defined by In the following, we often drop the last two arguments of x fl when no ambiguity arises. For concreteness, let us focus on the counterterm in (4.6). To be invariant under (4.9), it must be dressed to read [18] The significance of this replacement x → x fl soon becomes clear in the example of 1-loop bispectrum in section 4.5.2. Perturbatively expanding x fl around x, we get where all omitted space arguments are now x. This replacement essentially connects an infinite set of counterterms (and hence an infinite set of amputated diagrams) to one another. These are all diagrams in which the original diagram (3.8) is being displaced by the bulk motion produced by any number of ingoing soft lines attached at later times than the original attachment of More generally, whenever a soft line realizes the v i ∂ i part of the vertices in (4.4), it can be attributed to the shift from x to x fl of a lower order diagram with the soft line removed. These diagrams have the following characteristic property in perturbation theory: In the presence of a long-wavelength mode with constant amplitude δ 1 (k L ) and vanishing k L , they diverge as k L /k 2 L , namely, proportionally to the bulk velocity induced by the long mode. This will be referred to as infra-red (IR) singularity. The UV part of these diagrams are fully fixed in terms of that of the original (unshifted) diagram. The counter-terms associated to these diagrams too are fully fixed in terms of that of the original diagram. Therefore, they do not introduce any new parameters into the effective theory. Eliminating these diagrams by shifting x → x fl ensures that the new higher order counterterms needed to renormalize the remaining higher order 1PI diagrams depend only on locally observable quantities, again measured along the dark matter flow x fl . In other words there is no IR singularity, except in the argument x fl of the fields. 6

Soft outgoing momentum and vorticity
Another property of loop diagrams to be discussed is their softness in the outgoing momentum. In perturbation theory, it can be easily verified that all vertices are at least single soft: In any amputated loop diagram, the latest vertex is one in which two hard lines combine to give an outgoing soft line of momentum k (see e.g. (3.6)). Because of (4.15) the UV part of this diagram will contain at least one factor of k, and so do the corresponding counterterms. Rotational invariance implies that the full amputated diagram must be of second order in soft momenta since there is no singularity (i.e. negative power) in the ingoing soft momenta after taking care of the shift terms as explained in the previous section. Thus expressed in terms of Ψ a , the counterterms must start from second order in derivatives, one of them being an overall derivative because of the property (4.15) of the final vertex. An example is The ultraviolet part of C ab (η; k, η), calculated above agrees with this expectation. The above softness property of the counter-terms has a natural explanation based on the structure of fluid equations. The continuity equation relates the change in mass density to the divergence of the mass current (momentum density) 7 where π ≡ (1 + δ)v. (4.17) However, this relation between π and v depends on the definition of velocity. Fixing velocity at some scale Λ to be given by [π] Λ /[1 + δ] Λ and then coarse-graining to a smaller cutoff Λ ′ results in a velocity field that differs from Hence, when working with velocity field as the primary variable one must add a locally observable vector j c as a counterterm to the mass current. The absence of IR singularity ensures that a locally observable vector must be at least of first order in derivatives (in the same sense that (4.16) is of second order in derivatives). Upon taking the divergence, this leads to a counterterm ∇ · j c in the continuity equation that is second order in derivative counting, and because of the overall derivative it contains at least one factor of the outgoing momentum. Similarly, integrating out the short wavelength modes can introduce an effective force f c in the Euler equation. This is again a locally observable vector that incorporates the counterterms. After taking the divergence to obtain an equation for θ it leads to a counterterm ∇ · (f c /(1 + δ)) with an overall derivative.
It is worth noting that since the original system of equations (2.1) do not generate any vorticity, to cancel the cutoff dependence of just SPT loops, f c /(1 + δ) can be taken to be curl-free. In other words, the counterterms to the θ equation c a=2,··· are not only a total divergence but they are the Laplacian of locally observable scalars. We will verify this expectation in the explicit calculation of quadratic counterterms (see appendix A.2). Of course, the finite pieces of the counterterms do not have to obey this rule, and in principle the short scale dynamics can generate vorticity (see section 6.1 and appendix B.3 for comments on the non-uniqueness of the definition of velocity).
The property (4.15) and the fact that all counterterms have one overall derivative ensures that if the outgoing momentum k is much smaller than the soft ingoing momenta {k i } the solution for δ(k) and θ(k) is going to be single-soft in k. In fact, as we will discuss in section 5, momentum conservation implies that δ must not only be single-soft but rather double-soft in k in the limit k ≪ k i . As a result, only a subset of locally observable forces -those that respect this propertyare allowed to be added as counterterms to the system.

The list of counterterms
Using the above rules we can write down an over-complete list of counterterms which can be added to the equations of motion. The softness in the outgoing line implies that they must be of the form ∂ i c i . c i encodes the memory of all incoming soft lines that are attached to the hard lines of the amputated diagrams (3.6). Therefore, each field in c i (x, η) is integrated from η 0 to η against a general time-dependent memory function. Any attached ingoing line is either a locally measurable quantity made of the soft mode such as ∂ i ∂ j φ or ∂ i v j (and their convective time derivatives (4.10)), or it is a velocity, which is not a local observable. As argued above the latter is necessary to perturbatively dress the time integrals at constant x to those along the fluid flow. In the former case, since both ∂ i ∂ j φ and ∂ i v j have two indices, rotational invariance implies that to construct a vector c i one needs to have one additional derivative.
We proceed as follows. For simplicity, suppose there is zero vorticity; this allows us to write v = ∇φ v , and introduce Next, as an intermediate step we construct n th order two-index tensors c ij made of various contractions of the local observables Π ij a ≡ ∂ i ∂ j φ a . Using matrix notation, the first few orders are where 1 ij = δ ij , and as an example, Each Π a is evaluated at a different time which is integrated against a temporal kernel until the final time. An over-complete (see below) list of vectors c i consists of all possible ways of acting by ∂ j on any of Π a appearing in c ij . Finally, one takes the divergence ∂ i c i . One second-order example is where the derivative in ∂ j Π a 1 should, in principle, be with respect to x fl . However, the difference between this and the derivative with respect to x can be expressed in terms of higher order terms that are included in our list -in other words So we take it to be ∂/∂x j in what follows. Note also that the arbitrary time-kernel makes the use of convective time-derivatives redundant. Unfortunately the above procedure generates an over-complete list of counter-terms, as it allows us to introduce arbitrary observable vectors f c and j c as counter-terms. On the other hand, as we will argue in section 5, only a subset of these vectors -those which ensure double softness of δ when k out ≪ {k i } in -are permitted by momentum conservation. An allowed list can be obtained by restricting the kernels to those which arise from perturbation theory. This elects certain combination of these vectors that are necessary to cancel the cutoff dependence of the loops. These combinations inherit the double softness property of the leading order perturbation theory which is discussed in section 5.2. However, in this way all amputated diagrams with fixed number of ingoing soft and hard lines combine to give only one new counter-term whose coefficient can be allowed to have a cutoff independent piece. While this is sufficient to cancel the cutoff dependence of the loops, in order to produce a complete basis that can give an effective description for all possible microscopic scenarios one needs to proceed to diagrams with increasingly large number of hard modes.
We will next verify the adequacy of this formalism in cancelling the cutoff dependence in the explicit example of the 1-loop bispectrum. Several properties of the counterterms become apparent in this example. In section 6, we show how this rather cumbersome and nonlocal formulation can be reduced to a much more tractable and local one in any practical calculation by giving up systematic renormalization. Moreover, using the results of section 5, we will provide an alternative list of local counterterms that is complete but not over-complete.

Example: One-loop bispectrum
In this section we discuss the renormalization of the non-stochastic part of the 1-loop bispectrum which is conventionally called B 114 . The 1PI and 1PR contributions are evaluated and discussed separately.

1PR diagrams
As expected from the general argument at the beginning of this section, it is easily seen that the 1PR diagrams are automatically renormalized by the insertion of the linear counterterm in (4.6) in a tree-level diagram. The diagram corresponds to taking the expectation value We denoted the UV part of this loop diagram by C U V and it will be renormalized by adding the contribution of the diagram where the cross represents the counterterm in (4.6). The expression for this diagram is given by an identical expression as (4.25) with C → k 2 c, which have equal and opposite cutoff dependences. The other 1PR diagram has an expression (4.28) Similarly, the UV part of this diagram C U V is renormalized by the diagram containing the same counterterm in place of the loop: which has the same expression as (4.28) with C → k 2 c. It is worth emphasizing here that the temporal non-locality of counterterm c a a ′ (η; η ′ , k) is crucial for it to cancel the cutoff dependences of all 1PR loop integrals. This accommodates the possibility that the response of short modes should strongly depend on the time-dependence of the long modes. If only local in time counterterm are used, the UV divergences in 1PR diagrams are not generically canceled by the lower order counterterms, spoiling the systematic character of our renormalization procedure. This is expected and in agreement with the general argument presented in section 4.1 about the long memory effect.

1PI diagrams: shift terms
As discussed in section 4.2, some of the 1PI diagrams contributing to B 114 , whose Ψ play the role of shifting the temporally nonlocal lower order counterterms to the fluid position (here 2 indicates the velocity field). It is implicit in the above diagrams that the lines marked by 2 realize v in v · ∇Ψ a vertices, and that η 2 > η 1 . These correspond to contributions that diverge as 1/k 2 when k 2 → 0. We expect them to be renormalized by including the first shift term in (4.13) in the counterterm for the 1-loop power spectrum (4.12): That this cancellation of IR-singularities must happen is almost obvious since any diagram that has a potential singularity (i.e. any diagram where the soft line attached to the loop is a v i ∂ i ) can be identified with a contribution to (4.32). An explicit check is provided in appendix A.1.

1PI diagrams: new counterterms
After the subtraction of shift terms from 1PI diagrams, what remains must be renormalized with the genuinely new higher order counterterms. These are made of locally observable operators integrated over the history of the short modes but evaluated at x (their displacement to x fl would correspond to higher order counterterms). At this level one must add a second order term to the equation of motion, corresponding to a cubic effective vertex. From (4.20), there are 18 possibilities for quadratic terms in c ij . Note that even if a 1 = a 2 the two soft legs of the amputated diagram can be embedded in different ways in a connected diagram and be integrated against different timekernels, therefore no symmetrization is needed. When taking derivatives to construct counterterms ∂ i c i , the number of operators proliferates but some are degenerate. It follows that, for fixed a 1,2 and η 1,2 , possible c i vectors are Thus there is a total of 15 independent, temporally nonlocal quadratic counterterms. It is shown in appendix A.2 that they are sufficient to fully renormalize the 1PI diagrams. In fact, as expected from the structure of the leading vertices (4.4) and the requirement of momentum conservation (see the end of section 4.4), only a subset of these counterterms are allowed by the symmetries. The correct prescription will be discussed in section 6. Also, because SPT does not generate vorticity (see subsection 4.3), an even smaller subset is needed to cancel UV-dependencies.
In summary, we have demonstrated that the effective theory of large scale structure can be systematically renormalized, at the expense of introducing a plethora of temporally nonlocal counterterms. After discussing the implications of locality and momentum conservation in the next section, in section 6 we explore a more practical approach.

Momentum conservation, locality and double softness
As alluded to in section 4.3 momentum conservation and locality of the short scale dynamics ensure that the time evolution of short wavelength fluctuations can only lead to the generation of a longer wavelength perturbation suppressed at least by k 2 out [24]. In the following we will first review the heuristic argument of Peebles, and next show that double softness follows from the equations of motion of the fluid system.

A general argument for double softness
First consider an almost homogeneous initial densityρ 0 with small short wavelength fluctuations of scale 1/q, but no long wavelength fluctuation (δ(k, η 0 ) = 0). Due to gravitational instability, these initial fluctuations collapse and form a clumpy distribution of matter ρ(x). We ask what is the typical size of a long wavelength fluctuation δ(k, η) after this process? Approximating the final distribution by a set of point particles of mass m n at location x n , we have whereρ is the mean final density. If the short-scale dynamics was turned off, the local mass distribution would be uniform; we denote it by a tilde:ρ(x) =ρ. There is a nonzero δ(k, η) because this would-be uniform distribution has clustered. in this hypothetical uniform universe, to keep track of the matter that ended up in each clump, we divide the space into (possibly overlapping) regions R n with densityρ n (x). All of the mass in this region falls into the n th clump: m n = Rn d 3 yρ n (x n + y).

(5.2)
If there is no overlapρ n (x) is the same asρ(x) (and hence uniform) within R n . But if two regions R n 1 and R n 2 overlap, then in the overlapping regionρ n 1 (x) +ρ n 2 (x) =ρ (and similarly for more overlaps). The typical size of each region is the extent to which mass elements are displaced from their initial position, namely, of order 1/q. Now we can write where d i n and Q ij n are, respectively, the first and second moments of the mass distribution in R n as measured from x n : Q ij n = Rn d 3 y y i y jρ (x n + y).

(5.4)
The expansion of the exponential is justified because of the smallness of the size of the region compared to 1/k. Comparison with (5.1) implies The second and higher order terms on the r.h.s. satisfy our general expectation about the dependence on k. As for the first term, the sum can be related to the center of mass position of each region: The difference (x CM n − x n ) can only be caused by momentum transfer among different regions. If initially there existed only modes of wavelength ∼ 1/q, there is no coherent momentum transfer over much longer distances, and the sum n m n (x CM n − x n ) over any patch of size much larger than 1/q is negligible. Hence, the sum in (5.6) is nonzero to the extent that the exponential varies over any such patch, and therefore it must be suppressed at least by one factor of k. This shows that all terms in (5.5) are suppressed by at least two powers of k. 8 In the absence of long wavelength initial perturbations the location of the clumps x n is effectively random at large scales, and the sums in (5.5) become identical to random-walks. Correlating two such contributions and averaging over the short modes results in This is the generic asymptotic behavior of stochastic contributions to the long-wavelength power.
In the presence of initial long wavelength perturbations the above argument can be modified as follows. Now the clumpy universe should not be compared with a homogeneous universe, but with one in which the long modes are kept but the short modes are absent. Hence,ρ(x) is smooth but nonuniform. All other steps remain the same, though one should bear in mind that x n and x CM n both depend on and can get largely displaced by an initial long mode δ 1 (k 1 ); the leading dependence which goes as δ 1 (k 1 )/k 1 cancels in the difference in (5.6). Nevertheless because of the remaining dependence of both x n and the shape of each region R n on the initial long modes, (5.5) has nonzero correlation with the initial δ 1 (k 1 ). The non-stochastic diagram (3.1) (equivalently (3.8)) is one perturbative analog of such contribution.
Finally, a similar argument implies that the momentum density π(k, η) = 1 ρ n m n v n e ik·xn (5.8) starts from order k.

Double softness in perturbation theory
The double softness property is not manifest in the formulation (2.1) in terms of δ and v. However, it becomes manifest if we eliminate v in favor of the momentum density π = (1 + δ)v, in terms of which the continuity equation becomes linear (4.17) and the Euler equation becomeṡ Writing shows that the non-linear terms in this equation all have an overall derivative. Moreover, taking the divergence of this equation and using ∂ i π i = −δ to obtain an equation for δ makes all interaction vertices explicitly second derivative. In appendix B it is shown that these equations evolve ingoing hard lines into a single soft answer for π(k out ) and a double soft answer for δ(k out ) when k out ≪ {k i } in . The same conclusion is shown to hold when there are ingoing soft lines with momenta of the order of k out in addition to the hard lines. It then becomes clear that the softness property is preserved if counterterms of the form ∂ j τ ij are added to equation (5.9), where τ ij is made of local observables. This ensures locality and momentum conservation.

A practical remedy
So far we learned that a systematic renormalization of the two field (δ, θ) formulation of cosmic fluid is possible if one uses nonlocal in time counterterms and add corrections to both continuity and Euler equations. However, there are three reasons to seek an alternative approach: I One may not be interested in the velocity statistics, and even if she is, there are various ways to define velocity field and the one used in section 4 might not agree with the velocity field obtained from an actual survey or simulation.
II More importantly, the momentum conservation constraint is not straightforwardly implemented.
III The nonlocal in time counterterms introduce additional complications, since instead of a single time-dependent coefficient one needs to consider a multi-variable function of insertion times for every new counterterm.
We first show in section 6.1 that one can use the momentum formulation to eliminate velocity. In this language, only the Euler equation receives corrections and not the continuity equation. Also, now there exists a straightforward prescription to generate all counterterms and only the ones allowed by the symmetries, including the conservation of total momentum. The observed velocity on the other hand must be treated similar to biased tracers; it can differ from the one calculated in any perturbative scheme by a set of locally observable counterterms whose coefficients must be empirically determined. This resolves the first two problems. Then in section 6.2 we show that all nonlocal in time counterterms can be reduced perturbatively to a set of local operators, addressing problem III. This allows us to give a complete and more manageable list of counterterms.

Formulation in terms of the mass-weighted velocity
If we use the density δ and momentum π as the primary variables, the continuity equation reduces to the linear equation (4.17), and does not require any counterterm (it is not modified by smoothing) [2,18]. If the correlators of δ (and therefore ofδ) are renormalized, then so are the correlators of ∇ · π and vice versa. The same can be said for φ and δ, given that the Poisson equation is also linear. The problem then reduces to solving the only non-linear equation, namely the vector Euler equation (5.9). If one is interested in the correlators for δ, it is convenient to decompose the momentum into a divergence (µ ≡ ∇ · π) and a vorticity (ν ≡ ∇ × π), as discussed in appendix B (see (B.1), (B.2) and (B.3), see also [31]) . The Fourier space version of these equations is collected in appendix D. Notice that all nonlinear terms on the right hand side of the Euler equation (5.9) are total derivatives of some tensor that is not proportional to δ ij and hence all contribute to the curl equation. In particular, all non-trivial solutions have ν, δ = 0. Using the continuity equation (4.17), one can get rid of µ and solve the remaining equations for δ and ν. When using momentum, one can therefore compute all correlators of δ without including any counterterm in the continuity equation. We will see now that this is still the case when using a specific definition of the velocity.
In Standard Perturbation Theory (SPT), where τ ij = 0, there exist a convenient change of variables for which the curl of the vector Euler equation is trivially satisfied even for non-trivial solutions δ = 0. From (5.9), it is clear that a necessary condition for this to happen is that the curl is not sourced by the term (1 + δ)∂ i φ. Upon dividing the whole equation (5.9) by (1 + δ), this term reduces to the gradient of a scalar which disappears after taking the curl. This suggest to use the variable v π ≡ π/(1 + δ), which, as it is well-known, indeed decouples the vorticity ∇ × v π from δ in SPT. We will refer to v π as "mass-weighted" velocity.
More explicitly, let us define the divergence and curl of the velocity (we raise and lower all indices with δ ij ) from which it follows ∂ i w i π = 0. The equations of motion then take the schematic form (see appendix D for the explicit form) ∂ τ δ + θ π = −αθ π δ + α w · w π δ , (6.2) ∂ τ θ π + Hθ π + 3 2 H 2 Ω m δ = −βθ π θ π + β w · w π θ π + β ww ij w i π w j π + ∂ · ∂τ 1 + δ (6.3) ∂ τ w π + H w π = γ w ij w j π θ π + α w w π · w π + γ ww ijl w l π w j π + ∂ × ∂τ 1 + δ , (6.4) We can now perturbatively solve for θ π and w π in terms of δ using (6.2) and (6.4). This yield, order by order, a single differential equation for δ. Let us see explicitly how this works to quadratic order. The vorticity is sourced by ∂ ×[∂τ /(1 + δ)] and therefore starts quadratic in δ perturbations. We can hence neglect it in the following formulae, keeping in mind that it will contribute to δ at the next order, O(δ 3 1 ). We now solve the continuity equation 5) and substitute it into the divergence of the Euler equation to obtain a nonlinear equation purely in terms of δ:δ The linear part of this equation is a second order differential equation which gives a retarded Green's function, using which one can perturbatively solve for the evolution of δ(τ, x) in terms of the initial condition δ 1 (x), without ever mentioning θ π . Obviously the resulting equation for δ now has interactions of all orders in δ. On the other hand, the original system (2.1) contains only cubic interactions and quadratic mixings between intermediate θ and δ lines. What the above procedure of solving for θ π in terms of δ means is that all intermediate θ lines are contracted to points, thereby generating arbitrarily high order interactions of δ. This formulation justifies a common approach in the EFT literature which is to use δ and θ as variables but to consider counterterms only in the Euler equation while keeping the continuity equation intact. It is equivalent to renormalizing the equation (6.6), as any conventional nonlinear theory of a single field δ, and to redefine the composite operator θ π to absorb the counterterm ∇ · j c of the continuity equation, so that the redefined θ π is given by the same relation (6.5) at all scales. However, doing so can in principle introduce vorticity because ∇ · j c is not necessarily Laplacian of a locally observable scalar. Therefore, the counterterms of the (vectorial) Euler equation needed to renormalize the SPT loops are no longer pure gradients. This was the case in [21], where the quadratic counterterms ∂ j (∂ i ∂ k φ∂ j ∂ k φ), which is not the pure gradient of a locally observable scalar, was needed to renormalize one-loop bispectrum of δ.
In summary, one can use the standard perturbation theory approach to calculate various correlation functions of δ. The counterterms that are allowed to be added to the vector Euler equation are now fixed by momentum conservation to be of the form where in order to do systematic renormalization τ ij must be allowed to depend on observable quantities along the flow, and hence be nonlocal in time.
As for the statistics of the observed velocity field, the difference with that of perturbation theory (say v π ) can be parameterized by a set of counterterms v obs = v π + v ct , (6.8) where v ct is a locally observable vector. It has a non-stochastic piece which by symmetry must start from first order in derivative counting, and a stochastic piece which according to the arguments of appendix B.3 is generically O(k 0 ), leading to O(k 2 ) contribution to the vorticity power spectrum. Hence, the treatment of velocity is similar to the biased tracers. In particular the most general v ct can be expressed locally in time [8].

Local formulation
It was argued in the last section that a generic counterterm can be constructed out of a tensor τ ij . The tensor is made of observable quantities integrated along the flow: where we switched to the conformal time τ . In perturbation theory, any observable can be solved in powers of the initial field δ 1 (x), in particular, where D(τ ) is the growth factor and we are using the accurate approximation of replacing a(τ ) → D(τ ) in going from matter dominance to ΛCDM. This expansion contains displacement terms to take into account the bulk flow induced by long wavelength initial perturbations. As such, some of the second and higher order terms in the sum become singular in the presence of a finite δ 1 (k 1 ) with vanishing momentum k 1 → 0. These terms are responsible for shifting the argument of the fields on the right hand side (r.h.s.) to their initial (Lagrangian) fluid position After shifting the arguments x → z there remains no infra-red singularity on the r.h.s. as expected on symmetry grounds. We denote the operators thus obtained Lagrangian operators. In terms of them [21,8] Π(x, τ ) = D(τ )Π (1) lgr (z) + D 2 (τ )Π (2) lgr (z) + · · · (6.12) Substituting this expansion in (6.9), makes it possible to perform all time integrations to obtain an expansion of various products of Lagrangian quantities with coefficients that generically depend on (and only on) the final time τ . The whole expression is acted on by two Eulerian derivatives: It was shown in [8] that there is a one-to-one correspondence between the Lagrangian quantities and the locally measurable (Eulerian) operators. Hence, we can re-express (6.13) in terms of those Eulerian operators. A basis for the latter can be obtained by combining various convective time derivatives (4.10) of ∂ i ∂ j φ. Starting from the convective derivatives can be combined in such a way that a set of operators with increasing order in perturbation theory is obtained. Suppressing the indices and the argument (x, τ ) these higher order operators can be computed recursively: where prime denotes convective derivative with respect to log D (= η in a matter dominated universe). At leading order these operators match those that appear in the perturbative expansion of ∂ i ∂ j φ:Π (n) (x, τ ) = D n (τ )Π (n) lgr (z) + · · · (6.16) At lowest derivative level an Eulerian basis for two-index tensors can be constructed as follows. First, a two index tensor is made out of a (matrix) product of a (possibly empty) sequence ofΠ (n) 's. For instance, δ ij ,Π lj are two examples. This tensor can be multiplied by an arbitrary scalar constructed out of traces of products ofΠ ij . The trace Tr[Π (n) ] with n > 1 must be excluded since it is expressible in terms of lower order operators. The first few orders are (in matrix notation): Note that the τ ij matrix must be symmetrized. To obtain a list of counterterms, one substitutes any such tensor into the τ ij of (6.7). Let us close with a few clarifications.
• This is a complete basis at lowest derivative level. In particular, all locally observable operators made of the velocity field are also included. For instance, in terms of a rescaled velocity potential Clearly, at higher orders ∂ i ∂ j φ and ∂ i ∂ j φ v will be insufficient for a local description of counterterms and higher order convective time derivatives will be necessary. However, given that at n th order we need up toΠ (n−1) , the first example appears at the quartic level.
• The tensorsΠ (n) (x, τ ), and hence the counterterms in (6.7), are easily calculable in momentum space by combining the F and G kernels of SPT. For instance, where δ n (k) is the n th order SPT solution, and and so on.
• One should bear in mind that although ∂ i ∂ j φ(x, τ ) is a total derivative, neither its convective time derivatives nor the z-space matrices Π (n) ij (z) are necessarily total spatial derivatives of any kind.
• Although, the basis (6.17) for tensors is non-redundant, taking the derivatives leads to degeneracies. For instance, at second order there are 5 tensors but only 3 indepedent counterterms.
• Any tensor τ ij leads to an infinite series of counterterms when the factor (1+δ) −1 is expanded in (6.7). The relative coefficients of these counterterms are fixed so as to insure double softness of δ when the ingoing fields are hard (a fact that is verified in appendix B using the momentum equations).

Conclusions
In this paper, we have introduced a systematic renormalization procedure for the perturbative description of Large Scale Structures, which is analogous to the one commonly used in quantum field theory. One main advantage is that new counterterms can be straightforwardly read off amputated time-evolution diagrams, without the complicated mixing with lower-order counterterms that plagues the standard treatment. In this approach, correlators of both δ and θ are renormalized at each order. We have also provided for the first time an explicit prescription for the allowed counterterms to all orders in perturbation theory. In passing, we have justified the commonly used mass-weighted volcity formulation in which counterterms appear exclusively in the Euler equation. 9 We have also presented a rigorous derivation of the double-softness of δ and ∇ · π in perturbation theory, namely the well-known fact that δ(k) ∼ ∇π(k) ∼ k 2 on large scales, k → 0.

A Second order counter terms at 1-loop order
This section is devoted to a detailed calculation of second order counterterms in the leading loop order. This is sufficient for a 1-loop calculation of the bispectrum. However, the resulting counterterms renormalize all reducible diagrams where the vertex corrections of Fig. 1a and 1b appear as a sub-diagram, which is the goal of systematic renormalization. For this the counterterms are unavoidably nonlocal in time since the ingoing lines into this sub-diagram can be of different orders and hence have different time-dependence. This is the practical manifestation of the "long memory effect" discussed in section 4.1.

A.1 Shift terms
As discussed in section 4.1 promoting the time integration to be along the flow identifies several higher order counterterms as the displacement a lower order one. These shift terms can be isolated by looking at the infra-red limit when one of the ingoing momenta becomes much smaller than the other. In other words, the two vertices v i ∂ i θ and v i ∂ i δ with soft velocities lead to IR-divergences as 1/k in the momentum of Ψ 2 ∼ θ = ∇.v. In the case of 1-loop second-order counterterms, the first order shift of a first order counterterm (4.13) automatically renormalize these "IR-divergent" cutoff dependence. It is shown in section 4.5.2 that these non-original counterterms are , the above correction to the equation of motion can be written in Fourier space as Now let us examine the IR-limit of the corrections to the vertex function. Using the following identity between step functions one readily finds that IR-limit of diagram.1 contribution to the vertex correction as where the symmetry factor 2 comes from (1 ↔ 2) or more explicitly (k 1 ↔ k 2 , η 1 ↔ η 2 , a 1 ↔ a 2 ). And Again the third line is the soft limit amplitude of the inverted version of Diagram.2. Adding up these contributions together one recover the IR-divergent vertex correction already anticipated by promoting time integrals in propagator correction integral, C U V a , to an integral along x fl , which leads to (A.8).

A.2 The new counterterms
As discussed in 4.5.3, after the subtraction of shift terms from 1PI diagrams, what remains must be renormalized with the new higher order counterterms. For every second order term (counter-term) in the equation of motion there is an associated cubic vertex in Feynman diagrams so we use 1-loop second order counterterms and 1-loop cubic vertex correction interchangeably. Using Eq. (6.2), all second order counterterms are of the following form The general structure of F i kernels can be found using Eq. (4.20) as As already emphasized, although the basis 4.20 for tensors is non-redundant, taking the derivatives leads to degeneracies. It can be shown that the above structures are linearly dependent:F 6 −F 5 = F 3 − F 4 . Hence, we eliminateF 6 −F 5 and define Having found all possible momentum dependencies, the next step would be listing the time kernels. There is a unique basis for all possible time kernels of 2-nd order counterterms in the equation of motion. We denote them byκ i given bŷ Every time kernel has three indices associated with indices of the one outgoing and two ingoing legs respectively. As every index can be either 1 or 2 there are 2 3 = 8 distinct counter terms which are listed below Figure 2: Hard ingoing lines combining into a soft line.
in which Θ 12 = θ(η − η 1 )θ(η 1 − η 2 ) and Θ 21 = θ(η − η 2 )θ(η 2 − η 1 ). Moreover i K i 112 F i and i K i 212 F i could be simply read off from i K i 121 F i and i K i 221 F i , respectively, by the following replacementsκ 2 ↔κ 3 ,κ 4 ↔κ 5 , Note that the counterterms for the Euler equation, which are the ones with the first index 2, are all proportional to which in real space translate into ∇ 2 δ 2 and ∇ 2 (∂ i ∂ j φ) 2 , respectively. Hence, as expected they are Laplacians of locally observable scalar quantities and the corresponding force in the Euler equation is curl-free.

B Softness properties of perturbation theory
In this appendix we give a perturbative and iterative proof of the single softness of momentum π(k) and double softness of density contrast δ(k) which result from evolution of short wavelength initial modes -also known as decoupling. The proof is broken into several steps.

B.2 Soft and hard ingoing modes
In the presence of both soft and hard initial modes, double softness can be established by dividing the quadratic vertices of the perturbation theory into different categories. First, every vertex that results in a soft outgoing field either combines two short modes or two long modes, e.g.
[∂φ∂φ] l = ∂φ s ∂φ s + ∂φ l ∂φ l . (B.9) In the presence of initial soft modes we need to make a further distinction of whether a long mode is purely made of those initial long modes (indicated as ∂φ l;l ), or if it contains initial hard modes in addition to a (possibly empty) set of initial long modes (indicated as ∂φ l;sl ). Thus, the first vertex in (B.2) is decomposed schematically as ∂∂[∂φ s ∂φ s + ∂φ l;l ∂φ l;sl + ∂φ l;sl ∂φ l;sl ]. (B.10) Note that ∂φ l;l ∂φ l;l is excluded since we are only interested in the case where there are some hard initial modes. Figure 3 shows one example for any of the above three terms.
1. The first term in the above equation is similar what we considered in point 1 of the previous section, where the two external derivatives ensured double softness. However, there is an important difference: the short modes can now contain initial long modes (see figure 3a). One may therefore worry that ∂φ s may contain negative powers of k l due to the displacement of hard modes by the soft modes. That is, terms in perturbative evolution of the form v l;l ·∇δ s where v l;l = O(k −1 l ). These IR singularities indeed exist in any short wavelength quantity that evolves from initial hard and soft modes. They correspond to the motion of short modes with the bulk flow of the long modes and their presence is ensured by the equivalence principle. However, as shown more explicitly in appendix C whenever a collection of hard modes combine into a soft mode, adding any number of initial soft modes to such diagram does not lead to any negative power of k l . The singular displacement of hard lines with momenta {q i } by a soft line of momentum k 1 that are proportional to k 1 · q i /k 2 1 add up to k 1 · k total /k 2 1 − 1 = O(k 0 l ) at any moment. Hence, the first term in (B.10) leads to a double soft outgoing δ l (and as before a single soft π l ).
2. The second and third terms can also be shown to lead to at least double soft results using induction. First suppose there is no initial long modes in ∂φ l;sl . Then by the argument of the last section ∂φ l;s = O(k l ). And ∂φ l;l = O(k −1 l ) from (B.5). So the resulting outgoing δ l;sl = O(k 2 l ) because of the overall derivatives. One can then continue the induction provided that the same argument applies to the ∂∂(vπ) vertex of (B.2).
3. Decomposing vπ into short and long wavelength contributions we get [vπ] l = v s π s + v l;sl π l;l + v l;sl π l;sl + v l;l π l;sl . (B.11) Assuming that up to some order in perturbations and using the fact that v l;l , π l;l = O(k −1 l ) and δ l;l = O(k 0 l ), we see that the last two terms on the r.h.s. of (B.11) preserve these properties (when considering the vertex ∂∂(vπ)). The second term clearly violates them and the first term does so too: expanding the last term is seen to be O(k −1 l ), leading to O(k l ) contribution to δ l from the vertex ∂∂(v s π s ). However, the sum of the two dangerous terms preserves the condition (B.12). Substituting π l;l = (1 + δ) l;l v l;l (B.14) in the second term on the r.h.s. of (B.11) and using the fact that (1 + δ) l;sl = δ l;sl (since the subscript s indicates a nonzero set of initial short modes) we obtain for the sum of potentially dangerous terms Note finally that all of the above arguments generalize to the case where counterterms of the form ∂ j τ ij , with τ ij a locally observable tensor which by definition can never have an IR singularity (corresponding to a negative scaling with k l ), are added to the Euler equation (5.9). This seems to be the most general way of modifying the equations of motions that is still compatible with locality and momentum conservation. Intuitively, the short scale dynamics can generate an effective force on a volume element of the fluid only due to its momentum flux into that volume.

B.3 Velocity field
The final point to be discussed is the softness property of the mass weighted velocity field defined as v = π/(1 + δ). Using mass conservation we get where we neglected vorticity and set v = −H∇φ v . From these equations it follows that doublesoftness of a soft outgoing δ evolved purely from initial hard modes implies double softness of θ.
To show this one needs to argue that the second term in (B.16) which naively is only single soft is actually double-soft when the perturbative solution for δ and θ are substituted. This follows from the fact that the second vertex on the r.h.s. of (B.17) is automatically double soft and if δ is to be double soft the first term on the r.h.s. must also be double soft by itself. 11 If we use δ and v as primary variables -as we did in our systematic renormalization -and add counter-terms proportionally to SPT cutoff dependence, the resulting renormalized velocity inherits both curl-freedom and single-softness from SPT kernels. This renormalized velocity of course differs from v π . On the other hand, if we use v π and add counterterms as ∂ j τ ij to the momentum equation (5.9), they modify (B.17) as ∂ i [(1 + δ) −1 ∂ j τ ij ]. Expanding this leads to vertices which are not explicitly second derivative. For instance, taking τ ij = δ ij (∂ m ∂ n φ) 2 we get an associated cubic vertex Taking all ingoing modes to be hard, this vertex gives rise to an outgoing δ which is only single soft. In this formulation, the full result becomes double soft (as it must according to the arguments of sections B.1 and B.2) because the above single soft term cancels with the one arising from evolving a second order solution coming from the leading part of the counterterm ∇ 2 (∂ m ∂ n φ) 2 to third order via the vertex ∂ i (δv i ). Hence, the double softness of δ does not imply that of θ π anymore because the second term on the r.h.s. of (B.16) is generically only single soft after the addition of counterterms. Note that one can redefine velocity to absorb these single soft terms, but the new velocity is no longer given by π/(1 + δ).

C Cancellation of IR singularities
Consider the following time evolution diagram with ingoing momenta {k 1 , q 1 , q 2 } and outgoing momentum k, and such that k 1 ∼ k ≪ q 1,2 .
The potential IR singularity arises in this diagram from the vertex v · ∇Ψ a , where v realizes the soft k 1 mode and Ψ a a hard one, leading to k 1 · q i /k 2 1 . Focusing only on this vertex we have for the outgoing field (C.2) 11 One should also consider the case when the hard modes combine into two soft modes which later combine into a final δ. In this case too the second term on the r.h.s. of (B.17) is double soft because at least v l;s = O(k 0 l ). Now using Ψ d (q 2 , η 1 ) = g df (η 1 , η 2 ) Ψ f (q 2 , η 2 ), (C. 3) and the fact that γ bcd (k, q 1 + k 1 , q 2 ) = γ bcd (k, q 1 , q 2 ) + O(k 1 ) in the limit k 1 ≪ q i , we obtain f (q 2 , η 2 ). (C.4) The expression in the square brackets is k 1 · (k − k 1 )/k 2 1 which is O(k 0 ), thus the potential IR singularities coming from the v · ∇Ψ a vertex cancel in this diagram when we sum over the two possible attachments of the soft line. The full diagram will therefore have the same softness property as the one without attaching the soft line.
The above cancellation arises from the fact that the displacement produced by the long k 1 mode is universal, and hence when the individual hard momenta add up to a soft momentum k, the displacement terms also add up as effectively moving a long wavelength mode. Note that the cancellation happens independently at any time-slice. Next we show this to hold more generally whenever one or several ingoing soft modes are added to a diagram that combines a set of hard modes into a soft mode.
Let a i (q i , η 1 )} whose contribution to the outgoing field is A aa 1 ··· (η; η 1 , {q i }) i Ψ (n i ) a i (q i , η 1 ), (C. 6) where A a··· is a function that contains all momentum dependences arising from the vertices after η 1 and all propagators. Now add the soft line and sum over all possible attachments of v · ∇ vertex at η 1 i A aa 1 ··· (η, η 1 , {q 1 , · · · , q i + k 1 , · · · }) k 1 · q i k 2 1 Ψ (m) 2 (k 1 , η 1 ) i Ψ (n i ) a i (q i , η 1 ), (C.7) and use the fact that A a··· (η, η 1 , {q 1 , · · · , q i + k 1 , · · · } is regular when k 1 ≪ q i to get The sum gives k 1 · (k − k 1 )/k 2 1 = O(k 0 l ). Hence, the drifts caused by a long wavelength mode add up at any moment to the overall drift of the total (stripped) diagram whose momentum is k − k 1 . This proves that if we start from a composite operator such as ∂φ s ∂φ s , which combines hard modes into a soft outgoing line, and perturbatively add soft modes they will not introduce any negative power of k l .

D Equations in Fourier space
In this appendix we provide explicit formulae for the equations of motion of momentum and velocity in Fourier space.
The equations of motion using momentum are ∂ τ δ + µ = 0 , (D.1) Several definitions are in order. ξ δδ comes from the term δ∂ i φ where α was given in (2.9) and the arrow refers to symmetrization 12 The other interactions come from ∂ i π i π k /ρ l and are as follows.

E Glossary
In this paper, we have extensively used the same nomenclature as in quantum and statistical field theory to stress the similarity with the problem at hand and to make the discussion more intuitive. However, when needed we have adapted and refined some of the definitions. For the convenience of the reader we have collected below those technical terms used in this paper whose meaning differs in some important or subtle way from the common usage. 13 We used the magic vector identity 1-Particle Reducible (1PR) A loop diagram for a correlation function that can be divided into two disconnected parts by cutting an internal line that is not immediately connected to a contraction or an external line, e.g. (3.4)

1-Particle Irreducible
Any loop diagram that is not 1PR, e.g. (3.1) Time-evolution diagram It shows the evolution of a set of initial fields into one higher order final field at a later time (3.6)

Hard line
In a diagram, a line representing the time evolution of perturbation δ(k) with a large momentum k → ∞, e.g. thick lines in (3.6)

Soft line
In a diagram, a line representing the time evolution of perturbation δ(k) with a small momentum, e.g. thin lines in (3.6)

Amputated diagram
A time evolution diagram in which the soft lines are immediately connected to the hard lines, e.g. (3.11)

Non-stochastic
The counterterm of a time-evolution diagram in which all of hard lines are paired and contracted with each other, e.g. (3.8)

Stochastic
The counterterm of a time-evolution diagram in which one or more hard ingoing lines remain unpaired, e.g. (3.9) Double softness A field is double soft iff it scales as k 2 as k → 0, e.g. δ and ∇π

Systematic
Renormalization is systematic if new counterterm are needed only for 1PI diagrams, while divergences in 1PR diagrams are already canceled by lower order counterterms