\section{Introduction}
\noindent Murphy's Law is not a law in the formal sense; instead it is an interesting satirical statement: \emph{If something can go wrong then it will definitely go wrong} \cite{bloch2003murphy}. The origin of the statement is unknown (even Murphy himself!), yet it is not very hard to find it relevant in various situations, either in its rudimentary form or in the form of some derivative. Popular science often relates this statement to the Second Law of Thermodynamics, which states that the entropy of the universe will always increase with time \cite{handscombe2004entropy, cengel2007introduction}. The analogy between these two seemingly very dissimilar statements seems to crop up from the idea that left to themselves, things tend to get disorganized with time and Murphy comes into the picture, albeit in a modified form, stating that bad will eventually turn worse. In the realm of philosophy these statements may seem connected, however from a scientific viewpoint they are clearly unrelated, as (\emph{i}) one is a universal law, whereas the other, a rather ambiguous statement; (\emph{ii}) the second-law strictly talks about energy dissipation in physical processes and their irreversibility, whereas Murphy's statement (`law' is a misnomer in this case) is just a popular imagination, and does not even remotely relate to any of the physical or mathematical concepts to rely upon. Surprisingly, given these discrepancies, we still find it very interesting when the law is formulated as a mathematically equivalent statement. In this paper, we try to \emph{actually} relate Murphy's statement with the Second Law of Thermodynamics, and draw inspirations from energy and entropy. We check the correctness of the statement from a scientific standpoint, by first constructing a mathematically equivalent statement using concepts from mathematical logic and probability theory. Next, we understand the physical essence of the terms `rightness' and `wrongness' from the context of energy and time, using the Second Law and Principle of Least Action. Finally, we relate our findings to check the truthfulness of Murphy's statement. We observe that what initially seemed to be completely unrelated statements - the Second Law and Murphy's statement - \emph{indeed} have some interesting tales to tell.
\section{Methodology}
\subsection{Murphy and Logic}
Murphy's statement is more of a psychological statement having philosophical connotations of `rightness' and `wrongness'. A statement like this is not only hard to evaluate but even more challenging to prove (or disprove). The first step in any proof strategy lies in the rigorous formulation of the statement which is under scrutiny, and hence it is necessary that we try to formulate this rather anecdotal philosophy into a logical one. Doing so lets us set limits to the premise of the statement, as well as to its domain of application. However, in order to proceed we need to make some underlying assumptions. The idea of using propositional logic to build a mathematically coherent statement, even before attempting to prove or disprove it, is a necessity. Whereas, assigning an additional element of probability to the statement(s) is actually non-intuitive, because (\emph{i}) Murphy's speaks about occurrence of events and in mathematics, occurrence is always associated with \emph{chance}; (\emph{ii}) the presence of the intensifier, \emph{`definitely'}, in the linguistic construction of the statement itself, further strengthens our intuition of associating the statement(s) with finite values of probability. Logically speaking, Murphy's statement can be broadly divided into two distinct components: \textbf{1}. (\emph{If}) something can go wrong, and \textbf{2}. (\emph{then}) it will definitely go wrong. In propositional logic, this can be viewed as a causal chain of statements, $p \rightarrow q$, where $p$ is the first component, $q$ the second, and $p$ causes $q$. Such an argument is known as \emph{modus ponnens} \cite{tarski1993introduction} which asserts that the occurrence of $q$ is dependent upon the occurrence of $p$, or if $p$ happens, $q$ definitely happens. It can be represented mathematically as $\frac{p \rightarrow q, \enspace p}{\therefore, \enspace q}$. It was noted earlier that the above statements $p \rightarrow q$ have an inherent probability associated with them. The presence of an intensifier, `definitely' further affirms the occurrence of the latter half of the statement, irrespective of the magnitude of the chances of occurring of the former half. Mathematically, if the probability of occurrence of $p$ is \emph{finite} (remember, the probability that $\neg p$ occurs is also finite) then $q$ occurs \emph{almost surely}, \emph{i.e.}, with probability one (note that, even though we are interested in the happening of something which is not desired). What we earlier stated as a \emph{modus ponnens}, interestingly transforms into \emph{modus tollens} once we associate finite probabilities to the occurrence of the events (\emph{modus ponnens} is the law of affirming by affirming, therefore we affirm $p$ in order to justify $q$, and Murphy's statement affirms $q$ by associating finite chances of occurrence of $p$) \cite{tarski1993introduction}. According to the Law of Contrapositive: (\emph{modus tollens}), $\frac{p \rightarrow q, \enspace \neg q}{\therefore, \enspace \neg p}$, or if probability of occurrence of $p$ is \emph{finite} then $q$ occurs \emph{almost surely} and $q$ did not occur \emph{almost surely} enforces the fact that $p$ did not have a finite chance of occurring in the first place, or $\neg q \Rightarrow \neg p$ \cite{billingsley2008probability}. This is the form of Murphy's statement that we will adhere to throughout this paper. Also, in the following sections of this paper, we will analyze Murphy's statement using fundamental laws of physics, the Least Action Principle and the Second Law of Thermodynamics. It may seem puzzling at first as to how energy and entropy come into the picture, but a deeper inspection will reveal the intimate connection that Murphy's statement holds with respect to a system, process and entropy.
\subsection{Murphy, Thermodynamics and Action}
\noindent The connection between Murphy's statement and thermodynamics is the central theme of this paper. As proposed earlier, when we talk about Murphy's statement we refer to the probability of occurrence of any event in space and time. We thus make a transition from a system-theoretic perspective to a process-theoretic one, and once we do that we inevitably include the Second Law into the picture. Where, thermodynamics is a branch of science dealing with energy interactions with systems, energy conversion and energy degradation; it has also given birth to one of the most fundamental laws of nature, known as the Second Law or the Entropy Principle \cite{cengel2007introduction, clausius1854veranderte, planck2013treatise, lieb1999physics}. The Second Law deals with processes in nature, and acts as a standard test to justify the occurrence or nonoccurrence of any phenomenon. Every process (phenomenon) in nature can be visualized as a flow of energy. When this flow is spontaneous, \emph{i.e.}, the driving force behind the flow is the presence of gradients (like, energy, temperature, pressure, concentration, etc.), then such a flow is realized with an equivalent increase in entropy \cite{sharma2007natural, kaila2008natural, gronholm2007natural, pernu2012natural}. The flow of energy during any arbitrary process when coupled with the dimension of time gives rise to the notion of \emph{action}, formally represented as $\mathcal{A} = \int_{A \to B} (T-V)dt = \int_{A \to B} (\delta E)dt = \int_{A \to B} \mathcal{L}(p_{i},q_{i})dt$, where the pair, $(q_{i},p_i)$ is the generalized position-momentum pair, $t$ is the time, $\delta E$ is the change in energy along a specific path (in this case $A \to B$), $T$ and $V$ are the kinetic and potential energies respectively and $\mathcal{L} (p_i , q_i)$, the Lagrangian \cite{goldstein1957classical, de1751essais, euler1952methodus, lagrange1853mecanique}. There are several interpretations of the Action Principle, for example, the laws of physics when formulated in terms of the Action Principle identify any natural process as the one in which energy differences are leveled off in the least possible time (Maupertius' formulation) \cite{kaila2008natural, pernu2012natural}; among all existing possible paths, the path of least action is the one along which a natural process must proceed (Hamiltonian formulation) \cite{georgiev2002least, georgiev2012quantitative, georgiev2014mechanism, georgiev2014spontaneous, chatterjee2012action, chatterjee2013principle, chatterjee2011certain, atanu}. On the energy landscape, there is a unique trajectory that directs (energy) flows along the lines of steepest descent. A classical example of this principle can be observed in case of the Brachistochrone problem. The Principle of Least Action, thus dictates the existence of a unique path out of numerous possible trajectories (see figure 1). While some are thermodynamically feasible, others are not. From a thermodynamic perspective decrease in action is equivalent to an increase in entropy (we restrict ourselves to non self-organizing systems and processes), or along a trajectory, both the quantities can be represented along mutually orthogonal directions but with opposite signs, $\mathcal{L}=-T\mathcal{S}$ \cite{atanu, lucia2008probability, lucia2012maximum}. Murphy's statement vaguely mentions processes, but stresses more on their qualitative aspect of `rightness' and `wrongness'. It does not talk about the system, per se, undergoing the process, nor does it go into the details of the mechanism underlying the process; it basically concentrates only on the qualitative aspects of the outcomes. In order to test the correctness of the statement, it is absolutely necessary to look into the details of the process, as focusing on the finer details of the process will enable us to justify the outcomes.
\subsection{Defining a Probability Measure}
\noindent A deeper inspection is needed when we talk about processes and trajectories (pathways to achieve a particular outcome during a process). When we talk about Murphy's statement, we must invariably refer to processes and the trajectories (pathways) along which these processes occur (flow). In order to model a physical phenomenon into a thermodynamic one, we must explicitly distinguish between a system, surroundings and the set of processes that the system undergoes. In the algebraic formulation of variational thermodynamics, these terms are defined in a detailed, mathematically rigorous yet an abstract fashion \cite{atanu, lucia2008probability, lucia2012maximum, lucia1995mathematical, caratheodory1909untersuchungen, truesdell1984rational}. Roughly speaking, a system is a set of all those elements which are of current interest, and there exists a boundary that separates the system from its surroundings, where surroundings (relative to the system) refer to everything else in the universe except the system and the system boundary. A process is a change of state of a subset of system elements or the system itself in time. Mathematically, a process can be represented as a map that captures the changes incurred by the system or a subset of a system. We define a system undergoing a process by a tuple, ($\Omega , \Pi$), where the elements of $\Omega$ are the system elements/constituents/agents, and the elements of $\Pi$ are processes \cite{atanu}. We identify an event as a \emph{spontaneous} process in space-time if it proceeds along the path that reduces action and increases entropy \cite{atanu, lucia2008probability, lucia2012maximum}. The Principle of Least Action is that tool which helps us in finding that particular trajectory out of a set of uncountable possibilities along which an event must proceed. Using the standard notation, we define action for an arbitrary element in a system undergoing an arbitrary process, as $\mathcal{A}_{ij}=\int_{j}p_i dq_i$, where $i \in \mathcal{I}$, the set of all system elements in $\Omega$, and $j \in \mathcal{J}$, the set of all arbitrary processes, $\pi_j\in \Pi$ \cite{atanu}. For any arbitrary process, $\pi_j$ there exist numerous possibilities or trajectories, $\gamma_k$. We define the set $\Gamma^j$ to be the set of all such trajectories, $\gamma_k$ for the corresponding process, $\pi_j$ $\in \Pi$ and $\forall k \in \mathcal{K}$ ($\mathcal{I}$, $\mathcal{J}$, and $\mathcal{K}$ are index sets). It is now crucial to define what exactly is meant by `right' and `wrong' events as this is the central idea of Murphy's statement. We define two subsets of trajectories, $R$ and $W$, such that both $(R, W) \in \Gamma^j$. The trajectories which are thermodynamically feasible form the set, $R$, and the ones which are not, form the set, $W$. Therefore, the elements of $\Gamma^j$ are either labeled by $\gamma^{w}_k$ or $\gamma^{r}_k$ (for representation purpose only). The right trajectories are denoted by $\gamma^r_{k} \in R \subset \Gamma^j$, and the wrong ones by $\gamma^w_{k} \in \Gamma^j \setminus R \subset \Gamma^j = W$. Clearly, the subsets, $R$ and $W$ partition the set, $\Gamma^j$, so $R \cup W = \Gamma^j$ and $R \cap W = \phi$ (see figure 2). It is not very hard to see that the tuple, ($\Gamma^j , \mathcal{G}$) forms a $\sigma$-algebra\footnote{Let $\mathcal{G}$ be a collection of sets of $\Gamma^j$, then ($W$, $R$) $\in \mathcal{G}$, ($W^c$, $R^c$) $\in \mathcal{G}$, $\phi \in \mathcal{G}$, $\Gamma^j \in \mathcal{G}$, and finally $R\cup W \in \mathcal{G}$ and $R \cap W \in \mathcal{G}$. So, ($\Gamma^j , \mathcal{G}$) is a $\sigma$-algebra.} where $\mathcal{G}$ is a collection of sets of $\Gamma^j$ \cite{billingsley2008probability}. We can define a probability measure, $\mathcal{P} (\cdot)$, such that $\mathcal{P}(\Gamma^j)=1$ and $\mathcal{P}(\phi)=0$, then the triple, ($\Gamma^j ,\mathcal{G}, \mathcal{P}$) forms a probability space.
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.4\textwidth]{path.png}
\end{center}
\caption{Figure shows numerous trajectories of a particle $\omega_i\in\Omega$, and the action $\mathcal{A}_{ijk}$ (for $i^{th}$ particle) undergoing $j^{th}$ process along $k^{th}$ trajectory, $\gamma_k$ \cite{atanu}.}
\end{figure}
\section{Mathematical Proof}
As stated above, every element in $\Gamma^j$ is a possible trajectory, and it is not surprising that the set is unbounded because there exist infinitely many possible trajectories or sequence of trajectories for a given outcome (Feynman's path integral formulation argument) \cite{feynman2013feynman}. Furthermore, we associate some \emph{finite} (however small or large) probability to every trajectory (trajectories are basically the worldlines of a particular outcome). Thus, without any loss of generality, we can assume that the set, $W$, is \emph{countably} infinite, whereas the set, $R$, is finite. It is even more interesting to note that the set, $W$ is, in fact, \emph{uncountable}. Since the elements of the set $W$ represent all those trajectories which are thermodynamically infeasible, we can make any arbitrary combination of $\gamma^w_{k}$ (or a mix of $\gamma^w_{k}$ and $\gamma^r_{k}$, because it doesn't matter if a process is thermodynamically feasible in a later stage if it is infeasible in the very first instant!), and always come up with a new infeasible trajectory, say $\gamma^{w}_{k1}, \gamma^w_{k2}, ...$. As the set $W$ is infinite, we come up with an argument similar to Cantor's Diagonalization argument for uncountability of Reals \cite{halmos1960naive}, and prove that $W$ is uncountable because for every increment in $k$, we can always find a new element, $\gamma^w_{kt} \in W$ based on the combination of the diagonal elements (in this case, $\gamma^w_{kt}=\gamma^w_{1}\gamma^w_{3}\gamma^w_{5} ... \dotso \gamma^w_{\theta}$) and this will happen infinitely many times, \emph{i.e.},
\begin{eqnarray*}
\gamma^w_{k1}&=&\underline{\gamma^w_{1}}\gamma^w_{2}\gamma^w_{3}...\\
\gamma^w_{k2}&=&\gamma^w_{2}\underline{\gamma^w_{3}}\gamma^w_{4}...\\
\gamma^w_{k3}&=&\gamma^w_{3}\gamma^w_{2}\underline{\gamma^w_{5}}...\\
\vdots\\
\gamma^w_{kt}&=&\gamma^w_{9}\gamma^w_{7}\gamma^w_{6} ... \dotso \underline{\gamma^w_{\theta}}
\end{eqnarray*}
For the physical explanation of the above argument, we consider the case of a chemical reaction where substrates are converted into products. In order that the reaction is thermodynamically feasible, a minimum energy barrier (energy hill) called Activation Energy, $E_a$, has to be surpassed \cite{atkins2005elements}. Those trajectories along which this energy threshold is successfully realized can be visualized as the `right' trajectories and vice-versa. Even the so called `right' trajectories have further thermodynamic constraints, so basically, we can refine our intuition and come up with a finite number of possibilities, or maybe in some exceptional cases, \emph{atmost} countable `right' possibilities and \emph{always} uncountable `wrong' possibilities.
\subsection{Transformed Probability Space and Lebesgue Measure}
Let $\mathcal{T}_\mathcal{A}(\cdot)$ be an action map such that, $\mathcal{T}_\mathcal{A}=\frac{\mathcal{A}_{ijk}}{\max_k\mathcal{A}_{ijk}}$, in this case, $\mathcal{T}_\mathcal{A}(\gamma_k)=\frac{\mathcal{A}_{k}}{\max_k\mathcal{A}_{k}}=\frac{\int_{j,\gamma_k}p_i dq_i}{\max_k\int_{j,\gamma_k}p_i dq_i}$, and $\mathcal{T}_\mathcal{A}: \Gamma^j \rightarrow [0,1]$. The function, $\mathcal{T}_\mathcal{A}(\cdot)$, takes every trajectory, $\gamma_k$ from the set $\Gamma^j$, and calculates the action along that trajectory to the one along which the magnitude of action is maximum, which forms the transformed set, $\Gamma^{j\star}$. In fact, the map, $\mathcal{T}_\mathcal{A} (\cdot)$ is a homeomorphism\footnote{The inverse image of $\Gamma^{j\star}$ are the trajectories in $\Gamma^j$. The action map is similar to the pay-off function in Game Theory \cite{anttila2011natural} and Sigma-Profiles in multi-agent complex system modeling \cite{gershenson2011sigma}.}, $\mathcal{T}_\mathcal{A} (\cdot):\Gamma^j\rightarrow\Gamma^{j\star}$ (see figure 2). From our previous arguments, there are certain trajectories which are absolutely infeasible (although having a finite probability, however small), $\mathcal{A}_k \sim \max_k\mathcal{A}_k$ and few feasible ones, $\mathcal{A}_k \ll \max_k\mathcal{A}_k$. Thus, the transformed sample space, $\Gamma^{j\star}$ would take up all the values in the compact set, $\Gamma^{j\star}=[0,1] \subset \mathcal{R}$. Let $\mathcal{G}^\star$ be a collection of subsets of $\Gamma^{j\star}$, then the tuple $(\Gamma^{j\star},\mathcal{G}^\star)$ forms a $\sigma$-algebra. The interval, $[0,1]$ is Borel-measurable and the probability measure, $\lambda(\cdot)$ in this case is the Lebesgue measure; the triple, $(\Gamma^{j\star},\mathcal{G}^\star,\lambda)$ forms a probability space with $\lambda(\Gamma^{j\star})=1$, and $\lambda(\phi)=0$. From our previous argument based on the countability and uncountabilty of `right' and `wrong' trajectories, we can infer that the `right' trajectories in $\Gamma^{j}$ will cluster around $\{0\} \in [0,1]$, and would be countable in number, whereas the `wrong' trajectories would cluster around $\{1\} \in [0,1]$, and would be uncountably many. The subset, $R$ in the transformed sample space, $\Gamma^{j\star}$ is denoted by $R^{\star}$. Similarly, the subset $W$ is denoted by $W^{\star}$, and $\{0\}\in R^\star$ and $\{1\}\in W^\star$. On assigning the probability measure to the sets, $R^\star$ and $W^\star$, we observe $\lambda(R^{\star})=0$, since $R^{\star}$ is countable, and $\lambda(W^{\star})=1$, as $W^{\star}$ is uncountable \cite{billingsley2008probability}. What we just proved above is that no matter what, the probability that a `wrong' trajectory is chosen is one, whereas the probability of choosing a `right' trajectory is zero, which in fact, supports Murphy's statement!
\begin{figure}[b]
\begin{center}
\includegraphics[width=0.7\textwidth]{right.png}
\end{center}
\caption{Figure shows a homeomorphic relation between the two sets, $\Gamma^j$ and $\Gamma^{j\star}$.}
\end{figure}
There is a fallacy in the above argument. When we say there is a finite probability of any event to happen in some way, there is \emph{also} a finite probability that the event might happen in some other way, that we and also Murphy's statement ignore. As with any natural process, it can happen in infinitely many ways (as seen above), and the \emph{few} so-called `correct' trajectories get `lost' in the multitude of the uncountably many `wrong' ones, in spite of having finite chances of occurring.
\subsection{Dis-proving Murphy's Statement}
\noindent We define two sequences, $\{\gamma^{r}_{n}\} \in R$ and $\{\gamma^{w}_{n}\} \in W$. The sequence, $\{\gamma^{r}_{n}\}$ lies in the neighborhood of $\{0\} \in R^\star$, and the sequence, $\{\gamma^{w}_{n}\}$, in the neighborhood of $\{1\} \in W^\star$ ($\mathcal{T}_\mathcal{A}(\{\gamma^r_n\}\rightarrow\gamma^r) = 0$ and $\mathcal{T}_\mathcal{A}(\{\gamma^w_n\}\rightarrow\gamma^w)= 1$). From our initial assumption, $\mathcal{P}(\gamma^{r}_{k})$ and $\mathcal{P}(\gamma^{w}_{k})$ are both finite for \emph{every} $k$ in the respective sequences, $\{\gamma^{r}_n\}$ and $\{\gamma^{w}_n\}$. Based on our definition for `right' and `wrong' events, the terms in the sequences, $\gamma^{w}_{k}$ and $\gamma^{r}_{k}$ will again be infinitely many. Following our intuitive idea, we can easily foresee the probability, for every element in the sequence, $\gamma^{w}_{k}$ takes the form, $\mathcal{P}(\{\gamma^{w}_{n}\}) \sim n^{-\alpha}, \alpha \geq 1$, for every $\gamma^{w}_{k} \in \{\gamma^{w}_{n}\}$ ($n^{-\alpha}$ is not a distribution function, strictly because for some processes, the probability might decay, as $n^{-2}$, while for others, as $n^{-4}$, and so on). We can see that however small the probability for any process to occur along a specific trajectory, $\gamma^{w}_{k}$ be, it is finite. Since the probability for every trajectory is decaying with an increment in $k$, the summation of the probabilities, $\sum^{\infty}{\mathcal{P}(\gamma^{w}_{k}) < \infty}$ (because as we go nearer to the trajectories in the neighborhood of $\{1\}\in W^\star$, their probability will decrease at an even faster rate, making the probability of the sum of infinitely many such trajectories finite). From the Borel-Cantelli lemma, we can say that the probability of occurrence of these trajectories infinitely many times will be zero, and there are infinitely many of them \cite{billingsley2008probability}. So, the probability of their occurrence is zero. Conversely, for the sequence, $\{\gamma^{r}_{n}\}$, the probability distribution takes the form of $\mathcal{P}(\{\gamma^{r}_{n}\}) \sim n^{-\beta}$, where $\beta \sim 1$ for some $n$, while the probability is \emph{exactly} $1$ for finitely many trajectories. Thus, the summation of probabilities in this case, $\sum^{\infty}{\mathcal{P}(\gamma^{r}_{k})}$ diverges because $n^{-1}$ diverges. From the converse of the lemma, the probability of occurrence of these trajectories infinitely often will be one, and there are infinitely many of them and that such a trajectory (from the set, $R$) will be chosen \emph{almost surely}. We reflect on our findings based on the logic behind Murphy's statement, $\frac{p \rightarrow q, \enspace \neg q}{\therefore, \enspace \neg p}$ and observe that although `wrong' trajectories have finite probability, they do not occur, which is a contradiction to the proposition that, $\neg q \Rightarrow \neg p$ (with finite probabilities).

\section{Discussion and Conclusion}
\noindent Murphy's statement often serves as an anecdote in many real-life circumstances. Due to its affinity towards the negative aspects of any situation, it has often been misinterpreted to hold a resemblance with the Second Law of Thermodynamics. In this paper, we actually tried to relate both these statements and found quite striking results. As we mentioned earlier, Murphy's statement holds several philosophical connotation,s and appears as a generalized take on various things happening around us all the time. In order to test its scientific validity, we need to construct it as a chain of mathematically consistent statements and derive some logical inference from them. Any proof is invalid if the statements to be proved or disproved are themselves logically incoherent and mathematically not sound. Once this is achieved (in section 2.1), we need to outline a proper proof strategy. This is where the complexity associated with Murphy's statement increases, the reasons being: (\emph{i}) it deals with physical, observable processes, (\emph{ii}) the outcomes of the processes could be anything, and depending upon external perturbations the outcomes and even the processes may change their course of action, and (\emph{iii}) every outcome is plausible, \emph{i.e.}, every process has a finite chance of occurring and along numerous pathways leading to numerous finite possible outcomes. In order to deal with these uncertainties, we made assumptions like: any process or any outcome always has some energy associated with it, and any observable outcome can be achieved in several ways, which we call trajectories or pathways. In order to deal with a complex situation like this, we look for inspiration in the most fundamental laws of nature, and thereby propose that those outcomes, which are thermodynamically feasible shall be achieved along those trajectories that minimize action and maximize entropy. Using mathematical tools from probability theory and logic, we claim to disprove Murphy's statement in this paper. During the course of disproving the statement, we clearly showed the fallacy which is often the mis-interpreted truth behind the statement (section 3.1). We saw that there always exists an equivalent finite chance of an event to occur in different ways when it occurs in a certain way, which Murphy's statement generally ignores. We would further like to add that the intensifier, `definitely' in the construction of Murphy's statement enforces the fact that if something can go wrong, it will go wrong for sure. But we have seen earlier that the probability of something happening along the wrong way, as well as along the right way is finite, which raises a situation of \emph{post-hoc} fallacy in the statement. This type of logical fallacy is known as \emph{post-hoc ergo propter hoc} \cite{damer2012attacking}. A similar situation is observed in the case of the Law of Large Numbers and the Central Limit Theorem. According to the Law of Large Numbers (LLN), given sufficiently large number of trials, the average of the results thus obtained, shall converge (almost surely) to the expected value or mean of the experiment. Similarly, the Central Limit Theorem (CLT) states that the arithmetic mean of a sufficiently large number of iterates of independent random variables (each with a well-defined mean and variance) will be approximately normally distributed, regardless of the underlying distribution \cite{billingsley2008probability}. Both CLT and LLN seem to miss out on a crucial piece of information, \emph{i.e}, the information of the \emph{individual} outcomes or the \emph{individual} distributions respectively! Thus, Murphy's statement seems to be a far-fetched generalization of CLT and LLN in daily life. However, we must not forget that Murphy's statement has a physical essence associated with it (\emph{process}), and therefore needs to be analyzed on a different scale. Since it also deals with large number of outcomes and events, a statistical approach does sound good but a statistical formulation often fails to capture certain fine-tuned characteristics of the system. Our methodology, thus holds strong promises in the study of complex systems, where every motion of a particle counts and little perturbations may drive the system from an organized state to a state of chaos \cite{bak1997nature, bak1987self, heylighen2001cybernetics}. During the formulation of Murphy's statement, we also focused more on the process and the trajectories, and less on the system and its constituting elements. It is not necessary for the system or a system element to be a physical particle, rather it is immaterial and irrespective of the physical definition (or topology) of the system for our reasoning to hold. By a process, we refer to the change in state of the system, which is the same as the change in outcome(s) due to the system undergoing a process in time. As far as our assumption -  any process is associated with energy flows on the energy landscape - holds good, the central idea of the paper will hold true. Further, we mathematically prove that nature orders the occurrence of processes and their outcomes based on energy and entropy considerations. Thus a reinterpretation of Murphy's statement would read as: irrespective of all the possible outcomes and their qualitative aspects of `rightness' and `wrongness' anything that shall happen will definitely happen along the way that must minimize action.

\section{Acknowledgment}
\noindent The author would like to thank the reviewers for their valuable feedback highlighting the merits and drawbacks of this paper. Their constructive criticism and suggestions have made this paper much more informative and interesting to read. The author would also like to thank his colleagues, Sheetal Surana and Catherine Nipps, for their inputs and lengthy discussions on the topic of this paper.


\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.