\section{Introduction}
\label{sec:intro}

\hskip 0.8cm Over the last decade, we have accumulated a good deal of observational evidence that 
the large scale structure (LSS)   originated from seed fluctuations in the very early universe, 
which grew  via gravitational instability at later times. 
In this picture, galaxies form within Dark Matter (DM) halos \cite{whiterees}, and the  goal of 
the theory of bias \cite{kaiser84} is to relate observable properties of tracers such as the density 
contrast of galaxies to the underlying matter distribution and, ultimately, to the initial conditions.

DM  halos are, to a large extent, the building blocks of what we observe. Knowing  how DM  halos 
form and evolve under the action of gravity is crucial to understand  the environment that 
harbors galaxies and clusters. Understanding the clustering properties of halos
represents, therefore, a  key ingredient towards an accurate description of galaxy clustering 
statistics. 

The aim of this paper is to write down an Effective Boltzmann Equation (EBE) for the single phase 
space density of halo centers, which can be used to derive {\it ensemble average} halo clustering
statistics. Therefore, this must be an effective description, in which the halo distribution 
should be thought of as being the mean-field distribution associated with a given realization of the DM  
distribution. 
This is the correct interpretation of the halo overabundance defined through a perturbative bias expansion 
so long as the surveyed volume remains finite, as is the case of real or simulated data
\cite{dcss,Lazeyras:2015giz}.

Since the halo Boltzmann  equation will describe the dynamics of the halo mean-field, one should expect 
to recover exactly the Boltzmann equation for DM only if the halo centers are statistically unbiased 
relative to the coarse-grained DM density field. 
This is assumed to be true in the coupled-fluids approximation for the coevolution of DM and halos, 
which is largely used to compute the time evolution of bias by following the evolution over cosmic 
time and in Eulerian space of the halo progenitors - the so-called proto-halos - until their 
virialization \cite{fry,teg,elia,chan0}. 
In this framework, the momentum equations for halos and DM do not differ at all. However, this 
approximation is not valid in general, as one should expect that the force acting on large scale 
structure tracers is indeed biased relative to the force acting on a random field point. 
In other words, the sampling bias induced by the tracers propagates into their dynamical evolution.

Therefore, the EBE for halos will be different from that of DM.  
This point should not be underestimated: the DM Boltzmann equation is the starting point of
any analytical approach to the description of the  LSS (see e.g. \cite{davispeebles}).
Here we advocate that, since the ultimate goal is the understanding of the statistical 
properties of what we observe (that is galaxies), one ought  to adopt the EBE discussed below as the 
starting point for halo centers.

This is not the only reason for pursuing this project, though. 
LSS  probes such as peculiar velocities or redshift space distortions can be used to test the law of 
gravity through a measurement of peculiar velocities, i.e. deviations from pure Hubble flow 
\cite{vjd,kaiser87}. 
However, peculiar velocities are effectively measured at the position of halo centers. Therefore, even 
if the latter locally flow with the DM, velocity statistics may be biased if the halo peculiar velocities 
do not provide a fair sample of the matter flows 
(for recent discussions  see Refs. \cite{jennings,zzj1,jennings2,yzzj,jennings3,chiang}).
Thus far, measurements of the halo velocity power spectrum appear consistent with little or no 
statistical velocity bias \cite{pw,zzj1}. However, \cite{v1,jennings} reached somewhat different 
conclusions. Furthermore, the effect must vanish in the limit $k\to 0$ so that, at low redshift, 
scale-dependent nonlinearities complicate the interpretation of the results.

Peak theory predicts the existence of such a statistical effect which, in 
linear theory, manifests itself as a reduced halo velocity dispersion 
\cite{bbks} together with a $k$-dependent velocity bias in the usual Kaiser formula
\cite{v0} and in the time evolution of the linear peak bias \cite{dcss}.
The authors of Ref. \cite{v1} have recently detected through N-body simulations a 
large scale halo velocity bias $b_v(k)$ 

