{'fraction_non_alphanumeric': 0.04065600882028666, 'fraction_numerical': 0.013850606394707829, 'mean_word_length': 4.382382939267501, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Renormalized Halo Bias

This paper provides a systematic study of renormalization in models of halo biasing. Building on work of McDonald, we show that Eulerian biasing is only consistent with renormalization if non-local terms and higher-derivative contributions are included in the biasing model. We explicitly determine the complete list of required bias parameters for Gaussian initial conditions, up to quartic order in the dark matter density contrast and at leading order in derivatives. At quadratic order, this means including the gravitational tidal tensor, while at cubic order the velocity potential appears as an independent degree of freedom. Our study naturally leads to an effective theory of biasing in which the halo density is written as a double expansion in fluctuations and spatial derivatives. We show that the bias expansion can be organized in terms of Galileon operators which aren't renormalized at leading order in derivatives. Finally, we discuss how the renormalized bias parameters impact the statistics of halos.


Introduction
Over the next decade, large-scale structure (LSS) surveys will play an increasingly important role in the measurement of cosmological parameters and as a probe of initial conditions. In order to relate late-time observables to the physics of the early universe, several sources of secondary non-linearities need to be understood (see fig. 1). Reducing the theory error is essential if the full potential of future surveys is to be realized. 1 Non-linearities in the gravitational evolution can be characterized by numerical N-body simulations [3] and, on sufficiently large scales, by perturbation theory [4,5]. Less well understood are non-linearities in the biasing between the clustering of galaxies and the underlying dark matter density. The biasing problem is already visible in dark matter-only simulations, where it is reflected in the biasing of dark matter halos. On large scales, linear biasing has been shown to be a good approximation: where δ h and δ are the density contrasts of the halos and the dark matter, respectively, and the bias parameter b 1 is an unknown coefficient (to be fit to data). However, linear biasing is known to fail on small scales where non-linearities becomes important. One common procedure for describing halos beyond the linear biasing model is local Eulerian biasing [6] which assumes that the halo density contrast is a local function of the dark matter density, δ h (x, τ ) = F[δ(x, τ )]. Formally, we might write this relation as a Taylor expansion Local biasing is motivated both as a natural generalization of linear biasing and as a consequence of a number of semi-analytic models of halo formation. It is also often employed in data analysis [7][8][9][10][11][12]. However, the meaning of (1.2) is far from clear, as we need to define δ n (x, τ ) for n > 1. In particular, non-linear quantities like δ n receive contributions from all scales and are therefore not necessarily small, even on large scales. A common procedure is to smooth δ with a window function that removes power below some length scale Λ −1 . However, this only masks the problem. The scale Λ −1 is arbitrary and should not appear in physical quantities. 2 In [13], Mc-Donald showed how to reorganize the bias expansion in a way that makes this property manifest (see also [14,15]). In this paper, we revisit his idea of renormalized halo bias.
The key feature of renormalization is that, although δ n may be large (or depend on an unphysical smoothing scale Λ), all large (or Λ-dependent) contributions can be systematically removed by adding local counterterms We will find that consistent renormalization requires additional fields O beyond those appearing in the local Eulerian biasing model (1.2). For instance, at quadratic order, we have to add the tidal tensor (∇ i ∇ j Φ g ) 2 [14,16,17], where Φ g is the gravitational potential, while, at cubic order, the velocity potential Φ v has to be introduced as an independent degree of freedom. In addition, higher-derivative terms must be included to remove the sub-leading Λ-dependence. Finally, no stochastic parameters are needed to remove divergences in correlation functions at separated points, although one may choose to include them to remove delta-function localized divergences.
In general, we find that all the terms that are allowed by symmetry are required in order to define the renormalized fields 3 Renormalized halo biasing is naturally an effective theory, in which the halo density contrast is written as a double expansion in terms of powers of the fluctuations and their derivatives. In order to be well-defined, this expansion has to be organized in terms of the renormalized quantities: O are the renormalized bias parameters which by construction do not depend on the smoothing of the density field. It is to be expected that the biasing model (1.4) will contain all terms O consistent with the symmetries of the dark matter equations of motion, as was emphasized by McDonald and Roy [14] (see also [20]). 4 After renormalization, terms with higher number of fields and higher number of spatial derivatives are suppressed. To describe the halo statistics to a finite precision then only requires a finite number of terms in (1.4). How many terms need to be retained depends on the power spectrum of the initial conditions in the same way as in the effective field theory of large-scale structure (EFT-of-LSS) [21][22][23]. In principle, the different terms in the bias expansion can be distinguished by measuring correlation functions of the dark matter halos and fitting for the bias parameters b (R) O . For the gravitational tidal tensor 2 Semi-analytic models only add to the confusion, as they often identify Λ −1 with the Lagrangian size of the halo, which is a physical scale. 3 In this paper, we will only consider the case of Gaussian initial conditions. Additional terms can arise for non-Gaussian initial conditions [18]. We will discuss this elsewhere [19]. 4 The form of (1.4) was also argued to arise from the non-linear time evolution of local Lagrangian biasing [16].
this has been demonstrated in [17,24]. Here, we lay out the basic steps towards a systematic treatment of the effective theory of halo biasing.
The outline of the paper is as follows: In Section 2, we give a pedagogical introduction to renormalization in structure formation. We present the renormalization of δ 2 as an explicit example. In Section 3, we generalize this to all terms in the local Eulerian biasing model. We classify the terms that have to be added to the biasing model to make it consistent with renormalization. In Section 4, we show how these terms affect halo statistics at one-loop order. Our conclusions appear in Section 5. Four appendices contain technical details: In Appendix A, we discuss how the renormalized quantities depend on the renormalization scale. In Appendix B, we present the renormalization of δ n , for n > 2. In Appendix C, we prove a non-renormalization theorem for Galileon operators. In Appendix D, we give the details of the one-loop halo bispectrum.
We will use the following notation and conventions: Conformal time is denoted by τ and comoving coordinates by x. Overdots denote derivatives with respect to conformal time and H ≡ȧ/a is the comoving Hubble parameter. For momentum integrals, we use the shorthand 3 . We define the normalized inner product of two vectors as µ p,q ≡ p · q/pq and sometimes use σ 2 p,q ≡ µ 2 p,q − 1. A prime on a correlation function, . . . , indicates that the overall delta function is dropped. We will use P (q) ≡ δ for the linear power spectrum (which is often called P 11 ) and P mm (q) ≡ δ q δ q for the non-linear power spectrum. Sometimes, we write P i ≡ P (q i ). Renormalized operators are denoted by square brackets, [O]. We use two definitions for the bias parameters of δ n , namely b n and b δ n , which are related by b δ n ≡ b n /n!. When we present numerical results, the linear power spectrum is computed with CAMB [25], using a ΛCDM cosmology with the standard cosmological parameters, Ω Λ = 0.73, Ω b h 2 = 0.02, Ω m h 2 = 0.12, h = 0.7, and power law initial conditions, with n s = 0.96 and A s = 2.2 × 10 −9 .

Renormalization in Structure Formation
Before explaining the renormalization of local Eulerian biasing, we first review, in §2.1, dark matter perturbation theory in light of the recently developed EFT-of-LSS [21,26] (see also [22,23,[27][28][29][30][31]). We then discuss, in §2.2, the renormalization of the simplest non-linear term in the local Eulerian biasing model, namely δ 2 . The complete renormalization and further technical details will be presented in Section 3 and in the appendices.

Fluid Equations
On large scales, dark matter acts as a pressureless fluid. In the Newtonian approximation, the dark matter density contrast δ ≡ δρ/ρ and the dark matter velocity v satisfy the following evolution equations where Φ g is the gravitational potential. In a matter-dominated universe, Φ g is sourced by δ where Ω m = 1 in an Einstein-de Sitter (EdS) universe. If the velocity field is irrotational, it is described completely by its divergence where Φ v is the velocity potential. We will use Φ g and Φ v as our fundamental fields. As a consequence of the equivalence principle, the evolution equations are invariant under the boost symmetry x where T (τ ) is an arbitrary function of the conformal time τ and u is the boost velocity. The boost transformation (2.5) shifts the gradient of the gravitational potential and the velocity potential, Finally, in EdS, the fluid equations are invariant under a Lifshitz symmetry for a generic weight z. The potentials transform as For power law initial conditions, δ q (τ in ) δ q (τ in ) ∝ q n , the rescaled solutions have the same initial conditions iff 4z = n + 3. Further discussion of the symmetries of the Newtonian fluid equations can be found in [32,33]. Conserved tracers of the dark matter density satisfy the same equations and the same symmetries as the dark matter perturbations. Of course, in reality, the number of halos is not conserved, so the evolution equations will contain extra source terms. However, we will only require that the halo density contrast δ h is a scalar under the transformation (2.5). Any equation governing the time evolution of halos will only impact the time evolution of the bias parameters which is beyond the scope of this work.

Standard Perturbation Theory
We briefly review standard perturbation theory (SPT) in an Einstein-de Sitter universe [4]. It is convenient to write the equations of motion (2.1) and (2.2) in Fourier spacė where the left-hand sides capture the linear evolution and the right-hand sides contain non-linear convolutions (2.14) Eqs. (2.11) and (2.12) can be solved once the initial conditions for δ and θ have been specified. At sufficiently early times, the dark matter density contrast and the velocity divergence are small. Consequently, the equations of motion can be solved order by order in the initial conditions. In an Einstein-de Sitter universe, the solution can formally be written as 5 where only the growing mode has been kept. 6 In Fourier space, the fields δ (n) and θ (n) can be written as where the kernel functions F n and G n can be computed iteratively [4]. In particular, the firstorder kernels are just F 1 = G 1 = 1, while the second-order kernel are given by 7 with µ 12 ≡ q 1 · q 2 /q 1 q 2 . The initial conditions are encoded in the linear dark matter density contrast δ (1) q (τ in ). We will assume that the initial conditions are Gaussian, 8 so that the statistics of the initial fluctuations δ (1) is completely determined by its power spectrum (2.20) The computation of correlation functions can be organized using Feynman diagrams. Each field δ (n) is represented by a vertex with n external legs δ (1) : Although strictly speaking this form of the result is only valid for the EdS universe, it can approximately be extended to arbitrary values of Ωm and ΩΛ, by replacing the EdS growth factor a by the growth factor of the corresponding cosmology D(a) [4]. 6 Formally, this is equivalent to imposing the initial conditions at τin → 0. 7 We will use the symmetrized kernels, obtained by summing over all permutations of the momenta. 8 We will treat the non-Gaussian case in [19].
To compute an N -point correlation function in an EdS universe with Gaussian initial conditions, we use the following Feynman rules: 1. Draw every connected graph with N vertices.
3. To each propagator assign the time-evolved linear power spectrum, 9 where q is the magnitude of the momentum flowing in the propagator.

Integrate over each loop with measure of integration
taking into account the symmetry factor of the loop.

Effective Theory of Large-Scale Structure
As an example, we compute the dark matter power spectrum δ q δ q ≡ P mm (q) at one-loop order, i.e. at fourth order in the initial conditions δ (1) : where P nm ≡ δ (n) δ (m) and P 11 ≡ P . The first one-loop contribution, P 13 , has the following explicit form This is an integral over all comoving wavenumbers p, including those for which perturbations have already grown non-linear and are therefore outside the regime of validity of SPT. Consequently, one cannot take this first-order correction seriously, regardless of whether this integral diverges or not. To make sense of perturbation theory at the loop level, an effective field theory (EFT) approach has recently been developed [21,26] (see also [22,23,[27][28][29][30][31]). This approach describes the dark matter on scales larger than some cut-off scale Λ −1 k −1 NL , while systematically keeping track of the effects of short-distance fluctuations on long-wavelength observables through modifications to the Euler equation (2.2). These new "fluid" equations modify the SPT solution δ SPT to δ = δ SPT +δ , (2.28) and similarly for the velocity divergence θ. Here,δ is the new solution generated by the new parameters in the fluid equations. This new solution can also be written as an expansion in the initial conditions [23]δ The new kernel functionsF n depend explicitly on the cut-off Λ in such a way as to cancel the cutoff dependence arising from loop diagrams. We will use the following diagrammatic representatioñ . . =F n (q 1 , · · · , q n ; Λ) (2π) 3 δ D (q 1 + · · · + q n − q) . (2.31) The renormalized one-loop power spectrum can be written diagrammatically as where a counterterm proportional toδ (1) was added to cancel the cut-off dependence coming from P 13 . The second one-loop contribution, P 22 , can be written as the sum of a finite (Λ-independent) part and a divergent (Λ-dependent) part. However, the divergent part is proportional to integer powers of q 2 which correspond to delta-function localized contributions in position space. Such contact terms vanish when the two-point function is evaluated at separated points. We will ignore these contact terms and therefore our formulas will be correct 10 up to terms analytic in the momenta. 11

Renormalized Halo Bias
We are now in a position to illustrate the basic idea of renormalized halo biasing. We will first do this with a specific example, the renormalization of δ 2 , leaving a complete treatment to the next section. The results presented here have overlap with [13,14,34]. Our general approach to renormalization agrees with these works, although it differs in detail. One advantage of our approach is that it is systematic and can be carried out unambiguously order-by-order in perturbation theory. Moreover, we will clarify how the precise results depend on the choice of renormalization condition (and the renormalization scheme). 10 For correlation functions of conserved quantities like δ, some contact terms can be forbidden. We will not assume any conservation law for halos and therefore their correlation functions may include all possible contact terms. 11 An alternative approach is to introduce an additional stochastic variable that has delta-function localized correlation functions to remove these contact terms (see e.g. [13,26]).

Example: Renormalization of δ 2
Non-linear biasing contains terms with products of fields evaluated at the same point, e.g. δ n (x), for n > 1, in local Eulerian biasing (1.2). In field theory these terms are called composite operators and we will often use that terminology. Composite operators can lead to additional divergences that cannot be removed by the renormalization of the dark matter density contrast-i.e. by counterterms made out of (2.31). The simplest composite operator which appears in (1.2) is δ 2 .
In this section, we will investigate the renormalization of this object. In fact, as we show in Appendix B, this will be the essential building block for the one-loop renormalization of the other composite operators appearing in the local Eulerian biasing model, namely δ n , for n > 2.
The composite operator δ 2 can be made finite by defining a new renormalized operator [δ 2 ] as where the operators O are counterterms introduced to absorb the UV divergences that arise in correlation functions of δ 2 . We need to distinguish the divergences which arise from nonlinearities in the "external legs"-i.e in the dark matter contrast δ q i -and which are removed by the countertermsδ in the EFT-of-LSS, from the divergences which arise from contractions within the composite operator δ 2 and are renormalized by the counterterms in (2.33). To single out the internal divergences coming from the operator δ 2 , we replace the external legs δ q i by their linear approximations δ (1) q i and impose the renormalization conditions where the subscript "tree" indicates the tree-level result. Eq. (2.34) therefore implies that the counterterms are chosen to precisely cancel the loop divergences where the subscript "loop" refers to diagrams which contain one contraction between two different fields within the composite operator δ 2 . We see that, at one loop, the counterterms only need to be evaluated at tree level. We will discuss the renormalization conditions in more detail in §3.1.2. Here, we will simply illustrate them in an example.
We construct the renormalized operator [δ 2 ] order-by-order in perturbation theory. To facilitate the calculation, we introduce the following diagrammatic representation: where the "blob" indicates a convolution over the momenta.
• m=0.-The expectation value of δ 2 is where Λ is a momentum cut-off introduced to regulate the loop. Of course, this divergence (or cutoff-dependence) is removed by adding a constant counterterm Subtracting this tadpole contribution ensures that δ h = 0 at the loop level. We note that the constant counterterm does not affect correlation functions with m > 0, as it only contributes to disconnected graphs.
We see that the dependence on the cut-off can be removed by a counterterm proportional to δ, where we have introduced a "crossed circle" to denote this new counterterm.
While the first divergence can be absorb by a counterterm proportional to δ 2 , the second term cannot be removed by a counterterm which is local in δ. Instead, the counterterm which removes this divergence is whereΦ g is the rescaled gravitational potential. In order to keep the notation clean, we will drop the tilde and from now on denote the rescaled gravitational potential by Φ g . The renormalized operator [δ 2 ] then is where the ellipses refer to higher-derivative operators which are required to remove the subleading divergences of (2.42). Importantly, we see that the renormalization of the simplest operator in the local Eulerian biasing model requires the presence of an operator which is non-local in δ.

Comments
Before we move on to the more general treatment of renormalization, let us make a few comments: • The fact that renormalization requires non-local biasing is fundamentally different from the results of [16], which found that non-local bias is generated by time evolution. In order to derive the latter result, one needs to assume an equation for the time evolution of δ h . Instead, we are simply demanding that the predictions for correlation functions of halos are independent of how the theory is regulated. Since the choice of regulator is arbitrary (and hence unphysical), our results follow only from the assumption that halos are physical. Furthermore, since δ is itself a renormalized quantity, our counterterms are proportional to b n≥2 (see §3.2) and not b 1 as is found in [16].
• The operator δ 2 is a scalar under homogeneous boosts. Therefore, δ 2 can only be renormalized by operators which are also scalars under these transformations. From eq. (2.6), we see that the gravitational potential can only appear with (at least) two spatial derivatives acting on it, which is the case in (2.44). The gravitational potential itself only becomes an allowed operator if the symmetry is broken by the initial conditions [18] or if the equivalence principle is violated in the dynamics [20,35].
• At the order we have worked so far, the loop corrections have generated all operators which are scalars under boosts. In fact, this is quite general: loop corrections generate every operator consistent with the symmetries (here, the invariance under homogeneous boosts).
Consequently, the bias relation should not only include the local operators δ n , but every other operator allowed by these symmetries.
• It is important to note that each diagram in (2.41), when evaluated separately, has divergences which cannot be removed by a scalar operator. In particular, we find divergences such as q −2 i or q −2 i (q 1 · q 2 ) which can only be removed by the operators Φ g δ and ∇ i Φ g ∇ i δ. It is only when these three diagrams are summed that these undesirable divergences cancel.
• It is convenient to rewrite the counterterms in (2.44) as (2.46) As we will explain in §3.1.3, and prove in Appendix C, the operator (2.46) is part of a larger class of operators-the so-called Galileon operators [16,36]-which are not renormalized at leading order in derivatives. In other words, the divergences associated with these operators are only renormalized by derivative operators such as ∇ 2 δ.
• We have shown explicit results at leading order in derivatives. Keeping the subleading terms, we find that higher-derivative operators are also required by renormalization Note that the operator ∇ 2 δ does not appear at one-loop order. This is due to the fact that eq. (2.39) is an exact result, and hence there is no divergence that scales as q 2 . However, as we will explain in §3.1.3, the operator ∇ 2 δ is generated by the one-loop renormalization of the Galileon operator G 2 (Φ g ).
• The bias relation is an expansion both in small fluctuations and in derivatives. As we have explained, the expansion in small fluctuations makes sense only for renormalized operators. Similarly, the expansion in derivatives is only well-defined once the derivative terms have been appropriately renormalized. To be concrete, let us consider, at one-loop order, the operator (∇δ) 2 and its correlation with δ (1) : (2.49) We see that the one-loop contribution from the bare operator (∇δ) 2 does not vanish in the limit q → 0. It is precisely these non-vanishing contributions that are removed when we define the renormalized operator [(∇δ) 2 ]. Therefore, the derivative expansion is well-defined after the operators have been correctly renormalized.
• The velocity potential starts appearing at cubic order through the following operator is the rescaled velocity potential and the subscript "3" on Γ 3 was added to remind us that this is a cubic operator. For lack of better term, we will call this operator the "velocity tidal tensor". This operator is required to renormalize δ 2 at m = 3. In §4.4.2, we will see that the velocity potential leaves a distinct imprint in the halo trispectrum.

One-Loop Renormalization of Halo Biasing
We now generalize the δ 2 -example of the previous section to a complete treatment of renormalized halo biasing. In particular, we wish to extend our analysis to more general composite operators O (i.e. other products of two or more fields at coincident points). The goal is to define a biasing model that is manifestly independent of the way the theory is regulated and write the bias expansion in terms of renormalized operators [O] [37,38]. By construction, correlation functions of [O] will be finite modulo divergences that are renormalized by the EFT-of-LSS and divergences that correspond to contact terms. In §3.1, we describe the technical details associated with the construction of the basis of renormalized operators. In §3.2, we present explicit results for the renormalized bias parameters at one-loop order.

Renormalized Operators
The renormalized composite operators [O] are defined in terms of unrenormalized bare operators O by [37,38] [ where coefficients Z O, O depend on the cutoff Λ, but the operators [O] are independent of how the theory is regulated. The counterterms on the right-hand side of (3.1) typically contain every term consistent with the symmetries. We will first list the available terms, up to third order in perturbation theory and at lowest order in derivatives. As we will see, this analysis is complicated by the fact that some operators are related by the equations of motion and are therefore not independent. For related discussion, see [14,16,20].

Symmetries and Counterterms
In §2.1.1, we have seen that the gradients of the gravitational potential Φ g and the velocity potential Φ v shift by a vector under homogeneous boosts. If the initial conditions are also invariant under such boosts, the allowed counterterms should only involve operators which are scalars under these transformations. Consequently, the basic building block for constructing the renormalized theory is Although we will mostly work at leading order in derivatives, the results of this section can be generalized straightforwardly to include higher-derivative terms. Rotational invariance implies that the indices in (3.2) need to be contracted. At first order, the bias relation can therefore depend on while at second order, it may contain the following terms At a given order in perturbation theory, the operators in (3.3) and (3.4) are not all independent, but are related by the dark matter equations of motion: • 1st order.-At linear order, the (rescaled) gravitational and velocity potentials are equal Hence, there is just one independent operator in (3.3), which we choose to be δ = ∇ 2 Φ g .
v , we do not distinguish between the gravitational potential and the velocity potential in the list of quadratic operators in (3.4), which therefore only contains two independent operators: Moreover, it is easy to show from (2.18) and (2.19) that, at second order, the difference between the density contrast and the velocity divergence is the Galileon operator (2.46), Hence, although a priori eight operators are consistent with the symmetries, only three are independent after using the second-order equations of motion. We chose these independent operators to be δ , δ 2 , G 2 (Φ g ) . (3.7) • 3rd order.-At cubic order, the velocity potential appears for the first time as an independent degree of freedom. Indeed, a set of independent operator at this order is [16] δ , δ 2 , δ 3 , and Γ 3 contains the velocity tidal tensor G 2 (Φ v ).
Of course, this analysis can, in principle, be carried out to any order in perturbation theory. However, the set of independent operators listed in (3.8) will be sufficient for most of this paper (but see Appendix D). We will formulate renormalization in terms of this basis of operators.

Renormalization Conditions
We are interested in correlation functions of the form O q δ q 1 · · · δ q m . These correlation functions contain two types of divergences. Those associated with non-linearities in the external δ's (which are renormalized by the EFT-of-LSS) and those within the operator O. Since we will be interested in the latter, we may replace the external legs δ q i by their linear approximations δ (1) q i and restrict to (partially) one-particle irreducible ((p)1PI) diagrams (see fig. 2). In order to fix the finite contributions in the counterterms in (3.1), we impose renormalization conditions. We will require that the counterterms exactly cancel the divergences on large scales, i.e. in the limit where the external momenta vanish. In other words, we define the renormalized operators [O] through 12 where the subscript "(p)1PI" denotes (partially) 1PI diagrams. In Appendix A, we discuss the consequences of imposing the renormalization conditions at finite momentum, q i = µ. In particular, we show that the basis of renormalized operators that is orthogonal at a scale µ will mix and in general won't be an orthogonal basis at another scale µ . Related to this is the fact that the renormalized bias parameters (see §3.2) depend on the renormalization scale µ.
One may be concerned that the renormalization conditions (3.10) are not enough to ensure that every UV divergence which can appear in O q δ q 1 · · · δ q m has been removed. Indeed, the physical object we wish to compute is not the 1PI part of the correlation of O with the linear dark matter contrast, but the connected part of the correlation function O q δ q 1 · · · δ q m . However, since at any given order in perturbation theory the non-linear dark matter contrast δ can be written as a product of linear dark matter contrasts δ (1) , these (partially) 1PI diagrams are the building blocks of correlation functions with the non-linear dark matter density contrast. Any remaining loops are either renormalized by counterterms in the EFT-of-LSS or are finite up to contact terms.

Non-Renormalization of Galileon Operators
In eqs. (2.46) and (3.9), we defined the Galileon operators G 2 and G 3 . These definitions can be extended to n-th order Galileon operators, whose precise expressions can be found in Appendix C. Importantly, every Galileon operator G n (Φ) can be written as the second derivative of another operator T As a result, the correlation of G n with (δ (1) ) m scales as q 2 , provided that [Φ T isn't singular when any of the external momenta (or partial sums of the external momenta) go to zero. As we explain in Appendix C, the absence of such poles is guaranteed by symmetry considerations. Consequently, loops arising from Galileon operators are only renormalized by higher-derivative operators At leading order in derivatives, but to all orders in loops, the Galileon operators G n are not renormalized.

Renormalized Bias Parameters
In terms of the bare operators, the halo density contrast is q m contain divergences proportional to positive powers of the external momenta, additional renormalization conditions need to be imposed on derivatives of the (amputated) correlation functions with respect to these momenta.
If the expansion contains all operators consistent with the symmetries, then it can be re-written in the basis of renormalized operators O are the renormalized bias parameters. We can gain intuition for the form of the renomalized bias parameters from the results of §2.2. Consider the bare expansion at quadratic order, To write this in terms of [δ 2 ], we add and subtract the counterterms in eq. (2.45). This shifts b Notice that, at one-loop order, the operators δ 3 and δ 4 do not generate a divergence proportional to the Galileon operator G 2 , so the bare bias parameters b G 2 . This means that, even if the non-local bias parameter is set to zero at some scale Λ-as in the case of the local Eulerian biasing model-this will no longer be true at some other scale Λ . Consequently, the local Eulerian biasing model is not a consistent model beyond the tree-level approximation. Finally, we note that the dependence of b (R) G 2 on the bare bias parameters is quite different from the result of [16], where the time-evolved Eulerian bias parameter, b G 2 , was related to the linear Lagrangian bias parameter, b

Halo Statistics
Self-consistent renormalization has forced us to consider a biasing model of the form where the right-hand side is a double expansion in small fluctuations (Φ g and Φ v ) and spatial derivatives (∇ i ). After renormalization, the higher-derivative terms are suppressed by Λ , the physical scale of non-locality in halo formation. For dark matter halos, we expect 13 Λ −1 ≤ k −1 NL , but ultimately Λ should be determined from N-body simulations or observational data. The renormalized biasing model (4.1) should be viewed as an effective theory valid on scales larger than Λ −1 . As is typically for effective theories, only a finite number of terms need to be retained in (4.1) in order to describe halo statistics to a finite accuracy. In §4.1, we perform a simple power counting to estimate the relative sizes of the renormalized operators. In §4.2 and §4.4, we describe how the lowest-order bias parameters can be measured by fitting the predictions of (4.1) to a variety of halo correlation functions (either in N-body simulations or in observational data). Readers who are less interested in the technical details may find summaries of results in §4.3 and §4.5.

Power Counting
We begin with an estimate of the relative sizes of the renormalized operators. In the absence of derivatives of δ, we may simply count powers of the linear dark matter density contrast δ (1) . The relative contribution from higher-derivative terms depends on the initial statistics. 14 In [22,23], the contribution from derivatives were estimated using a power law ansatz for the dark matter spectrum where n is scale-dependent, varying from n −2.1 [23] near the non-linear scale (at z = 0) to n ≈ 1 on large scales. Since the higher-order biasing terms are mostly relevant near the non-linear scale, we will use n = −2 for our estimates. In that case, 1 2 (3+n) = 1 2 and each power of δ roughly adds a power of q 1/2 . To be conservative about the relevance of higher-derivative operators, we will use Λ ∼ k NL . In reality, we expect Λ > k NL and higher-derivative contributions will be more 13 Recall that dark matter particles have travelled less than the non-linear distance k −1 NL over the history of the universe [26]. 14 In the EFT-of-LSS, it has been shown that higher-derivative terms become important as one approaches the non-linear scale [21].
suppressed than what we estimate here. Adopting this power counting for suitably renormalized operators, we get where in the last equality we have used Λ ∼ k NL and 1 2 (3+n) ∼ 1 2 . At leading order in derivatives, we therefore have Shown here is the leading q-scaling in the long-wavelength limit. When these operators are inserted into correlation functions they may be further suppressed due to the fact that certain correlations vanish for Gaussian initial conditions. Higher-derivative operators enter at higher order in the q/k NL expansion: Of course, when the operators in (4.4)-(4.7) are inserted in the bias expansion (4.1) their relative contributions will depend on the relative sizes of the renormalized bias parameters b (R) O . A large hierarchy between certain bias parameters can affect the estimates that we have performed here. Hence, our estimates should only be viewed as qualitative guidelines, and a more detailed treatment (analogous to [23]) is clearly required for comparisons with real data.

Two-Point Statistics
In this section, we compute the halo-matter power spectrum, P hm (q) ≡ (δ h ) q δ q , and the halo-halo power spectrum, P hh (q) ≡ (δ h ) q (δ h ) q , at one loop (i.e. to fourth order in δ (1) ). We will present results only at leading order in derivatives, but it will be clear how higher-derivative operators would be included. The (renormalized) operators which give non-vanishing one-loop contributions to P hm and P hh then are: (4.8) The non-linear dark matter power spectrum, P mm (q) ≡ δ q δ q , is renormalized in the EFT-of-LSS, cf. eq. (2.32). The contributions from the remaining operators [O] will be discussed in this section.

Halo-Matter
Diagrammatically, the correlation between [O] and δ is We have separated [O] q δ q into F-terms, which contain a contraction between the two legs of the composite operators, and I-terms, which only contain contractions with the external leg. The halo-matter power spectrum (4.10) can then be written as [14] P hm (q) = b (R) where we have defined b (R) where σ 2 p,q ≡ (p · q/pq) 2 − 1. The two non-zero F-terms are proportional to each other (4.14) The two operators G 2 and Γ 3 therefore contribute degenerate F-terms to the halo-matter power spectrum. Since the F-terms are proportional to the linear power spectrum they be interpreted as a scale-dependent contribution to the linear bias b (R) In the case of a scaling universe with P (q) ∝ q n , the finite part of F can be computed analytically using dimensional regularization [22] f where c(n) is a coefficient which depends on the scaling of the power spectrum. In our universe, the integrals in (4.12) and (4.13) are convergent as the cut-off Λ is taken to infinity. Of course, by sending Λ to infinity one is including contributions from scales which are outside the regime of validity of the effective description. However, it is easy to see that these finite errors can be absorbed into the bias parameters of higher-derivative operators. After renormalization, we expect the higher-derivative operators to be suppressed by a (momentum) scale that is larger than the scale of variation of the functions in (4.12) and (4.13).
• I-terms.-The functions I [O] in (4.10) are Although the integrals in (4.17) and (4.18) receive contributions from all scales, they are finite up to contact terms. More precisely, these integrals can be written as the sum of a finite term and a cut-off dependent term, Importantly, the functions J [O] (q) are analytic in q 2 , i.e. they can be written as expansions in powers of (q 2 ) n . In position space, these terms become derivatives of delta functions and therefore disappear when correlation functions are evaluated at separated points. The cut-off dependent parts are therefore contact terms which can be safely discarded and only the physical finite terms I [O] (q) are kept. In a scaling universe, the finite parts can again be computed in dimensional regularization where d [O] (n) are coefficients which depend on the scaling n of the power spectrum. In the real universe, the two functions I [δ 2 ] and I [G 2 ] are approximately proportional to each other, as illustrated in fig. 3. The operators δ 2 and G 2 therefore contribute degenerate I-terms.

Halo-Halo
A similar one-loop calculation for the halo-halo power spectrum gives [14] P hh (q) = b  where the functions As before, these integrals are finite up to contact terms. In a scaling universe, these finite parts simply are (4.28) In the real universe, the two functions I [δ 2 ,G 2 ] and I [G 2 ,G 2 ] are approximately proportional to each other, as illustrated in fig. 4. The function I [δ 2 ,δ 2 ] is divergent, but the finite piece can be evaluated by subtracting the q = 0 divergent part, This expression does not remove (possibly divergent) contact terms of the form q 2n for integer n ≥ 1, which are derivatives of delta functions in position space. These terms can be removed by hand without altering correlation functions at finite separation. This freedom to alter contact terms can be made explicit by introducing a stochastic bias parameter (see e.g. [13]). We will simply ignore them with the understanding that our expressions are correct up to such contact  terms. For the real universe, these contact terms are known to be highly suppressed in the quasi-linear regime [21,26].

Summary: Effective Bias Parameters
We have arrived at the following predictions for the halo power spectra L includes contributions from δ, G 2 and Γ 3 , while the sums over I-terms are only over the operators δ 2 and G 2 . Fitting (4.31) and (4.32) to data over a sufficiently wide range of momenta, in principle, allows the effective bias parameters b In that case, measurements of the power spectra only determine two linear combinations of these three bias coefficients. The real universe is sufficiently close to a scaling universe (in the regimes of interest), that these problems may persist. In fact, taking into account the near-degeneracy of the I-terms, we can write the halo-matter power spectrum as where i(q) ≡ I(q)/P (q) and The functions f (q) and i(q) are sufficiently different that it should be possible to measure the effective bias parameters b F and b I independently. However, a degeneracy between the parameters b (R) 2 the halo-halo power spectrum,

Higher-Point Statistics
An alternative to break the degeneracy between the different bias contributions is to consult higher-point statistics. In this section, we present an analysis of the bispectrum 36) and the trispectrum T hmmm (q, q 1 , q 2 , q 3 ) ≡ (δ h ) q δ q 1 δ q 2 δ q 3 . (4.37) The bispectrum will be evaluated at one loop, while the trispectrum will only be presented at tree level. Details of the bispectrum calculation can be found in Appendix D.

Three-Point Function
At tree level, the bispectrum (4.36) gets contributions from δ, δ 2 and G 2 : where µ 12 ≡ q 1 · q 2 /q 1 q 2 and B mmm (q, q 1 , q 2 ) ≡ δ q δ q 1 δ q 2 is the dark matter bispectrum. We see that the contributions from δ 2 and G 2 come with distinct momentum dependences, which has been used previously to fit for the parameters b At one-loop, operators up to quartic order contribute to (4.36). However, composite operators with four legs, such as δ 4 , do not contribute after renormalization. We therefore only have to consider operators with up to three legs. In total there are eleven such operators, including all the operators in (4.8). These operators are discussed in detail in Appendix D. Here, we collectively call them O. The one-loop diagrams can be organized into two classes of contributions: • F-terms.-Diagrams whose only contractions are between the internal legs of the composite operator are proportional to the dark matter power spectra where the functions f B (q 1 , q 2 ) are defined explicitly in Appendix D.
• I-terms.-Diagrams with at least one contraction with an external leg contain non-trivial convolutions between the dark matter power spectra and the kernel functions of standard perturbation theory, Six of the eleven operators have such contributions, which are presented explicitly in Appendix D.
In the end, the one-loop bispectrum can be written as Due to the complexity of the final answer, the one-loop contributions to the bispectrum are probably of limited use in determining the effective bias parameters. However, knowledge of the functional form of the one-loop corrections gives us a handle on the expected theoretical error in the tree-level ansatz (4.38).

Four-Point Function
We have seen that the one-loop power spectra in §4.2 involve contributions proportional to b (R) Γ 3 . This parameter can be measured from the tree-level trispectrum (4.37), which receives contributions from every operator up to cubic order If the linear and quadratic bias parameters are measured from the power spectrum and the bispectrum, then the trispectrum can be used to extract the cubic bias parameters O 3 . At leading order in perturbation theory, we have where the functions g [O 3 ] contain the specific momentum dependence associated with each operator in O 3 : We see that the different contributions can be distinguish on the basis of their unique momentum dependences in the trispectrum. In particular, it is possible to extract information about the velocity tidal tensor Γ 3 (Φ g , Φ v ) and measure the associated bias parameter b

Conclusions
In this paper, we have shown explicitly how renormalization forces us to treat biasing models as effective theories, i.e. as a double expansion in terms of fluctuations in the dark matter density and velocity and derivatives thereof. Consistently removing the short-scale physics in the local Eulerian biasing model [6] doesn't just renormalize the bias parameters, but also generates nonlocal terms and higher-derivative contributions. In order for the theory to become independent of the unphysical regulator of composite operators, all terms consistent with the symmetries have to be included in the biasing model. At lowest order, this means adding the gravitational tidal tensor, while at cubic order, the velocity potential appears as an independent degree of freedom.
In the process, we have clarified a few technical aspects of the renormalization procedure proposed by McDonald [13]. We organized the renormalization of composite operators in a convenient diagrammatic representation and derived the building blocks for the one-loop renormalization of the halo power spectrum and bispectrum. We proved that Galileon operators aren't renormalized at leading order in derivatives. Finally, we showed explicitly how the definition of the renormalized theory depends on the renormalization scale and how the terms in the bias expansion mix as this scale is varied. This scale-dependence of the renormalized halo bias is relevant for interpreting recent N-body results [24,34,39].
Our work motivates several future directions: • Most importantly, it remains to be quantified how many terms need to be kept in the effective theory to achieve a given target of precision in the predictions for the statistics of halos. The importance of including the gravitational tidal tensor has been established through N-body simulations in [16,17]. However, a completely systematic exploration of the effective theory of halo biasing has not yet been performed. For example, the scale Λ which suppresses higher-derivative terms has not been measured and the errors that arise from truncating the effective theory have not been quantified.
• A key motivation for understanding non-linear biasing is primordial non-Gaussianity. One may hope that the non-linearities arising from the biasing are sufficiently distinct, so that the primordial signals can be extracted from the shape information of the correlation functions. To analyze this self-consistently, we must consider the effects of these non-Gaussian contributions on the renormalization of the biasing model [18,19].
• We have focused on Eulerian biasing, but it has been suggested that biasing in Lagrangian space may have some advantages. The approach taken here can likely be adapted to the Lagrangian EFT-of-LSS [28]. Which scheme is more useful may depend on the type of observable that is considered.
• We have shown how renormalization forces biasing to be non-local in space, but the theory has remained local in time. Recently, the necessity of non-locality in time has been emphasized in the EFT-of-LSS [27,29]. This deserves further consideration in the present context.
• Ultimately, galaxies are observed in redshift space. Redshift space distortions can add extra non-linear contributions to the observed galaxy correlation functions. This should be taken into account.
• We have seen that a biasing model in terms of just the dark matter density is not consistent and that the velocity potential has to be added as independent degree of freedom. It would be interesting to explore the physical effects of the velocity potential more widely, e.g. beyond perturbation theory.

A Scale Dependence and Operator Mixing
In the Wilsonian renormalization scheme that we have adopted in this paper, the bare bias parameter will depend on the cutoff, 15 b (0) O (Λ), in a way that is dictated by the renormalization group equations db Although the Λ-dependence of bare quantities isn't physical, it is often an indication that the corresponding physical quantities depend on the scale at which they are measured. This is typically reflected in the Callan-Symanzik equation applied to correlation functions. In this appendix, we will show that this intuition is correct, although the use of the Callan-Symanzik equation will not be necessary. 16 In Section 3, we defined the renormalization conditions in the long-wavelength limit (q → 0), where non-linearities are negligible (δ → δ (1) ). At finite separation, 17 |x| = 2π/µ, the basis of renormalized operators, [O I ], will therefore, in general, not be orthogonal where Γ IJ (µ, k NL ) is a real positive-definite matrix. 18 Of course, we can diagonalize Γ IJ in terms of a new basis of operators, The basis of operators that is orthogonal at a scale µ will mix and in general won't be an orthogonal basis at another scale µ . In order for the halo density contrast to be independent of µ, the renormalized bias parameters need to be µ-dependent and the bias expansion is 15 It would be interesting to relate this Λ-dependence to the dependence of scatter-plot bias parameters on the smoothing scale that is observed in N-body simulations. 16 Renormalization group flow is usually most useful for logarithmic divergences. For a scaling universe with P (q) ∝ q n , the variance σ 2 (Λ) is logarithmically divergent only for n = −3, in which case there are also infrared divergences. For this reason, we will find this language less useful. 17 We chose to work in position space to avoid contributions that are delta-function localized. However, this means that our definition of δ h is only correct up to contact terms. The above procedure can also be applied in momentum space with q = µ, if one is careful to identify only the terms that are non-analytic in q. 18 Here, we are ignoring "descendants" of O, i.e. operators which are total derivatives of O. Since these terms are fully correlated with O, they cannot be diagonalized. 19 A similar approach to measuring bias parameters was taken in [39], although it wasn't applied to the renormalized theory.
What are we to make of the µ-dependence of the bias parameters? Let us first consider the case where the correlation functions of δ are purely Gaussian. In this case, we can define [δ n ] = H n (δ) [40], where H n are Hermite polynomials, such that where ξ(|x|) ≡ δ (1) (x)δ (1) (0) . In terms of the Hermite polynomials, the bias coefficients b (R) n are independent of µ. This basis is diagonal at all scales. The µ-dependence of the bias coefficients must therefore arise from non-linear evolution. As a simple example, let us consider the mixing between the operators [δ] and [δ 2 ] as a function of separation. The detailed form of this mixing was computed in §2.2.1. In a scaling universe, where P ∝ q n , dimensional analysis shows that to fourth order in δ (1) , we have where ε ≡ (µ/k NL ) 3+n . To order ε 2 , the orthogonal basis of operators is As a result, we get the following relationship between the bias parameters measured at µ and µ (dropping mixing with higher-derivative terms) We see that for n −3, the scale dependence of the bias parameters vanishes rapidly as we approach the linear regime, µ k NL . Of course, there can still be important mixing with higherderivative terms even for large values of n. However, the measurement of higher-derivative bias terms is more complicated as the bias parameters cannot always be diagonalized.

B Renormalization of Higher-Order Operators
In this appendix, we describe in more detail the renormalization of operators of third order and higher.

B.1 Renormalization of δ n
First, we show that the one-loop renormalization of δ n , with n > 2, is completely determined by the one-loop renormalization of δ 2 . More precisely, at one-loop, we have where the superscript "loop" indicates that two (and only two) linear δ (1) 's are contracted inside of δ n , while the superscript "tree" denotes that the δ's of δ n−2 are contracted only with the external legs δ (1) . The proof of the statement (B.1) is quite straightforward. Consider the correlation of δ n with E external legs. There are n(n−1) 2 ways of contracting two of the δ's in δ n to form a loop. The external legs are then contracted with either the δ 2 forming the loop or with the n − 2 remaining δ's. The following is a diagram in which r ≤ E external legs are contracted with δ 2 , while the rest are contracted with δ n−2 : The full correlation function is then obtained by summing over r : which implies (B.1). This result relates (at one loop) the counterterms of δ 2 and those of δ n For example, we can infer the one-loop counterterms for the renormalization of δ 3 and δ 4 directly from (2.44), where we have only shown terms up to order (δ (1) ) 4 .

B.2 Renormalization of G 2 δ
In the main text, we argued that the one-loop renormalization of the linear bias requires the renormalization of all cubic operators. However, the renormalization of the operator δ 3 is related to the renormalization of δ 2 (see §B.1) and the Galileon operators aren't renormalized at leading order in derivatives (see Appendix C), so the only operator which remains to be renormalized explicitly is G 2 δ. Since the renormalization of this operator is identical to the renormalization of δ 2 , we just state the final result

B.3 Renormalization of Quartic Operators
Finally, we renormalize the quartic operators relevant for the discussion in §3.2 up to m = 2. The operator δ 4 has already been renormalized in §B.1, while the Galileon operators are only renormalized by higher-derivative operators. The renormalization of G 2 δ 2 follows straightforwardly from the identity The remaining operators which need to be renormalized are G 3 δ, Γ 3 δ and [G 2 ] 2 . The counterterms required to cancel the divergences arising from these operators are:

C Non-Renormalization of Galileon Operators
In this appendix, we prove a non-renormalization theorem for the Galileon operators G n . In particular, we will show these operators aren't renormalized at leading order in derivatives. This result is similar to the non-renormalization theorem of Galileons in modified gravity [41,42].

C.1 Non-Renormalization of G 2
At zeroth order in derivatives, the quadratic Galileon operator 2 does not get renormalized. Loops arising from ∇ i ∇ j Φ∇ i ∇ j Φ are exactly canceled by (∇ 2 Φ) 2 and we get To prove this result, we only need to show that First, let us rewrite G 2 (Φ) in a more convenient way which, in Fourier space, becomes To complete the proof, we only need to demonstrate that correlation functions of [ΦD ij Φ] q with the linear density contrast (δ (1) ) m are't singular in the limit of vanishing the external momenta q i . In other words, we need to show that q m (p)1PI = P 1 · · · P m + · · · . (C.5) We first consider the case m = 1 and construct a proof by contradiction. Let us assume that the correlation function has the following singularity in the limit q 1 → 0, and show that this leads to a contradiction. The q −2 1 divergence in (C.6) can only be absorbed if one introduces the potential Φ as a counterterm (C.7) We now show that this term is not consistent with the symmetries. Indeed, the renormalized operator [ΦD ij Φ] can be written as where O ij are operators which depend locally on Φ and its derivatives. These operators need not be invariant under a constant shift and a homogeneous boost, since ΦD ij Φ is not. However, we expect that the renormalized operators satisfy the same symmetries as the bare ones. In particular, shifting the potential Φ by a constant c on both sides of (C.8), we get Since the potential is not renormalized, i.e. [D ij Φ] = D ij Φ, we find Comparing eqs. (C.8) and (C.10), we get Since the operators O ij form a basis, this is satisfied, if and only if As a result, the parameter a(Λ) in (C.7) necessarily vanishes and (C.14) m > 1 Let us now consider the general case m > 1. We assume that in the limit |q 1 + · · · + q p | → 0, with p ≤ m, the correlation function behaves as The only counterterm which could remove this divergence is δ m−p Φ. By the same logic as before, such a counterterm violates the symmetries of the problem. The divergence in (C.15) therefore cannot be present and we have q m (p)1PI = P 1 · · · P m + · · · . (C.16)

C.2 Non-Renormalization of G n
Finally, we prove that the non-renormalization theorem holds for every Galileon operator G n . The proof proceeds by induction. We will assume that 17) and prove that this holds for the n-th order Galileon operator.
Definition.-The n-th order Galileon operator can be written as where the tensor η is defined as 19) and the sum runs over n! permutations σ and (−1) σ represents the signature of the permutation. From the definition, it is clear that can therefore be written as Counterterms of ΦT ] q with (δ (1) ) m is not singular in the limit where the external momenta, or a partial sum of these external momenta, vanish. To show this, we will prove that symmetries forbid the appearance of counterterms proportional to the operator δ ij Φ. As before, we write the renormalized operator as where again the operators O ij are local functions of the potential and its derivatives. Shifting the potential Φ by a constant c, we get and, hence, . We see that if there is a counterterm proportional to O ij = δ ij Φ, the renormalization of T (n−1) ij would contain an operator proportional to δ ij . As we will show next, this cannot be the case.

D One-Loop Bispectrum
In this appendix, we give explicit expressions for the functions F B = f B P 1 P 2 and I B , which appear in the one-loop expression of the halo-matter-matter bispectrum Operators up to quartic order may contribute to the one-loop bispectrum. Classified by their number of legs, these operators can be grouped into three categories where we have introduced the quartic operators and the third-order Galileon operator with ∇ ij ≡ ∇ i ∇ j and (α, β, γ) ∈ {g, v}. Furthermore, at quartic order, the velocity divergence θ ≡ ∇ 2 Φ v becomes an independent degree of freedom which gives non-trivial contributions to the bispectrum at one-loop. This is captured by the operator [16] which vanishes up to third order.