We consider a cloud of dust with constant density $\rho$, select a ball of radius $a$ and consider the motion of a particle with mass $m$ on its surface. Its velocity is $\dot{a}$. From the energy conservation law \[E= \frac{m\dot{a}^2}{2} -\frac{4\pi G \rho a^2 m}{3}\] and \[\left(\frac{\dot{a}}{a}\right)^2 = \frac{8\pi G}{3}\rho - \frac{k}{a^2},\] where $k = 2E/m$. Therefore the first Friedman equation corresponds to the energy conservation law in the Newtonian interpretation.

Consider a uniform ball of radius $a$ with density $\rho$. Newton's law of motion for a point particle of mass $m$ on the surface reads: \[m\ddot a =- \frac{Gm}{a^2}\left(\frac{4\pi a^3\rho }{3} \right) \quad\Rightarrow\quad \frac{\ddot a}{a} = - \frac{4\pi G}{3}\rho.\] The difference with the second Friedman equation is in the absence of pressure $p$, which is a purely relativistic effect. We stress that Newtonian gravity can only lead to decelerated expansion.

Taking into account the homogeneity of the Universe in the continuity equation for the ideal liquid \[\frac{\partial \rho (t)}{\partial t} + \mbox{div}\left( \rho V \right) = 0\] and using the Hubble's law $V = HR$ to obtain the conservation law $\dot \rho + 3H\rho = 0$.

The considered Newtonian equations take the following form: equation for gravitational potential \[\Delta \Phi = 4\pi G \rho,\] continuity equation \[\dot{\rho}+\mbox{div}(\rho\vec{v})=0,\] equation of motion \[\dot{\vec{v}} + \left(\vec{v}\cdotp\nabla\right)\vec{v}= -\frac{\nabla p}{\rho}- \nabla \Phi.\] The stationarity and uniformity conditions for the Universe imply that $\vec{v}=0$, while $\rho$ and $p$ are constant in time and space. Then immediate substitution into the equation of gravitational potential allows one to see that \[\Phi = \frac{2}{3}\pi G \rho \left(\vec{r}\cdotp\vec{r}\right) +\vec{C}\vec{r} + K\] represents a solution, where $\vec{C}$ and $K$ are arbitrary integration constants. The continuity equation also reduces to identity. But the equations of motion (reduced to the condition $\nabla\Phi=0$) can be satisfied by no choice of $\vec{C}$ and $K$. Therefore there are no solutions of the system with required properties.

Let us consider the mass distribution of a collection of particles of masses $m_i$ at coordinate positions $\vec{x}_i$. We neglet nongravitational interaction, which assumes the fluctuations from homogeneity are appreciable only on scales small compared to the Hubble length, and peculiar motion is nonrelativistic. In this limit the Lagrangian for the system is \[L = \frac{1}{2}\sum\limits_i {{m_i}{a^2}\dot{\vec{x}_i^2} - MV}. \] The sum is over the kinetic energy of the particles in some large comoving volume. The total mass in this volume is $M=\sum_{i} m_{i}$. The gravitational potential energy per unit mass is \[V = - \frac{1}{2}\frac{G{a^5}}{M} \int {{d^3}{x_1}} {d^3}{x_2}\frac{\left[ {\rho \left( {x_1} \right)-{\rho _b}} \right] \left[ {\rho \left( {x_2} \right) - {\rho _b}} \right]} {x_{12}}.\] It is determined by the difference between the local mass density $\rho(\vec{x})$ and its mean value $\rho_b$. The momentum of the $i$-th particle \[{\vec p_i} = \frac{\partial L} {\partial \dot{\vec{x}}_i} = {m_i}{a^2}{\dot{\vec{x}}_i}\] and the Hamiltonian is \[H = \sum\limits_i {\frac{\vec p_i^2}{2{m_i}{a^2}}} + MV = M(K + V).\] We introduce the dimensionless mass autocorrelation function \[\xi \left( {\left| {{{\vec r}_1} - \vec r} \right|} \right) = \frac{\left[ {\rho(x_1) - {\rho _b}} \right] \left[\rho(x_2) -\rho_b\right]}{\rho _b^2}.\] With the change of variables to $\vec x = {\vec x_1} - {\vec x_2}$ and using $M = {a^3}\int {{d^3}x{\rho _b}}$, for the potential energy we obtain \begin{align*} V &= - \frac{1}{2}G{\rho _b}{a^2} \int {{d^3}x\frac{\xi (x)}{x}} = - 2\pi G{\rho _b}{J_2},\\ J_2&\equiv \int\limits_0^\infty {\xi (r)dr} \end{align*} Then from \[\frac{dH}{dt} = \frac{\partial H}{\partial t}\] where the partial derivative is computed at fixed particle coordinates and momenta. At fixed $\vec{x}_i$ and $\vec{p}_i$ we have $K\sim a^{-1}$ and $V \sim {\rho _b}{a^2} \propto {a^{ - 1}}$. So the last equation is reduced to \[\frac{dE}{dt}= \frac{d}{dt}(K + V) = - \frac{\dot a}{a}(2K + V).\] This equation is known as the cosmic energy or the Layzer-Irvine equation$^*$. $^*$It was derived independently in W.~M.~Irvine, Ph.D. thesis, Harvard University, (1961) and D.~Layzer, Astrophys. J. 138, 174 (1963).

Let a system's potential energy be a homogeneous function of coordinates of order $k$. Then if its motion is restricted to a final volume, the time averaged values of its kinetic and potential energies are related through the virial theorem \[2\left\langle K \right\rangle = k\left\langle V \right\rangle.\] For gravity $k=-1$, so \[2\left\langle K \right\rangle = - \left\langle V \right\rangle.\] The virial equation works only for the averages over large time intervals, when the system is in equilibrium. In this case $K$ and $V$ are nearly independent on time and the Layzer-Irvene equation becomes \[2K + V = 0,\] which is the Newtonian virial theorem.