\begin{equation}}   \def\ee   {\end{equation}
\label{biasv}
b_v(k)=\bigl(1-R_v^2k^2\bigr),
\ee
which appears to remain constant throughout time until virialization. Here, 
$R_v$ is the typical scale of the halo velocity bias

\begin{equation}}   \def\ee   {\end{equation}
R_v^2=\frac{\sigma_0^2}{\sigma_1^2},\,\,\,\,\sigma_j^2=\int\frac{\d^3 k}{(2\pi)^3}\, k^{2j} 
P(k) \tilde{W}^2(kR),
\ee 
$P(k)$ is the DM power spectrum and $\tilde{W}(x)$  is a spherically symmetric, time-independent 
kernel with filtering scale $R$ equal to the Lagrangian halo radius.
This result confirms that, 
{\it when the purpose is to compute statistical quantities such as power spectra, etc}, 
the coupled-fluid approximation requires a modification of the momentum conservation equation for 
halos  \cite{v1} as   the  force felt by the halos  is  biased relative to the force acting on a 
random DM particle.

The EBE  provides a  way of describing the evolution equation for biased tracers. 
As first suggested in Ref. \cite{v2}, it encodes the  statistical velocity bias which originates 
from the fluctuations  around the ensemble average of the gravitational force acting on the halo 
and it predicts the right value (\ref{biasv}) observed at the linear level. 
Therefore, the EBE explains at  more fundamental level why and how one needs to modify the momentum 
fluid equation  for the coupled-fluids of halos and DM, as suggested by Ref. \cite{v1}. 
In addition, it allows us to compute the halo velocity bias at higher orders and to derive additional 
stochastic components that describe in a {\it mean-field} way (i.e. at the statistical level) how 
biasing influence the dynamics of halo centers.
 
The paper is organised as follows: in Sec.~\S\ref{sec:boltzmann equation}, we motivate our approach 
and introduce the various averaging we perform throughout the calculation.
In Sec.~\S\ref{sec:boltzmann eq LHS} we perform a detailed computation of the left-hand side 
of the EBE adopting  a general method based on the path integral approach for computing conditional 
averages. 
Section~\S\ref{sec:boltzmann eq RHS} is devoted to a similar discussion regarding the right-hand 
side of the EBE.
We then conclude with some final remarks in Sec.~\S\ref{sec:conclusions}.


\section{From the Klimontovich-Dupree equation to the Boltzmann equation}
\label{sec:boltzmann equation}

\hskip 0.8cm
Let us first explain in more details what we mean by halo ``mean-field''.
In a local bias approximation (e.g. \cite{Fry:1992vr}), one can express density fluctuations in the halo mean-field as the series expansion 
\begin{equation}
\delta_\text{h}({\vec r},t)=b_1 \delta_R({\vec r},t) + \frac{1}{2} b_2 \big(\delta_R^2({\vec r},t)-\sigma_0^2(t)\big) 
+ \dots
\end{equation}
where the filtering scale $R$ is proportional to the halo mass \cite{v0}, the bias parameters $b_1$,
$b_2$ etc. are {\it ensemble average} whereas $\delta_R({\vec r},t)$ is {\it one} particular realisation of 
density fluctuations in the large scale structure at cosmic time $t$. The filtering reflects the fact
that we are not interested in the internal structure of the halos, only in their center-of-mass 
position and kinematics. From the knowledge of this halo mean-field, we can compute halo clustering 
statistics, such as the halo 2-point correlation, as
\begin{equation}
\xi_\text{h}({\vec x}_1-{\vec x}_2,t) = b_1^2 \big\langle \delta_R({\vec x}_1,t) \delta_R({\vec x}_2,t)\big\rangle + \dots
\end{equation}
where $\langle\dots\rangle$ is an average over realisations of the large scale structure.


Our goal is to write down an EBE  for the single particle  phase space density $f_\text{h}$ of halo 
centers, such that its zeroth moment yields $\delta_\text{h}({\vec r},t)$. 
This EBE must be  different from the Boltzmann equation of DM: the latter is recovered  only if the 
halos are statistically unbiased relative to the coarse-grained DM density field (which they are not).
Our starting point is  the Klimontovich density \cite{klim}
\begin{equation}}   \def\ee   {\end{equation}
f_{\rm K}({\vec r},{\vec v},t)=\sum_i\delta_{\rm D}\left[{\vec r}-{\vec r}_i(t)\right]\delta_{\rm D}\left[{\vec v}-{\vec v}_i(t)\right],
\ee
which is the single particle phase space density for DM. We have assumed that all the elementary DM particles have equal mass $m$, and we normalize our units so that $m=1$. 

The corresponding equation of motion for such a phase space density is the Klimontovich-Dupree equation
\begin{equation}}   \def\ee   {\end{equation}
\label{kd}
\frac{\partial f_{\rm K}}{\partial t}+\vec{v}\cdot \frac{\partial f_{\rm K}}{\partial{\vec r}}
-\vec{\nabla}\Phi_{\rm K}\cdot\frac{\partial f_{\rm K}}{\partial{\vec v}}=0,
\ee
where 
\begin{equation}}   \def\ee   {\end{equation}\label{phik}
\vec{\nabla}\Phi_{\rm K}({\vec r},t)=G_{\rm N}\int\d^3r'\d^3v'f_{\rm K}(\vec{r}^{\, \prime},\vec{v}^{\, \prime},t)\frac{({\vec r}-\vec{r}^{\, \prime})}{\left|{\vec r}-\vec{r}^{\, \prime}\right|^3}
\ee
is the corresponding gravitational force. 
At this stage we have not learned much, as the Klimontovich density encodes the trajectories of all the particles in the system, which is much more information than we can handle and want to keep track of (in particular, we are not interested in the internal properties of virialized structures). Moreover, eq.~\eqref{kd} describes one particular realization of a universe, which may not resemble our own Universe. 
Therefore, we go one step further and average over realizations of universes that have the same large-scale structure (i.e. similar coarse-grained properties). In this way we can construct macrostates by averaging over a statistical ensemble of microstates with similar phase space density in (small) volumes containing a sufficient amount of particles.\\ 
Let us explain in more detail how this ensemble average is performed with the help of fig.~\ref{fig:micro}.
\begin{figure}[h!]
\centering
\begin{tikzpicture}[node distance=1cm, auto,]

\node[rectangle, rounded corners,draw=black, very thick,
text width=8.7cm,minimum height=8cm,
align=center] (microbox) {
  { \textbf{(a)} $f_\text{K}(\vec r,\vec v,t)$}\\
  \vspace{1em}		{\pgfuseimage{micro}}\\
  \vspace{1em}		{\hspace{4em} smearing}\\
  \vspace{1em}		{\pgfuseimage{smearedSmall}}\\
  \vspace{1em}		{ \textbf{(b)} $P_{\vec r}(\rho,\vec v,t)$}
};
\node[above=0.cm of microbox] (microtitle) {{\large {Microstates}}};
\node[above=-3.4 of microbox] (c1){};
\node[above=-4.9cm of microbox] (c2){};
\node[] {} (c1.south) edge[freccia,bend right=30](c2.north);

\node[rectangle, rounded corners,draw=black, very thick,
text width=3.7cm,minimum height=8cm,
align=center, inner sep=5pt,right=2.2cm of microbox] (macrobox) {
  { \textbf{(c)} $f(\vec r,\vec v,t)$}\\
  \vspace{1em}		{\pgfuseimage{macroSmall}} 
};
\node[above=0.cm of macrobox] (macrotitle) {{\large {Macrostate}}};

\node[right=-0.2cm of microbox](r){};
\node[above=2cm of r](rUp){};
\node[below=2cm of r](rDown){};
\node[left=-0.2cm of macrobox](l){};
\node[above=0.1cm of l](lUp){};
\node[below=0.1cm of l](lDown){};
\node[] {} (rUp.east) edge[freccia,bend right=45](lUp.west);
\node[] {} (rDown.east) edge[freccia,bend left=45](lDown.west);
\node[below right=0.cm of rUp]{{ \small $\langle \cdots\rangle_\text{micro}$}};
\node[above right=-0.1cm of rDown]{{\footnotesize $\int \hspace{-3pt}\mathrm d\rho \rho P_{\vec r}(\rho,\vec v,t)$}};
\end{tikzpicture}
\caption{A pictorial exemplification of the ensemble average described in the text. We consider all the configurations (\textit{microstates}) for the distribution of DM particles that yield to indistinguishable large scale structure configurations (a given \textit{macrostate}). 
\textbf{(a)} The distribution of DM particles is described in the phase space by the Klimontovich density $f_\text K({\vec r},{\vec v},t)$. Four different realizations are shown. 
\textbf{(b)} After a suitable smearing (the convolution with a window function) we can describe these microstates in terms of a smooth density $\rho$. In any point ${\vec r}$, $\rho$ can vary in a range of values among the microstates: this stocasticity is encoded in the probability distribution $P_{\vec r} (\rho,{\vec v},t)$. 
\textbf{(c)} The macrostate is described by the phase space density $f({\vec r},{\vec v},t)$, which follows from the ensemble of microstates through the ensemble average $\langle f_\text K ({\vec r},{\vec v},t)\rangle_\text{micro}$, or equivalently by computing the expected value of $\rho$ with respect to the probability density $P_{\vec r} (\rho,{\vec v},t)$. 
}
\label{fig:micro}
\end{figure}

The Klimontovich density describes a set of points in the phase space. Let us focus on a small box in the configuration space.
In the part \textbf{(a)} of fig.~\ref{fig:micro} we show four different realizations of this box, which yield the same LSS. Their ensemble average, which we denote by $\langle\cdots \rangle_{\rm micro}$, gives the mean value $f({\vec r},{\vec v},t)$:
\begin{equation}}   \def\ee   {\end{equation}
\label{f average}
\Big<f_{\rm K}({\vec r},{\vec v},t)\Big>_{\rm micro}=
\Big<\sum_i\delta_{\rm D}\left[{\vec r}-{\vec r}_i(t)\right]\delta_{\rm D}\left[{\vec v}-{\vec v}_i(t)\right]\Big>_{\rm micro}=f({\vec r},{\vec v},t)\;,
\ee

An alternative way of understanding this average is to smear out the discrete distribution of points into a continuous density with an associated probability $P_{{\vec r}}(\rho, {\vec v},t)$, where ${\vec r}$ is any position inside the box, as represented in the part \textbf{(b)} of fig.~\ref{fig:micro}. The range of values assumed by $\rho$ over the ensemble defines the probability density $P_{\vec r}(\rho,{\vec v},t)$ for having a velocity ${\vec v}$ and a density $\rho$ in a small volume centered at position ${\vec r}$, provided that the LSS is the one fixed in \textbf{(c)} in fig.~\ref{fig:micro}. \footnote{Note that this smearing should not be confused with the filtering procedure introduced in section \ref{sec:boltzmann eq LHS}. The latter defines a UV cutoff in momentum space denoted by $\Lambda$, while the former represents an average over microstates in small, real space volumes of size $\ell$, with the requirement that $\ell \ll 1/\Lambda$.}
In this case, $f({\vec r},{\vec v},t)$ is given by the mean value of $\rho$ with respect to the probability distribution $P_{\vec r}(\rho,{\vec v},t)$:
 \begin{equation}}   \def\ee   {\end{equation}
 \label{f integral rho}
 f({\vec r},{\vec v},t)=\int_0^\infty\,\d\rho\,\rho\, P_{{\vec r}}(\rho,{\vec v},t)=\int_0^\infty\,\d\rho\,\rho\, P_{{\vec r}}(\rho|{\vec v},t)\,P_{{\vec r}}({\vec v},t)=\langle\rho|{\vec v}\rangle\,P_{{\vec r}}({\vec v},t).
 \ee
We will drop the subscript ${\vec r}$ here and henceforth, as it is clear that $P(\rho,{\vec v},t)$ is defined at a given spatial location. 

The distribution function $f({\vec r},{\vec v},t)$ that characterizes the macrostates defines, upon taking moments of the velocity ${\vec v}$, the overdensity field, the bulk velocity, etc. of a given macrostate. For instance, the matter overdensity field $\delta({\vec r},t)$ is the zeroth moment
\begin{equation}}   \def\ee   {\end{equation}
\bar{\rho}\bigl(1+\delta({\vec r},t)\bigr) = \int\d^3v\, f({\vec r},{\vec v},t).
\ee
This distribution function satisfies the equation
\begin{equation}}   \def\ee   {\end{equation}
\label{akd}
\frac{\partial f}{\partial t}+\vec{v}\cdot \frac{\partial f}{\partial{\vec r}}
-\vec{\nabla}\Phi\cdot\frac{\partial f}{\partial{\vec v}}=\vec\nabla_{{\vec v}} \cdot\vec F,
\ee
where
\begin{equation}}   \def\ee   {\end{equation}\label{phi}
\vec{\nabla}\Phi=G_{\rm N}\int\d^3r'\d^3v'f(\vec{r}^{\, \prime},\vec{v}^{\, \prime},t)\frac{({\vec r}-\vec{r}^{\, \prime})}{\left|{\vec r}-\vec{r}^{\, \prime}\right|^3}
\ee
and 
\begin{align}
\label{force}
\vec{F}({\vec r},{\vec v},t) 
&={\rm Cov}\left[\vec\nabla\Phi_{\rm K}({\vec r},t),f_{\rm K}({\vec r},{\vec v},t)\right]_{\rm micro} \nonumber \\
&=  G_{\rm N}\int\d^3r'\, \d^3 v'\, \frac{({\vec r}-\vec{r}^{\, \prime})}{\left|{\vec r}-\vec{r}^{\, \prime}\right|^3} f_{\rm 2c}({\vec r},{\vec v},\vec{r}^{\, \prime},\vec{v}^{\, \prime},t)
\;,
\end{align}
where
\begin{equation}}   \def\ee   {\end{equation}
\label{irr}
f_{\rm 2c}({\vec r},{\vec v},\vec{r}^{\, \prime},\vec{v}^{\, \prime},t)=f_{\rm 2}({\vec r},{\vec v},\vec{r}^{\, \prime},\vec{v}^{\, \prime},t)-f({\vec r},{\vec v},t)f(\vec{r}^{\, \prime},\vec{v}^{\, \prime},t)
\ee
is the irreducible   two-particle correlation function. 
 
Equation (\ref{akd}) is our starting point towards a derivation of the EBE  for halos, as it enables us to apply constraints involving $\delta$ etc.~to define biased tracers of the large scale structure.


\section{Halos as extrema of the smoothed density field}
\label{sec:halos}
\hskip 0.8cm
In order to write down an EBE for the halo phase space density $f_\text{h}({\vec r},{\vec v},t)$, we follow the evolution over cosmic time, and in Eulerian space, of the so-called proto-halos until their virialization. The proto-halos are the progenitors of isolated DM halos, i.e. DM halos that are not contained in any larger halo.
While  their shape and topology change as a function of time (smaller substructures gradually merge to form the final halo), their centre of mass moves along a well-defined trajectory determined by the surrounding mass density field.
Hence, unlike virialized halos that experience merging, by construction proto-halos
always preserve their identity. Their total number is, therefore, conserved over time, such that we can write a continuity and Euler equation for their number density and velocity, respectively.

We will consider a specific model for the halo biasing in order to derive explicit expressions for the amplitude and scale-dependence of the effective corrections to the halo EBE.
Namely, we will hereafter designate by halos (or peaks) those clustered objects that are located at the points where $\vec{\nabla}\delta_R({\vec r})=\vec{0}$ and where the smoothed density contrast satisfies $\delta_R({\vec r})=\nu\sigma_0$ (we will occasionally write $\delta_\text{pk}\equiv\nu\sigma_0$).
The subscript $R$ represents quantities that have been smoothed on the (comoving) size $R$ of a halo. For instance,
\begin{equation}
\delta_R({\vec r},t)= \int \d^3 r' \delta(\vec{r}^{\, \prime},t)\, W\!\left(\frac 1R |{\vec r}-\vec{r}^{\, \prime}|\right)\;,
\end{equation} 
where $W(x)$ a spherically symmetric smoothing kernel. 
Moreover, $\nu$ is the peak height and $\sigma_0\equiv \sigma_0(t)$ is the variance of the DM field at time $t$ in a given region of size $R$. In general, the spectral moments of the smoothed fields are given by
\begin{equation}}   \def\ee   {\end{equation}
\label{sigma def}
\sigma_j^2=\int\frac{\d^3 k}{(2\pi)^3}\, k^{2j} 
P_{\rm }(k)\tilde{W}^2(kR),
\ee
where $\tilde{W}(x)$ is the Fourier transform of $W(x)$,  $P(k)=\big\langle |\delta({\vec k},t)|^2\big\rangle_\text{macro}$ is the matter power spectrum, the ensemble average is over the macrostates and the time dependence is implicit. We will occasionally refer to those conditions as the \textit{peak constraint}, although we enforce the extremum condition solely. However, notice that, in the high peak limit $\nu\gg 1$, nearly all extrema are maxima of $\delta_R$.

Last but not least importantly, the filtering of the density field reflects the fact that halos are extended objects. In our model, we will take this into account upon assuming that halo centers feel a force which is a smoothed version of the force acting on dark matter particles.


\section{The halo Boltzmann equation: the left-hand side}
\label{sec:boltzmann eq LHS}
\hskip 0.8cm

Our goal is to specialize the Boltzmann equation \eqref{akd} from the DM phase space density to the peak phase space density so as to describe the evolution of halos. For this purpose, $f_\text{pk}(\vr,\vv,t)$ will designate the peak phase space density -- i.e. the number density of peaks of velocity ${\vec v}$ at the position ${\vec r}$ and time $t$ --
 such that its zeroth moment corresponds to the halo mean-field $\delta_\text{h}({\vec r},t)$.
The density $f_\text{pk}$ can be written in terms of conditioned probabilities, as we will show in eq.~\eqref{fpk} of section \ref{sec:boltzmann eq RHS}. 
The peak constraint in ${\vec r}$ also affects the computation of the force on the right-hand side of Boltzmann equation, 
introducing conditioned probabilities. 
As we will discuss in section \ref{sec:boltzmann eq RHS}, this side of the Boltzmann equation eventually vanishes.

Let us focus now on the left-hand side of \eqref{akd}. There, the following operator appears: 
\begin{equation}}   \def\ee   {\end{equation}
\label{composite operator}
\vec{\nabla}\Phi(\vr,t)\cdot\partial_{\vec v} f (\vr,\vv,t) \,,
\ee
which is the product of two operators evaluated at the same point, that is, a \textit{composite operator}. 
It is well-known that, owing to the  local product of fields making up the composite operator, new ultraviolet divergences~\footnote{In the picture in which $\Lambda$ represents the smoothing scale at which the peak is defined, $\Lambda$ never truly diverges in CDM cosmologies since it is a physical scale. However, in the formalism we introduce here we wish to appeal to the reader's intuition of the QFT renormalisation scheme. Within this context, our terminology is similar to \cite{Assassi:2014fva}.} appear and depend on the ultraviolet cut-off $\Lambda\sim\mathcal{O}(1/R)$. 
The appearance of these divergences  requires the introduction of new counterterms. 
The physical reason why the term~\eqref{composite operator} changes in the equation for halos is that the force felt by peaks will be \textit{statistically} biased with respect to the DM case. The technical way by which we can reconstruct the form of $\vec \nabla\Phi\cdot\partial_{\vec v} f $ for the halo Boltzmann equation is through renormalization. 

Renormalized operators can be defined and generically expressed as linear combinations of all the bare operators of equal or lower canonical dimensionality (see, for instance, \cite{collins})
\begin{equation}}   \def\ee   {\end{equation}
\label{renorm definition [O]}
\left[{\cal O}_i(\vr,t)\right]=\sum_{j} Z_{ij}(\Lambda)\,{\cal O}_j(\vr,t).
\ee
Therefore, any composite operator can be expressed as a sum of operators allowed by the symmetries of the problem at hand.
The renormalized operator $[\vec \nabla\Phi \cdot \partial_{\vec v} f ]$ being a scalar field, it can mix with all the operators allowed by the extended Galilean symmetry \cite{con} and must remain a function of the velocity ${\vec v}$ and preserve the derivative $\partial_{v_i}$. 
We can classify the operators $\mathcal O_i$ that mix with $\vec \nabla\Phi \cdot \partial_{\vec v} f $ according to their number of derivatives in ${\vec v}$:
\begin{equation}
\label{composite op renorm}
\left[\vec \nabla\Phi \cdot \partial_{\vec v} f  \right]= \vec \nabla\Phi \cdot \partial_{\vec v} f   +
	\underset{\text{leading order in }\partial_{\vec v}}{\underbrace{Z_{01}\,\vec \nabla \delta\cdot \partial_{\vec v} f  + \dots}} +
	\underset{\text{higher order in }\partial_{{\vec v}}}{\underbrace{\tilde{Z}\vec \nabla\Phi \cdot \partial_{\vec v}^3 f + \dots}} \,.
\end{equation}

Once we have accomplished this classification, we select the leading operators according to the number of derivatives (we are considering large scales) and fields contained in them 
(adding additional fields brings to a result including more power spectra). 
The fields we are dealing with are $\Phi$ (which can only appear with derivatives, since the gravitational potential is not observable and  does not directly affect the dynamics) and the density contrast $\delta=\nabla^2\Phi/\alpha $, where $\alpha=3 H^2 \Omega_m/2$. We can then organise the $\mathcal O_i$ according to the number of extra couples of derivatives (in order to build a vector, we must add two derivatives to $\vec \nabla\Phi$ and contract two indices) and extra fields appearing in them. The result is shown in fig.~\ref{fig: operators}, where moving to the right means adding two derivatives, and moving downwards means adding a field. The appearance of derivative operators makes clear that the statistical bias manifests itself only through $k$-dependent corrections and, therefore, vanishes in the limit of small momenta.
\begin{figure}[h!]
\centering
\begin{tikzpicture}[node distance=1cm, auto,]
\node[operator] (o0) {$\mathcal O_0=\vec \nabla\Phi \cdot\partial_{\vec v} f $};
\node[operator, right=1.5cm of o0] (o1) {$\mathcal O_1=\vec \nabla \delta \cdot\partial_{\vec v} f $};
\node[operator, below=1.cm of o0] (o2) {$\mathcal O_2=\delta\vnabphi \cdot\partial_{\vec v} f $\\
$\mathcal O_3=(\vnabphi \hspace{-2pt} \cdot \hspace{-2pt} \vec\nabla)\vnabphi \cdot\partial_{\vec v} f $};
\node[operator,right=1.5cm of o1] (r) {$\cdots$};
\node[operator,right=1.5cm of o2] (c) {$\cdots$};
\node[operator,below=1.cm of o2] (l) {$\cdots$};
\node[right=-1mm of o0,align=center]{{\footnotesize add 2 der.}\\\vspace*{1em}};
\node[below=0.25cm of o0]{{\footnotesize add a field}\hspace*{5em}};
\node[] {} (o0.east) edge[freccia,bend left=0](o1.west);
\node[] {} (o0.east) edge[freccia,bend left=0](o1.west);
\node[] {} (o1.east) edge[freccia,bend left=0](r.west);
\node[] {} (o2.east) edge[freccia,bend left=0](c.west);
\node[] {} (o0.south) edge[freccia,bend right=0](o2.north);
\node[] {} (o1.south) edge[freccia,bend right=0](c.north);
\node[] {} (o2.south) edge[freccia,bend right=0](l.north);
\end{tikzpicture}
\caption{Operators allowed by rotational invariance that can mix with $\vec \nabla\Phi\cdot\partial_{\vec v} f $, classified according to their number of fields and derivatives. 
}
\label{fig: operators}
\end{figure}

The renormalization procedure requires a  suitable prescription. For our purpose, it is convenient to adopt the following renormalization condition: the correlators of the renormalized quantities should reproduce the correlators at the peak,
\begin{equation}}   \def\ee   {\end{equation}
\label{renorm prescription}
\Big<\ren{\mathcal O_1(\vr,t)}\cdot\ren{\mathcal O_2(\vr,t)}\cdots \ren{\mathcal O_n(\vr,t)}\Big>
=\Big<{\cal O}_1(\vr,t)\cdot{\cal O}_2(\vr,t)\cdots{\cal O}_n(\vr,t)\Big>_{\rm pk}.
\ee
Here and henceforth, the notation $\langle \cdots\rangle_\text{pk}$ will designate conditional averages at the peak position.
In our problem, it is precisely the  renormalization of the operator $\vec \nabla\Phi(\vr,t)\cdot \partial_{\vec v} f (\vr,\vv,t)$ that leads to differences in the clustering properties of halos and DM, as we will now see. Notice that, since we want to follow the evolution over cosmic time of the so-called proto-halos until their virialization, the renormalization condition has to be imposed at all times.


\subsection{Renormalization  at the linear order}
\hskip 0.8cm In this subsection, we study the renormalization of the composite operator $\vec \nabla\Phi\cdot \partial_{\vec v} f (\vr,\vv,t)$ at the linear order in perturbation theory.

We begin with the computation of the average of $\vec \nabla\Phi(\vr,t) \cdot \partial_{\vec v} f (\vr,\vv,t)$ at the peak position:
\begin{equation}}   \def\ee   {\end{equation}
\label{average composite op}
\Big<\vec \nabla\Phi(\vr,t) \cdot \partial_{\vec v} f (\vr,\vv,t) \Big>_{\rm pk}=\partial_{\vec v} f _{\rm pk}(\vr,\vv,t) \cdot \Big<\vec \nabla\Phi(\vr,t) \Big>_{\rm pk}\;.
\ee
This amounts to computing the correlators using the appropriate conditional probabilities.

Let us first operate at the linear level in perturbation theory such that, at early times, the DM  density $\delta_{\rm }$ and  the gravitational force $\vec\nabla\Phi$  are Gaussian variables. We can then apply the theorem stated, for instance,  in Refs.~\cite{bbks,bema}, which ensures that the conditional probability of zero-mean Gaussian variables $X_A$ and $X_B$ is itself a Gaussian variable with mean
\begin{equation}}   \def\ee   {\end{equation}
\Big<X_B\Big|X_A\Big> \equiv \frac{\Big<X_B \otimes X_A\Big>}{\Big<X_A \otimes X_A \Big>}X_A
\ee
and covariance matrix
\begin{equation}}   \def\ee   {\end{equation}
\label{covariance constrained}
{\rm C}\Big(X_B,X_B\Big) \equiv \Big<X_B \otimes X_B \Big> - \frac{\Big<X_B \otimes X_A\Big>}{\Big<X_A \otimes X_A \Big>}\Big<X_A \otimes X_B \Big>.
\ee
In our case, we identify $X_A$ with $\vec\nabla\delta_R$ and $X_B$ with $\vec\nabla\Phi_R$. As explained above, the gravitational force acting on the halo centers is smoothed to reflect the finite extent of the halos.
Since $\langle \vec \nabla\Phi\otimes \delta \rangle=\vec 0$ because of rotational invariance, we can simply compute $\langle \vec \nabla\Phi \rangle_\text{pk}$ as $\langle\vec \nabla\Phi_R|\vec \nabla \delta_R=\vec 0\rangle$.The mean shift of $\vec \nabla\Phi$ at the peak position is given by
\begin{equation}
\label{average grad phi}
\begin{aligned}
\Big<\vec \nabla\Phi\Big>_{\rm pk}&=
	\Big<\vec \nabla\Phi_R\Big |\vec \nabla \delta_R \Big> = \frac{\Big<\vec{\nabla}\Phi_R\cdot\vec{\nabla}\delta_R\Big>}{\Big<\left(\vec \nabla \delta_R \right)^2\Big>}\,\left.\vec{\nabla}\delta_R \right|_{\rm pk}
=-\alpha \frac{\Big<\delta^2_R \Big>}{\Big<\left(\vec{\nabla}\delta_R \right)^2\Big>}\,\left.\vec{\nabla}\delta_R \right|_{\rm pk} \\ 
& = -\alpha\,\frac{\sigma_0^2(\Lambda)}{\sigma_1^2(\Lambda)}\,\left.\vec{\nabla}\delta_R \right|_{\rm pk},
\end{aligned}
\end{equation}
where we have highlighted the dependence of the spectral moments $\sigma_j$ on the UV cut-off $\Lambda$. Since at the peak $\vec{\nabla}\delta_R$  is vanishing, no renormalization is needed.
We will hereafter drop the subscript $R$ whenever $\vec \nabla\Phi$ and $\delta$ appear because we will always refer to peaks, hence the smoothing is understood.
Notice also that, while the renormalisation condition is imposed at all times, the scale $\Lambda$ is fixed at initial time and does not evolve. In reality, $\Lambda$, which is related to the smoothing scale $R$ that define the Lagrangian extent of a halo, likely evolves with time (see e.g. \cite{Baldauf:2013hka}). For simplicity however, we will ignore this complication in the present work.
 
We wish now to compute the average of $\vec{\nabla}\Phi(\vr,t) \cdot \partial_{\vec v} f (\vr,\vv,t)$ with itself at the peak position. On applying Wick's theorem, the only non-vanishing term that remains is
\begin{multline}
\Big<\left(\vec \nabla\Phi(\vr,t) \cdot \partial_{\vec v} f (\vr,\vv,t)\right)^2\Big>_{\rm pk}=
\partial_{\vec v} f _{\rm pk}(\vr,\vv,t) \cdot \partial_{\vec v} f _{\rm pk}(\vr,\vv,t)\Big <\vec{\nabla}\Phi(\vr,t) \cdot\vec{\nabla}\Phi(\vr,t)\Big>_{\rm pk},
\end{multline}
where we used the fact that $\Big <\nabla_i\Phi(\vr,t) \cdot\nabla_j\Phi(\vr,t)\Big> \propto \delta_{ij}$.
This requires the knowledge of the covariance matrix for $\vec \nabla\Phi$ given the constraint at the peak. At the linear order the unconstrained correlators needed to evaluate eq. ~\eqref{covariance constrained} read 
\begin{equation}
\begin{aligned}
\Big<\vec{\nabla}\Phi\otimes\vec{\nabla}\Phi\Big>&=\frac{\alpha^2}{3}\sigma^2_{-1}(\Lambda)\,\mathbb{1}_{3\times 3},\\
\Big<\vec{\nabla}\Phi\otimes\vec{\nabla}\delta_{\rm }\Big>&=-\frac{\alpha}{3}\sigma^2_{0}(\Lambda)\,\mathbb{1}_{3\times 3},\\
\Big<\vec{\nabla}\delta_{\rm }\otimes\vec{\nabla}\Phi\Big>&=-\frac{\alpha}{3}\sigma^2_{0}(\Lambda)\,\mathbb{1}_{3\times 3},\\
\Big<\vec{\nabla}\delta_{\rm }\otimes\vec{\nabla}\delta_{\rm }\Big>&=\frac{1}{3}\sigma^2_{1}(\Lambda)\, \mathbb{1}_{3\times 3}.
\end{aligned}
\end{equation}

The covariance matrix for $\vec \nabla\Phi$ given the peak constraint is then equal to
\begin{equation}}   \def\ee   {\end{equation}
\label{cov}
{\rm C}\left(\vec{\nabla}\Phi,\vec{\nabla}\Phi\right)_\text{pk}=\frac{\alpha^2}{3}\left(\sigma^2_{-1}(\Lambda)-\frac{\sigma^4_{0}(\Lambda)}{\sigma^2_{1}(\Lambda)}\right)\,\mathbb{1}_{3\times 3},
\ee
and
\begin{equation}}   \def\ee   {\end{equation}
\label{cov1}
\Big <\vec{\nabla}\Phi(\vr,t) \cdot \vec{\nabla}\Phi(\vr,t)\Big>_{\rm pk}={\rm Tr}\, C\left(\vec{\nabla}\Phi,\vec{\nabla}\Phi\right)=\alpha^2\left(\sigma^2_{-1}(\Lambda)-\frac{\sigma^4_{0}(\Lambda)}{\sigma^2_{1}(\Lambda)}\right).
\ee
We recall now the prescription~\eqref{renorm prescription}: at the lowest order,  $[\vec \nabla\Phi \cdot \partial_{\vec v} f ]$ contains only the operators $\mathcal O_0$ and $\mathcal O_1$ (see fig.~\ref{fig: operators}). Hence, we learn that the renormalized force felt by the halos at the linear order reads
\begin{equation}}   \def\ee   {\end{equation}
\label{bb}
\left[\vec{\nabla}\Phi(\vr,t)\cdot \partial_{\vec v} f (\vr,\vv,t)\right]^{\rm first}=\left( \vec{\nabla}\Phi(\vr,t)+\alpha\frac{\sigma_0^2(\Lambda)}{\sigma_1^2(\Lambda)}\vec{\nabla}\delta_{\rm }(\vr,t)\right) \cdot \partial_{\vec v} f (\vr,\vv,t).
\ee
Although the mean shift of the gravitational force vanishes at the peak and thus there is no extra force at the peak-by-peak level, the gravitational force at the peak receives a correction when the statistical ensemble average is taken, as exemplified by eq.~\eqref{cov1}.
The extra effective force felt by the halos is purely of statistical origin and, therefore, does not violate the Equivalence Principle. 
Furthermore, the time-dependence of $Z_{01}$ is that of $\alpha \sigma_0^2/\sigma_1^2$, that is, the amplitude of the resulting $k^2$-correction in Fourier space {\it does not depend on time}. Therefore, this gravity bias does not decay with time, in agreement with \cite{dcss,v1}.
Note that the same result eq.~\eqref{bb} was found in Refs. \cite{bao,v0,v1,v2}, but our Boltzmann approach combined with a renormalization procedure puts it into a different perspective



\subsection{Renormalization at higher order with the path integral technique}
\label{sec:renorm ho}

\hspace{0.8cm} In this subsection, we will compute the renormalized operator $[\vec \nabla\Phi(\vr,t)\cdot\partial_{\vec v} f (\vr,\vv,t)]$ at the second order in perturbation theory. We will compute the correlators beyond the Gaussian approximation with the path integral technique (details can be found in appendix~\ref{sec:appendix A}). We will restrict ourselves to the four operators $\mathcal O_0,\,\mathcal O_1,\,\mathcal O_2,\,\mathcal O_3$ defined in fig.~\ref{fig: operators}.

Let us give a closer look to the renormalization prescription \eqref{renorm prescription}. We begin by observing that the expectation values of single operator $\langle\mathcal O_i\rangle$ are not helpful, because they are all proportional to expectation values of vector fields and thus vanish, as in~\eqref{average composite op} and \eqref{average grad phi}. Therefore, in order to determine the coefficients $Z_{ij}$ of \eqref{renorm definition [O]}, we write down the system given by the expectation values of products of two operators:
\begin{equation}
\label{renorm system}
\sum_{l,m} Z_{il} Z_{jm}\corr{l}{m} =\corrpk{i}{j}\,.
\end{equation}

We must now specify how we define the leading and subleading orders of our computation. On the left-hand side of eq.~\eqref{renorm system}, operating at second order the correlation functions $\corr{i}{j}$ contain at most two power spectra. Therefore, at leading order we compute the correlators by including only their Gaussian components and obtaining one power spectrum in the result; the next-to-leading order includes up to two power spectra. 
On the right-hand side, the leading order of the computation is obtained by considering again only the Gaussian component of the fields. 
When computing the correlators $\corrpk{i}{j}$, the linear order is equivalent to consider only the leading order in $\nu\sigma_0$.

The results for the correlators without the peak constraint are\footnote{see Appendix~\ref{sec:appendix B} for details about the calculation of these correlators, and for the definition of the finite quantities $\mathcal K_1$, $\mathcal K_2$)}
\begin{eqnarray}
  & \text{(leading order)} & \text{(next-to-leading order)}\nonumber\\
\corr{0}{0} = &\displaystyle \frac{\alpha^2}{3}\sigma^2_{-1}\delta_{ij}\\
\corr{0}{1} = & \displaystyle -\frac{\alpha}{3}\sigma^2_{0}\delta_{ij}\\
\corr{0}{2} = & & \frac{\alpha^2}{3}\mathcal K_1 \delta_{ij}\\
\corr{0}{3} = & & -\frac{\alpha^3}{6} \mathcal K_1 \delta_{ij} \\
\corr{1}{1} = & \displaystyle \frac{1}{3}\sigma^2_{1}\delta_{ij}\\
\corr{1}{2} = & & -\frac{\alpha}{3}\frac{17}{7}\sigma^4_{0}\delta_{ij}\\
\corr{1}{3} = & & - \frac{\alpha^2}{3} \mathcal K_2\, \delta_{ij} \\
\end{eqnarray}
where we are dropping the factors of $\partial_{\vec v} f$ to simplify the notation. The correlators calculated with the peak constraint read 
\begin{eqnarray}
\corrpk{0}{0} = &
  \displaystyle  \frac{\alpha^2}{3}\delta_{ij}\biggl[\left(\sigma^2_{-1} - \frac{\sigma^4_0}{\sigma^2_1}\right) & +
	(\nu\sigma_0)\left(\mathcal K_1\frac{1}{\sigma_0^2} -\frac{20}{21} \frac{\sigma_0^4}{\sigma_1^2}\right) \biggr] \\
\corrpk{0}{2}  = & &
  \frac{\alpha^2}{3} (\nu\sigma_0)\delta_{ij}	\left(\sigma^2_{-1} - \frac{\sigma^4_0}{\sigma^2_1}\right) \\
\corrpk{0}{3} = && 
   \frac{\alpha^3}{9} (\nu\sigma_0)\delta_{ij} \left(\sigma^2_{-1} - \frac{\sigma^4_0}{\sigma^2_1}\right)
\end{eqnarray}
We also notice that our renormalization prescription \eqref{renorm prescription} enforces $\corrpk{1}{i}=0\,\, \forall i$. This corresponds to a subsystem of \eqref{renorm system} (when $i=1$) whose right-hand side is 0, the solution of which is $Z_{1i}=0$ or $[\vec \nabla \delta\cdot\partial_{\vec v} f ]=0$. 
This is a consequence of our renormalization prescription: once we impose that  at the position ${\vec r}$ there is a peak defined by $\vec \nabla \delta({\vec r})=\vec 0$, the operator $\vec \nabla \delta$ after renormalization cannot be anything but zero.

We are now free to  impose the condition $\ren{\mathcal O_i}=\mathcal O_i +\dots$, i.e.  the diagonal entries satisfy the condition $Z_{ii}=1$ for $i\neq 1$.
Let us first discuss the system for the renormalization of the three operators $\{ \mathcal O_0,\mathcal O_1,\mathcal O_2\}$. We begin by inspecting the equation $\corrpk{0}{0} = Z_{0i}Z_{0j}\corr{i}{j}$, which contains all (and only) the coefficients $Z_{0i}$. We know that
\[ Z_{00}=1, \qquad 
Z_{01} = \alpha \frac{\sigma_0^2}{\sigma_1^2}+\mathcal O(\nu\sigma_0), \qquad  
Z_{02} = \mathcal O(\nu\sigma_0)
\,. \]
We should impose the equality $\corrpk{0}{0}=Z_{0i}Z_{0j}\corr{i}{j}$ separately between the leading order terms and the next-to-leading order ones. We  find that 
\[ Z_{01} = \alpha \frac{\sigma_0^2}{\sigma_1^2}+\mathcal O(\nu\sigma_0),\quad Z_{02}=(\nu\sigma_0)\left[\frac{1}{2 \sigma_0^2} + \frac{31 \sigma_0^4}{6 \left(-17\sigma_0^6 + 7 \mathcal K_1 \sigma_1^2\right)}\right]+\mathcal O\left((\nu\sigma_0)^2\right)\,.\]
We now look at the equations $\corrpk{0}{2}=Z_{0i}Z_{2j}\corr{i}{j}$ and $\corrpk{2}{2}=Z_{2i}Z_{2j}\corr{i}{j}$, which include the unknowns $Z_{20}$ and $Z_{21}$. These two equations, once we select the terms consistently with their order (for example, the connected contribution of $\corr{2}{2}=\langle \delta\vec \nabla\Phi \delta\vec \nabla\Phi  \rangle$ contains 3 power spectra and should not be put together with second order quantities) we get
\begin{align*}
Z_{20} =&\frac{\mathcal{K}_1 \sigma _1^2-\tfrac{17}{7} \sigma _0^6}{(\sigma _0^4-\sigma _{-1}^2 \sigma _1^2 )} + \underset{=\delta_\text{pk} Z_{00}}{\underbrace{\left(\nu\sigma_0\right)}} \,, \\
Z_{21} =&\alpha \left[ \sigma_0^2 \frac{\mathcal{K}_1 -\tfrac{17}{7} \sigma _0^2\sigma_{-1}^2}{(\sigma _0^4-\sigma _{-1}^2 \sigma _1^2 )} 
\pm \frac{1}{\sigma_1^2}\sqrt{-\frac{\left(\mathcal{K}_1 \sigma_1^2-\tfrac{17}{7} \sigma _0^6\right)^2}{(\sigma _0^4-\sigma _{-1}^2 \sigma _1^2 )}
+\sigma_0^2 \left[\left(\tfrac{17}{7}\right)^2 \sigma_0^6- \sigma_{-1}^2 \sigma_1^2\right]}
\right] +
\underset{=\delta_\text{pk} Z_{01}}{\underbrace{\alpha \frac{\sigma_0^2}{\sigma_1^2}\left(\nu\sigma_0\right)}} \,.
\end{align*}
This concludes the determination of the renormalization coefficients at the next-to-leading order for the set of operators $\{ \mathcal O_0,  \mathcal O_1,  \mathcal O_2\}$, keeping up to two power spectra in the (connected) correlators. 
We could now include the operator $\mathcal O_3$ in the discussion, since it contains the same number of fields $\Phi$ and derivatives acting on them as $\mathcal O_2$. 
But in this way we should compute correlation functions as $\corrpk{3}{3}$, whose connected part contains three power spectra, and thus goes beyond the perturbative order we are considering. 
The difference with respect to the operator $\mathcal O_2$ is that the correlator $\corrpk{2}{2}$, with the constraint $\delta=\nu\sigma_0$, gives $(\nu\sigma_0)^2$ times a correlator with one power spectrum. 
Now, $\sigma_0^2$ is of the same order of a power spectrum in a perturbative expansion (see eq. ~\ref{sigma def}), thus $\corrpk{2}{2}$ is of a lower order with respect to $\corrpk{3}{3}$. Therefore, at the perturbative order we are considering, $[\vec \nabla\Phi\cdot \partial_{\vec v} f ]$ does not mix with $\bigl((\vnabphi \hspace{-2pt} \cdot \hspace{-2pt} \vec\nabla)\vnabphi\bigr)\cdot \partial_{\vec v} f $.

We have therefore found that the next-to-leading order renormalization correction reads

\begin{eqnarray}
\label{bbNLO}
\left[\vec{\nabla}\Phi(\vr,t)\cdot \partial_{\vec v} f (\vr,\vv,t)\right]^{\rm NLO}&=&\left( \nu\sigma_0  (\Lambda)\left[\frac{1}{2 \sigma_0^2(\Lambda)}+ \right.\right.\nonumber\\
&+& \left.\left.\frac{31 \sigma_0^4(\Lambda)}{6 \left[-17\sigma_0^6(\Lambda) + 7 \mathcal K_1(\Lambda) \sigma_1^2(\Lambda)\right]}\right] \delta \vec \nabla\Phi\right) \cdot \partial_{\vec v} f (\vr,\vv,t).\nonumber\\
\end{eqnarray}

\vspace{3mm}

Let us summarize the main results of this section, as they are the central point of this work: we have renormalized the composite operator $[\vec \nabla\Phi\cdot \partial_{\vec v} f ]$ imposing the peak constraint. Already at first order, we have found that while the mean of $\vec \nabla\Phi$ is zero as in the unconstrained case, the variance has an additional term. The composite operator therefore renormalized in such a way as to generate the correct variance and account for this statistical effect, in \eqref{bb}. We then compute the next-to-leading order effect and find the corresponding renormalization, \eqref{bbNLO}.


\subsection{Non-renormalization of the DM Boltzmann equation}
\label{sec:renorm ho}

\hspace{0.8cm} 
It is crucial to notice that for the smooth component of DM (that is without peak constraint) no renormalization is needed and the corresponding Boltzmann equation is not altered. This is consistent with the fact that we have used standard perturbation theory results for DM to evaluate the statistical correlators at the peak locations.
Consistency can be explicitly checked upon noticing that, for DM, probabilities are not conditional. For  instance, we get at linear order
\begin{eqnarray}
\Big<\vec{\nabla}\Phi(\vr,t)\, f(\vec{r},\vec{v},t)\Big>&=& f(\vec{r},\vec{v},t)\,
\Big<\vec{\nabla}\Phi(\vr,t) \Big>=0, \nonumber\\
\Big<\vec{\nabla}\Phi(\vr,t) \, f(\vec{r},\vec{v}_1,t)\cdot \vec{\nabla}\Phi(\vr,t) \, f(\vec{r},\vec{v}_2,t)\Big>&=&\Big<f(\vec{r},\vec{v}_1,t) \, f(\vec{r},\vec{v}_2,t)\Big>\Big <\vec{\nabla}\Phi(\vr,t) \cdot \vec{\nabla}\Phi(\vr,t)\Big>\nonumber\\
&=&\frac{\alpha^2}{3}\sigma^2_{-1} \Big<f(\vec{r},\vec{v}_1,t)\,  f(\vec{r},\vec{v}_2,t)\Big>
\end{eqnarray}
and one simply finds 
\begin{equation}}   \def\ee   {\end{equation}
\left[\vec{\nabla}\Phi(\vr,t)\right]^{\rm linear}_{\rm }=\vec{\nabla}\Phi(\vr,t).
\ee
This remains true also at higher orders. 


\section{The halo Boltzmann equation: the right-hand side}
\label{sec:boltzmann eq RHS}
We will now deal with the right-hand side of the Boltzmann equation, see Eqs.~(\ref{force}) and (\ref{irr}). 
It contains a force term that depends on the irreducible two-particle correlation (which, in turns, depends on the irreducible three-particle correlation etc.).
It is the correlation  encoded by $f_{\rm 2c}$ that causes the non-conservation of the one-particle phase-space distribution function in  configuration space. The physical interpretation of this effect is that the gravitational interactions between particles  (which induce long-range correlations) lead to clustering and, therefore, non-conservation of the one-particle phase-space distribution. 

To calculate this extra  force, we follow Ref.~\cite{bema}. The procedure consists in relating  phase-space densities to probability distributions for the mass and velocity field, which are naturally specified by models of cosmological structure formation. 
As explained in sec.~\ref{sec:boltzmann equation}, we can express the phase space density $f({\vec r},{\vec v},t)$  in terms of the probability distribution $P(\rho,{\vec v},t)$. 
This approach can be extended to the irreducible two-particle correlation. Inspecting   eqs.~\eqref{force} and \eqref{irr}, we observe that
$f_{\rm 2c}({\vec r},{\vec v},\vec{r}^{\, \prime},\vec{v}^{\, \prime},t)$ is integrated over $\d^3 v'$ and we needs only the velocity at one point (the superscript $'$ indicates that quantities are evaluated at $\vec{r}^{\, \prime}$)
\begin{eqnarray}
\label{as}
f_{\rm 2c}({\vec r},{\vec v},\vec{r}^{\, \prime},t)&=&\int\d^3 v'\,f_{\rm 2c}({\vec r},{\vec v},\vec{r}^{\, \prime},\vec{v}^{\, \prime},t) \nonumber \\
&=&\int_0^\infty\,\d\rho\,\rho\, 
\int_0^\infty\,\d\rho'\,\rho'\,\left[ P(\rho,{\vec v},\rho',t)-P(\rho,{\vec v},t)P(\rho',t)\right]\nonumber\\
&=& \int_0^\infty\,\d\rho\,\rho\, 
\int_0^\infty\,\d\rho'\,\rho'\,\left[ P(\rho,\rho'|{\vec v},t)-P(\rho|{\vec v},t)P(\rho',t)\right]P({\vec v},t)\nonumber\\
&=&\left[
\langle\rho\rho'|{\vec v}\rangle-\langle\rho|{\vec v}\rangle\langle\rho'\rangle
\right]P(\vec{v}^{\, \prime},t)\;,
\end{eqnarray}
where, in the last equality, the averages are equal-time correlators.
Even though one derives these equations assuming a
single-valued velocity field in each realization, they are also
valid when there is a distribution of velocities at each point. We
simply interpret $\rho$  as the total mass density while ${\vec v}$  is a single bulk 
velocity. Therefore, these equations are fully general \cite{bema}.

Let us specialize to the case of halos. The phase-space distribution (\ref{f integral rho}) can be written using Bayes' theorem as 
\begin{equation}}   \def\ee   {\end{equation}
\label{fpk}
f_{\rm pk}({\vec r},{\vec v},t)=\int_0^\infty\,\d\rho\,\rho\, P({\vec v}|\rho,t)\,P(\rho,t)=\rho_{\rm pk}\, P({\vec v},t | {\rm pk})=
				 \bar{\rho}(1+\delta_{\rm pk})\; J_{{\vec v}\vec{\psi}}^{-3}\,P(\vec{\psi},t | {\rm pk})
\ee
where, in the last passage, we have introduced the quantity
\begin{equation}}   \def\ee   {\end{equation}
\vec{\psi}(\vr,t)=-\frac{1}{\alpha}\vec{\nabla}\Phi(\vr,t)=-\int\frac{\d^3 r'}{4\pi}\, \frac{({\vec r}-\vec{r}^{\, \prime})}{|{\vec r}-\vec{r}^{\, \prime}|^3}\delta_{\rm }(\vec{r}^{\, \prime},t),
\ee
and $J_{{\vec v}\vec{\psi}}$ is the Jacobian describing the passage from the variable ${\vec v}$ to the variable $\vec{\psi}$.

In a similar fashion, we can compute the irreducible two-point correlation function
\begin{eqnarray}
\label{ad}
f^{\rm pk}_{\rm 2c}({\vec r},{\vec v},\vec{r}^{\, \prime},t)&=& \int_0^\infty\,\d\rho\,\rho\, 
\int_0^\infty\,\d\rho'\,\rho'\,\left[ P({\vec v},\rho'|\rho,t)P(\rho,t)-P({\vec v}|\rho,t)P(\rho,t)P(\rho',t)\right]\nonumber\\
&=&\rho_{\rm pk}\int_0^\infty\,\d\rho'\,\rho'\, \left[P({\vec v},\rho',t|{\rm pk})-P({\vec v},t|{\rm pk})P(\rho',t)
\right]\nonumber\\
&=&\rho_{\rm pk}\, P({\vec v},t|{\rm pk})\,\int_0^\infty\,\d\rho'\,\rho'\, \left[P(\rho',t|{\vec v},{\rm pk})-P(\rho',t)
\right]\nonumber\\
&=&\bar{\rho}^2(1+\delta_{\rm pk})\; J_{{\vec v}\vec{\psi}}^{-3}\,P(\vec{\psi},t | {\rm pk})\left[\Big<\delta(\vec{r}^{\, \prime},t)\Big|\vec{\psi}(\vr,t),{\rm pk}\Big>-\Big<\delta(\vec{r}^{\, \prime},t)\Big>\right]
\nonumber\\
&=&\bar\rho\; f_{\rm pk}({\vec r},{\vec v},t)\,\Big<\delta(\vec{r}^{\, \prime},t)\Big|\vec{\psi}(\vr,t),{\rm pk}\Big>\;.
\end{eqnarray}
In the last equality, we have substituted $\big<\delta(\vec{r}^{\, \prime},t)\big>=0$ as this average vanishes in the absence of peak constraint.
The corresponding force felt by the halos is therefore
\begin{equation}}   \def\ee   {\end{equation}
\label{forcepeak}
\vec F_{\rm pk} = G_{\rm N}\bar\rho f_{\rm pk}({\vec r},{\vec v},t)\,  \int \d^3 r'\,\Big<\delta(\vec{r}^{\, \prime},t) \Big| \vec{\psi}(\vr,t),{\rm pk} \Big>\cdot \frac{{\vec r} - \vec{r}^{\, \prime}}{|{\vec r} - \vec{r}^{\, \prime}|^3}.
\ee
The computation of the force (\ref{forcepeak}) is reported in Appendix \ref{sec:appendix C} and we provide only the final result up to second-order
\begin{equation}
\label{eq:force peak result}
\begin{aligned}
\vec F_{\rm pk} ({\vec r},{\vec v},t) =& -4\pi\,G_{\rm N}\, \bar\rho\,f_{\rm pk}({\vec r}, {\vec v}, t)\,\Biggl[\vec\psi(\vr,t) +\left(\vec\psi(\vr,t) -\frac{\sigma^2_0}{\sigma^2_1} \vec\nabla\delta(\vr,t)\right)\frac{\delta(\vr,t)}{\sigma^4_0-\sigma^2_1\sigma^2_{-1}}\frac{\sigma_1^2}{\sigma_0^2}\\
	&\times \left(\Big<\delta(\vr,t)\vec{\psi}(\vr,t)\cdot\vec{\psi}(\vr,t)\Big>-\int\d \ell\,\frac{\ell_i}{\lvert\ell\lvert}\,\Big<\delta(\vec{r}^{\, \prime},t)\delta(\vr,t)\psi_i(\vr,t)\Big> \right)\\
	&+\left(\vec\psi_1 -\frac{\sigma^2_{-1}}{\sigma^2_0} \vec\nabla\delta(\vr,t)\right)\frac{\delta(\vr,t)}{\sigma^4_0-\sigma^2_1\sigma^2_{-1}} \biggl(\int\d \ell\,\frac{\ell_i}{\lvert\ell\lvert}\,\Big<\delta(\vec{r}^{\, \prime},t)\delta(\vr,t)\partial_i\delta(\vr,t)\Big> \\
	&-\Big<\delta(\vr,t)\vec{\psi}(\vr,t)\cdot\vec{\nabla}\delta(\vr,t)\Big> \biggr)\Biggr],
\end{aligned}
\end{equation}
where $\vec\ell=(\vec{r}^{\, \prime}-{\vec r})$. Using the results for the cumulants reported in \eqref{mu233}, \eqref{mu234}, \eqref{mu123}, \eqref{mu124} we find (notice that the last two lines in \eqref{eq:force peak result} cancel each other)
\begin{multline}
\label{eq:force peak result final}
\vec F_{\rm pk} ({\vec r},{\vec v},t) = -4\pi\,G_{\rm N}\, \bar\rho\,f_{\rm pk}({\vec r}, {\vec v}, t)\,\Biggl[\vec\psi(\vr,t)
  +\left(\vec\psi(\vr,t) -\frac{\sigma^2_0}{\sigma^2_1} \vec\nabla\delta(\vr,t)\right)\\
  \times\frac{\delta(\vr,t)}{\sigma^4_0-\sigma^2_1\sigma^2_{-1}}\frac{\sigma_1^2}{\sigma_0^2} \left(\frac 32 \mathcal K_1 -\frac{10}{21}\sigma_0^2 \sigma_{-1}^2\right)\Biggr].
\end{multline}
Two comments are in order. Firstly, we can easily recover the DM case away from the peaks, either by computing directly the force from expression (\ref{ad}) with aid of $\langle\delta(\vec{r}^{\, \prime},t)\vec{\psi}(\vr,t)\rangle$, or by eliminating the peak constraint taking $\sigma_0^2,\sigma_1^2\gg 1$. In both cases, we obtain
\begin{equation}
\vec F_{\rm dm} ({\vec r},{\vec v},t) = -4\pi\,G_{\rm N}\, \bar\rho\,f_{\rm dm}({\vec r}, {\vec v}, t)\,\vec\psi(\vr,t).
\end{equation}
Secondly, neither the force on the peaks nor on the DM contributes to the momentum equation once the integration of the force over momenta is considered,
\begin{equation}}   \def\ee   {\end{equation}
\int\d^3 v\, {\vec v}\, \vec\nabla_{\vec v}\cdot\vec F_{\rm pk,dm} ({\vec r},{\vec v},t) =\vec 0,
\ee
in agreement with the Equivalence Principle. Technically, this cancellation follows from the fact that integrating out the velocities is equivalent to an integration over the configurations of the field $\vec{\psi}$. As the phase-space distribution is proportional to $P(\vec{\psi},t)$, we obtain
\begin{equation}}   \def\ee   {\end{equation}
\int\d^3\psi\,P(\vec{\psi},t)\, \vec{\psi}=\big<\vec{\psi}\big>=\left\{
\begin{aligned}
& -\alpha\,\frac{\sigma_0^2(\Lambda)}{\sigma_1^2(\Lambda)}\,\left.\vec{\nabla}\delta_{\rm }\right|_{\rm pk}=\vec 0 && (\text{for peaks})\,, \\
& \ \vec 0 			& & 			 (\text{for DM})\,,
\end{aligned}\right.
\ee
where, for peaks, we made use of eq.~(\ref{average grad phi}). 
This equation seems to suggest that one would break the equivalence principle if the constraint were different from $\vec{\nabla}\delta=\vec 0$, as may be the case for other tracers of the large scale structure. However, isotropy and homogeneity would enforce a constraint on the modulus of $\vec{\nabla}\delta$ only, e.g. $|\vec{\nabla}\delta|\leq c$ for some constant $c>0$. Therefore, the integral over $\vec{\psi}$ would also vanish in those cases, and the Equivalence Principle would still be satisfied, as it should be.


\section{Discussion and conclusion}
\label{sec:conclusions}
The fine understanding of the large scale structure of our Universe is a central theme in cosmology. 
Our knowledge in the field is nowadays (and will be even more in the future) dominated by data.  
It is, therefore, essential to understand in detail the statistical properties of the observed galaxies. 
Achieving this goal automatically calls for a detailed investigation of the statistical properties of the 
DM halos, where galaxies are expected to reside.

In this paper we have derived an effective Boltzmann equation that describes the dynamics of the halo 
mean-field phase space. 
This equation is necessarily  different from the Boltzmann equation of DM as halos are not statistically 
unbiased relatively to the coarse-grained DM density field. 
Through a renormalization procedure of composite operators, we have obtained up to second-order in 
perturbation theory the statistically biased gravitational force felt by halos. 
We conclude that, at this order, the Boltzmann Equation for halos reads
\begin{multline}\label{final}
\frac{\partial f_{\rm pk}}{\partial t}+\vec{v}\cdot \frac{\partial f_{\rm pk}}{\partial{\vec r}} 
-\Biggl(\vec \nabla\Phi_R(\vr,t)+\alpha\frac{\sigma_0^2(R)}{\sigma_1^2(R)}\vec{\nabla}\delta_R(\vr,t) + \\
+ \nu\sigma_0(R)\left[\frac{1}{2 \sigma_0^2(R)} + \frac{31 \sigma_0^4(R)}{6 \left[-17\sigma_0^6(R) + 7 \mathcal K_1(R) \sigma_1^2(R)\right]}\right] \delta_R(\vr,t) \vec \nabla\Phi_R(\vr,t)\Biggr)\cdot\frac{\partial f_{\rm pk}}{\partial{\vec v}}=0 \;,
\end{multline}
where the extra contributions to $\vec \nabla\Phi_R(\vr,t)$ are purely of statistical origin: they do not 
arise because halos feel a different gravitational force than dark matter, as this would violate the equivalence principle, but rather as a result of imposing the peak constraint on their statistics\footnote{ An alternative, perhaps more intuitive, interpretation of this result is the following: instead of modifying the Boltzmann equation for halos with an effective force field, one may interpret the results of equations \eqref{bb} and \eqref{bbNLO} as the corrections to the probability distribution function (PDF) of the force acting on halos. In this picture, the Boltzmann equation for halos \eqref{final} would look exactly as the one of DM, equation \eqref{eq:dmfinal}, but with a conditional PDF of the form of the terms in bracket in equation \eqref{final}. }.
This should be contrasted to the DM Boltzmann equation
\begin{equation}}   \def\ee   {\end{equation}\label{eq:dmfinal}
\frac{\partial f_{\rm }}{\partial t}+\vec{v}\cdot \frac{\partial f_{\rm }}{\partial{\vec r}}
-\vec{\nabla}\Phi(\vr,t)\cdot\frac{\partial f_{\rm }}{\partial{\vec v}}=0 \;.
\ee
We emphasize again that $f_\text{pk}$ is the one-particle phase space density associated with the mean-field 
distribution of halos that corresponds to a particular $f$, that is, a particular realization of the large 
scale structure. 
The biased gravitational force experienced by the halo centers imprints a signature, among others, in the 
peculiar velocities of virialized halos. 
This generates the first-order statistical velocity bias which has been discussed in \cite{v0,dcss,v1} and 
measured in \cite{porciani,v1}.
However, tidal forces etc. will also be affected. Our approach can in principle be applied to work out these
corrections at any order in perturbation theory.

To give intuition about the magnitude of the statistical correction to the gravitational force in 
eq.~\eqref{final}, we plot in fig. \ref{xir} the amplitude of the second and third terms in parenthesis
\begin{eqnarray}
\xi_1(R) &=& \alpha\frac{\sigma_0^2(R)}{\sigma_1^2(R)},\\
\xi_2(R) &=& \nu\sigma_0(R)\left[\frac{1}{2 \sigma_0^2(R)} + \frac{31 \sigma_0^4(R)}{6 \left[-17\sigma_0^6(R) + 7 \mathcal K_1(R) \sigma_1^2(R)\right]}\right]
\end{eqnarray}
as a function of the Lagrangian scale $R = 1/\Lambda$ assuming the biasing is described by the peak 
constraint.
\begin{figure}[h!]
\centering
\includegraphics[width=0.8\textwidth]{vb.pdf}
\caption{$\xi_1(R)$ and $\alpha\xi_2(R)$ plotted here as a function of the Lagrangian scale $R = 1/ \Lambda$. Here we normalize $\xi_2(R)$ with a factor of $\alpha$ which relates $\delta$ with $\nabla \phi$ in \eqref{final}.}
\label{xir}
\end{figure}
The magnitude of the second-order statistical correction also rises with $R$, in agreement with the 
expectation that these effects should increase with halo mass.

Measurements from N-body simulations are challenging since the effect must vanish in the limit $k\to 0$. 
Furthermore, the sparsity of massive halos could bias the measurements when the purpose is to define 
volume-weighted statistics (see e.g. \cite{zzj2} for recent discussions). 
Note, however, that number-weighted statistics are most suited to extract these effects as they are the
quantities being observed.


There is some confusion in the literature regarding the nature of this effect. A physical bias will immediately 
arise as soon as there are at least two different fluids since, in this case, velocities do not need to share
the same value. This bias is effective on an object-by-object basis.
This is obviously the case of e.g. baryons and DM. Furthermore, this will also be the case of 
the halos and DM since, owing to the finite halo size, the force acting on the halo center-of-mass, $\vec \nabla\Phi_R$, 
is different from the force acting on DM particles, $\vec \nabla\Phi$. 
Such a coarse-graining procedure can introduce additional terms in the fluid equations, as was already pointed 
out in \cite{buchert}. This possibility was recently reconsidered in \cite{chan1}. 
If we ignore in Eq.(\ref{final}) all the statistical corrections but the term $\vec \nabla\Phi_R$, then our results 
are identical to his, except for the time constancy of our $R$. 
As shown in \cite{chan1}, a time-dependent filter $W({\vec r}/R,t)$ would yield extra terms in the moments of the 
Boltzmann equation owing to the time derivative $\partial f/\partial t$. 


However, \cite{chan1} did not include the statistical corrections, which are the focus of this work. In fact,
our results would hold even in the extreme case where no extra filtering is applied to the macrostate. This 
illustrates the fact that the statistical effects we are discussing are {\it not} caused by the finite extent 
of the halos.
Regarding the existence of these statistical effects, we note that, while a time-dependent window $W({\vec r}/R,t)$ 
mimicking the collapse of halos can induce $k$-dependent contributions at low redshift, corrections are small 
at high redshift \cite{chan1}. Therefore, this cannot be the explanation for the suppression of the proto-halo 
velocity power spectrum measured in \cite{porciani,v1}.


To conclude, note that halo number conservation implies that we follow the collapse of Lagrangian ``peak-patch'', 
as coined by \cite{bond/myers:1996}, into parent or isolated DM halos. Therefore, we do not explicitly track the 
merging history of these host halos. Although one could naively apply the extended Press-Schechter (EPS) approach
\cite{bond/cole/etal:1991} to determine the merger history, it would be very interesting to assess whether the 
statistical effects considered here also affect merger rates etc. However, it is not obvious how to do this within 
our current formulation, especially since the smoothing scale $R$ is linked to the parent halo mass and, thus, is 
not free unlike in EPS theory. Our assumption that the density peak retain their initial critical height (i.e. 
$\nu$ does not depend on time) is another limitation of our approach. A better treatment would consistently 
track the time evolution of both $R$ and $\nu$ such that the mass enclosed within $R$ always equates the host halo
mass. We leave all this for future work.

\acknowledgments
We acknowledge an anonymous referee for her/his careful reading and helpful suggestions.
The research of A.K. was implemented under the Aristeia Action of the Operational Programme Education and 
Lifelong Learning and is co-funded by the European Social Fund (ESF) and National Resources. 
A.R. is supported by the Swiss National Science Foundation (SNSF), project `The non-Gaussian Universe" 
(project number: 200021140236). 
M.B. and V.D. acknowledge support by the SNSF.


