Steady state in a gas of inelastic rough spheres heated by a uniform stochastic force

Let us consider a homogeneous dilute granular gas of identical hard spheres with mass $m$, diameter $\sigma$, and moment of inertia $I$, subject to a stochastic volume force $\mathbf{F}^{\text{wn}}$ (also called thermostat) with the properties of a Gaussian white noise.Williams and MacKintosh (1996); Williams (1996); Swift et al. (1998); van Noije and Ernst (1998); Montanero and Santos (2000); Gradenigo et al. (2011) This kind of forcing can model, for example, the energy input to grains immersed in a gas in turbulent flow.Ojha1 et al. (2004)

We denote the translational and angular (or rotational) particle velocities with $\mathbf{v}$ and $\bm{\omega}$ respectively. Binary collisions are characterized here by two material parameters, namely the coefficient of normal restitution $\alpha$ and the coefficient of tangential restitution $\beta$, which determine the shrinking of the normal and tangential component, respectively, of the relative velocity of the two surface points at contact. They are defined by (see, for instance, Refs. Kremer, 2010; S, 10; Kremer, Santos, and Garzó, 2014)

$$\widehat{\bm{\sigma}}\cdot\mathbf{u}^{\prime}=-\alpha(\widehat{\bm{\sigma}}\cdot\mathbf{u}),\quad\widehat{\bm{\sigma}}\times\mathbf{u}^{\prime}=-\beta(\widehat{\bm{\sigma}}\times\mathbf{u}).$$

(1)

Here, the primes denote post-collisional values, $\widehat{\bm{\sigma}}$ is the unit collision vector joining the centers of the two colliding spheres (and pointing from the center of particle 1 to the center of particle 2), and

$${\mathbf{u}=\mathbf{v}_{1}-\mathbf{v}_{2}-{\frac{\sigma}{2}}\widehat{\bm{\sigma}}\times(\bm{\omega}_{1}+\bm{\omega}_{2})}$$

(2)

is the relative velocity of the points of the spheres which are in contact during a binary encounter.

The coefficient of normal restitution $\alpha$ takes values between $0$ (completely inelastic collision) and $1$ (completely elastic collision), while the coefficient of tangential restitution takes values between $-1$ (completely smooth collision, implying that rotational velocities are unchanged) and $1$ (completely rough collision).Chapman and Cowling (1970); Kremer (2010) Equation (1) supplements the linear and angular momentum conservation laws to yield the collision rules,Kremer (2010); Santos, Kremer, and dos Santos (2011); Brilliantov and Pöschel (2004); Kremer, Santos, and Garzó (2014)

$${m\mathbf{v}_{1,2}^{\prime}=m\mathbf{v}_{1,2}\mp\mathbf{Q},\quad I\bm{\omega}_{1,2}^{\prime}=I\bm{\omega}_{1,2}-{\frac{\sigma}{2}}\widehat{\bm{\sigma}}\times\mathbf{Q},}$$

(3)

where the impulse exerted by particle 1 on particle 2 is given by

$${\mathbf{Q}={m\widetilde{\alpha}}(\widehat{\bm{\sigma}}\cdot\mathbf{u})\widehat{\bm{\sigma}}-m\widetilde{\beta}\widehat{\bm{\sigma}}\times(\widehat{\bm{\sigma}}\times\mathbf{u}).}$$

(4)

Here, we have introduced the abbreviations

$${\widetilde{\alpha}\equiv{\frac{1+\alpha}{2}},\quad\widetilde{\beta}\equiv{\frac{1+\beta}{2}}{\frac{\kappa}{\kappa+1}},\quad\kappa\equiv{\frac{4I}{m\sigma^{2}}}.}$$

(5)

While collisions (except if $\alpha=1$ and $\beta=\pm 1$) dissipate kinetic energy (translational plus rotational), the homogeneous stochastic force $\mathbf{F}^{\text{wn}}$ injects translational kinetic energy to grains. It has properties of a Gaussian white noise, i.e.,

	$\displaystyle\langle{\bf F}_{i}^{\text{wn}}(t)\rangle=$	$\displaystyle{\bf 0},$		(6a)
	$\displaystyle\langle{\bf F}_{i}^{\text{wn}}(t){\bf F}_{j}^{\text{wn}}(t^{\prime})\rangle=$	$\displaystyle\mathsf{I}m^{2}\chi_{0}^{2}\delta_{ij}\delta(t-t^{\prime}),$		(6b)

where indexes $i,j$ refer to particles, $\mathsf{I}$ is the $3\times 3$ unit matrix, and $\chi_{0}^{2}$ measures the characteristic strength of the stochastic force.

In homogeneous states, the Boltzmann equation corresponding to the stochastic external force $\mathbf{F}^{\text{wn}}$ becomesvan Noije and Ernst (1998); Montanero and Santos (2000)

$$\frac{\partial f(\mathbf{v},\bm{\omega};t)}{\partial t}-\frac{\chi_{0}^{2}}{2}\left(\frac{\partial}{\partial\mathbf{v}}\right)^{2}f(\mathbf{v},\bm{\omega};t)={J[\mathbf{v},\bm{\omega}|f(t)]},$$

(7)

where $f(\mathbf{v},\bm{\omega};t)$ is the velocity distribution function and $J[\mathbf{v},\bm{\omega}|f]$ is the collision integral in the (inelastic) Boltzmann equation for rough spheres, which accounts for the collision rules (3) and (4).

Given a certain function $A(\mathbf{v},\bm{\omega})$, we define its average as

$$\langle A(t)\rangle=\frac{1}{n}\int d\mathbf{v}\int d\bm{\omega}\,A(\mathbf{v},\bm{\omega})f(\mathbf{v},\bm{\omega};t),\quad n=\int d\mathbf{v}\int d\bm{\omega}\,f(\mathbf{v},\bm{\omega};t),$$

(8)

where $n$ is the number density. By Galilean invariance, we choose $\langle\mathbf{v}\rangle=\mathbf{0}$. Moreover, we assume isotropic states, so that $\langle\bm{\omega}\rangle=\mathbf{0}$.Kremer, Santos, and Garzó (2014) Thus, the basic quantities are the translational ($T_{t}$), rotational ($T_{r}$), and total ($T$) granular temperatures:

$${T_{t}={\frac{m}{3}}\langle v^{2}\rangle,\quad T_{r}={\frac{I}{3}}\langle\omega^{2}\rangle,\quad\quad T=\frac{T_{t}+T_{r}}{2}=T_{t}\frac{1+\theta}{2}.}$$

(9)

where the temperature ratio $\theta\equiv T_{r}/T_{t}$ is an important quantity in this system because, as we will see, its steady-state value is independent of the level of forcing ($\chi_{0}^{2}$).

The evolution equations for the granular temperatures $T_{t}$ and $T_{r}$ are obtained by just multiplying both sides of Eq. (7) by the particle translational and rotational kinetic energies, respectively, and integrating over all velocity values. This gives

$$\partial_{t}T_{t}-{m\chi_{0}^{2}}=-\xi_{t}T_{t},\quad{\partial_{t}T_{r}=-\xi_{r}T_{r}},\quad\partial_{t}T-\frac{m\chi_{0}^{2}}{2}=-\zeta T$$

(10)

Here, the translational ($\xi_{t}$) and rotational ($\xi_{r}$) production rates and the cooling rate ($\zeta$) are defined by the collision integralsSantos, Kremer, and dos Santos (2011); Kremer, Santos, and Garzó (2014)

		$\displaystyle\xi_{t}=-{\frac{m}{3nT_{t}}}\int d\mathbf{v}\int d\bm{\omega}\,v^{2}{J[\mathbf{v},\bm{\omega}|f]},\quad\xi_{r}=-{\frac{I}{3nT_{r}}}\int d\mathbf{v}\int d\bm{\omega}\,\omega^{2}{J[\mathbf{v},\bm{\omega}|f]},$		(11a)
		$\displaystyle\zeta=\frac{\xi_{t}T_{t}+\xi_{r}T_{r}}{2T}=\frac{\xi_{t}+\xi_{r}{\theta}}{1+\theta}.$		(11b)

For convenience, let us rewrite the Boltzmann equation (7) in dimensionless form. For this, we first define the reduced velocity distribution function

$$\phi(\mathbf{c},\mathbf{w};\tau)\equiv\frac{1}{n}\left(\frac{4T_{t}(t)T_{r}(t)}{mI}\right)^{3/2}f(\mathbf{v},\bm{\omega};t),$$

(12)

where the reduced particle velocities are

$$\mathbf{c}(t)\equiv\frac{\mathbf{v}}{\sqrt{2T_{t}(t)/m}},\quad\mathbf{w}(t)\equiv\frac{\bm{\omega}}{\sqrt{2T_{r}(t)/I}}.$$

(13)

In addition, time is measured by the parameter $\tau$ defined by

$${\tau(t)=\int_{0}^{t}dt^{\prime}\,\nu(t^{\prime}),\quad\nu(t)=2n\sigma^{2}\sqrt{\pi T_{t}(t)/m},}$$

(14)

where $\nu(t)$ is the (time-dependent) collision frequency. Thus, $\tau(t)$ measures the accumulated number of collisions per particle up to time $t$. With this time scale, we can measure relaxation times to the steady state that are physically meaningful, since they coincide with the aging times to hydrodynamics.Vega Reyes, Santos, and Kremer (2014)

The moments of the reduced distribution function $\phi(\mathbf{c},\mathbf{w};\tau)$ are

$$M_{pq}^{(r)}(\tau)\equiv\langle c^{p}w^{q}(\mathbf{c}\cdot\mathbf{w})^{r}\rangle=\int d\mathbf{c}\int d\mathbf{w}\,c^{p}w^{q}(\mathbf{c}\cdot\mathbf{w})^{r}\phi(\mathbf{c},\mathbf{w};\tau),$$

(15)

with $p,q,r=\text{even}$, by symmetry. Note that $M_{pq}^{(r)}$ is a moment of degree $p+q+2r$ and $M_{0,0}^{(0)}=1$, $M_{2,0}^{(0)}=M_{0,2}^{(0)}=\frac{3}{2}$. The dimensionless measure of the noise intensity is defined by

$$\gamma(\tau)\equiv\frac{3}{2}\frac{\chi_{0}^{2}}{\nu(t)T_{t}(t)/m}=\frac{3\chi_{0}^{2}}{4\sqrt{\pi}n\sigma^{2}[T_{t}(t)/m]^{3/2}}.$$

(16)

Notice that, although $\chi_{0}^{2}$ is a constant, its dimensionless counterpart $\gamma$ is time-dependent due to the scaling with the temperature $T_{t}$. Of course, once a steady state is reached, $\gamma$ becomes constant too. Note also that $\gamma^{-2/3}$ can be seen as the translational temperature $T_{t}$ in units of a reference temperature $T_{\text{ref}}=m(3\chi_{0}^{2}/4\sqrt{\pi}n\sigma^{2})^{2/3}$ defined from the white noise parameter $\chi_{0}^{2}$.

Let us finally define the reduced collisional moments

$$\mu_{pq}^{(r)}(\tau)\equiv-\int d\mathbf{c}\int d\mathbf{w}\,c^{p}w^{q}(\mathbf{c}\cdot\mathbf{w})^{r}{\mathcal{J}[\mathbf{c},\mathbf{w}|\phi(\tau)]},$$

(17)

where

$$\mathcal{J}[\mathbf{c},\mathbf{w}|\phi(\tau)]=\frac{1}{n\nu(t)}\left[\frac{4T_{t}(t)T_{r}(t)}{mI}\right]^{3/2}J[\mathbf{v},\bm{\omega}|f(t)]$$

(18)

is the dimensionless form of the collision integral. From the definitions (11a) and (17), it is easy to see that $\xi_{t}=\frac{2}{3}\nu\mu_{20}^{(0)}$, $\xi_{r}=\frac{2}{3}\nu\mu_{02}^{(0)}$. The evolution equations for the temperature ratio $\theta(\tau)$ and the reduced noise strength $\gamma(\tau)$ can be obtained from Eq. (10) as

$$\frac{\partial\ln\theta(\tau)}{\partial\tau}=\frac{2}{3}\left[\mu_{20}^{(0)}(\tau)-\mu_{02}^{(0)}(\tau)-\gamma(\tau)\right],\quad\frac{\partial\ln\gamma(\tau)}{\partial\tau}=\mu_{20}^{(0)}(\tau)-\gamma(\tau).$$

(19)

Note that, since $\mathbf{v}$ and $\bm{\omega}$ are scaled each with a different temperature in Eq. (13), the collision rules (3) and (4), when expressed in terms of $\mathbf{c}$ and $\mathbf{w}$, depend parametrically on $\theta$.Santos, Kremer, and dos Santos (2011) This dependence is obviously transferred to the collisional moments (17).

Taking into account the previous dimensionless quantities, and making use again of the temperature evolution equations (10), one can obtain the dimensionless version of the Boltzmann equation (7) as

	$\displaystyle\frac{\partial\phi(\mathbf{c},\mathbf{w};\tau)}{\partial\tau}$	$\displaystyle+\frac{\mu_{20}^{(0)}(\tau)-\gamma(\tau)}{3}\frac{\partial}{\partial\mathbf{c}}\cdot\left[\mathbf{c}\phi(\mathbf{c},\mathbf{w};\tau)\right]+\frac{\mu_{02}^{(0)}(\tau)}{3}\frac{\partial}{\partial\mathbf{w}}\cdot\left[\mathbf{w}\phi(\mathbf{c},\mathbf{w};\tau)\right]$
		$\displaystyle-\frac{\gamma(\tau)}{6}\left(\frac{\partial}{\partial\mathbf{c}}\right)^{2}\phi(\mathbf{c},\mathbf{w};\tau)={\mathcal{J}[\mathbf{c},\mathbf{w}|\phi(\tau)]}.$		(20)

Multiplying both sides by $c^{p}w^{q}(\mathbf{c}\cdot\mathbf{w})^{r}$ and integrating over $\mathbf{c}$ and $\mathbf{w}$, we get the moment hierarchy

	$\displaystyle\frac{\partial\ln M_{pq}^{(r)}(\tau)}{\partial\tau}=$	$\displaystyle\frac{p+r}{3}\left[\mu_{20}^{(0)}(\tau)-\gamma(\tau)\right]+\frac{q+r}{3}{\mu_{02}^{(0)}}(\tau)-\frac{\mu_{pq}^{(r)}(\tau)}{M_{pq}^{(r)}(\tau)}$
		$\displaystyle+\frac{\gamma(\tau)}{6}\frac{p(p+1+2r)M_{p-2,q}^{(r)}(\tau)+r(r-1)M_{p,q+2}^{(r-2)}(\tau)}{M_{pq}^{(r)}(\tau)}.$		(21)

This equation must be supplemented with Eq. (19). As will be seen later, the steady-state distribution function is independent of the precise value of the stochastic force intensity $\chi_{0}^{2}$ or, equivalently, of the initial value $\gamma(0)$, although of course the transient states are dependent on it.

In principle, the reduced distribution function $\phi(\mathbf{c},\mathbf{w})$ would be a function of six velocity components. On the other hand, since we restrict ourselves to isotropic states, $\phi$ is actually a function of three scalar quantities: the norms $c^{2}$ and $w^{2}$, and the squared dot product $(\mathbf{c},\mathbf{w})^{2}$.Santos, Kremer, and dos Santos (2011) Thus, it makes sense to represent $\phi(\mathbf{c},\mathbf{w})$ as an expansion in a complete set $\{\Psi_{jk}^{(\ell)}(\mathbf{c},\mathbf{w})\}$ of orthogonal polynomials,

$$\phi(\mathbf{c},\mathbf{w};\tau)=\phi_{M}(c,w)\sum_{j=0}^{\infty}\sum_{k=0}^{\infty}\sum_{\ell=0}^{\infty}a_{jk}^{(\ell)}(\tau)\Psi_{jk}^{(\ell)}(\mathbf{c},\mathbf{w}),\quad\phi_{M}(c,w)\equiv\pi^{-3}e^{-c^{2}-w^{2}},$$

(22)

where

$$\Psi_{jk}^{(\ell)}(\mathbf{c},\mathbf{w})\equiv L_{j}^{(2\ell+\frac{1}{2})}(c^{2})L_{k}^{(2\ell+\frac{1}{2})}(w^{2})\left(c^{2}w^{2}\right)^{\ell}P_{2\ell}(\cos\vartheta),\quad\cos\vartheta\equiv\frac{\mathbf{c}\cdot\mathbf{w}}{cw},$$

(23)

where $L_{j}^{(2\ell+\frac{1}{2})}(x)$ and $P_{2\ell}(x)$ are the Laguerre and Legendre polynomials, respectively.Abramowitz and Stegun (1972) Defining the scalar product of two arbitrary isotropic (real) functions $\Phi_{1}(\mathbf{c},\mathbf{w})$ and $\Phi_{2}(\mathbf{c},\mathbf{w})$ as

$$\langle\Phi_{1}|\Phi_{2}\rangle\equiv\int d\mathbf{c}\int d\mathbf{w}\,\phi_{M}(c,w)\Phi_{1}(\mathbf{c},\mathbf{w})\Phi_{2}(\mathbf{c},\mathbf{w}),$$

(24)

one can obtain the orthogonality relation

$$\langle\Psi_{jk}^{(\ell)}|\Psi_{j^{\prime}k^{\prime}}^{(\ell^{\prime})}\rangle=\mathcal{N}_{jk}^{(\ell)}\delta_{j,j^{\prime}}\delta_{k,k^{\prime}}\delta_{\ell,\ell^{\prime}},\quad\mathcal{N}_{jk}^{(\ell)}\equiv\frac{\Gamma(j+2\ell+\frac{3}{2})\Gamma(k+2\ell+\frac{3}{2})}{(\pi/4)(4\ell+1)j!k!}.$$

(25)

As a consequence, the coefficients in the expansion (22) are

	$\displaystyle a_{jk}^{(\ell)}=$	$\displaystyle\frac{\langle\Psi_{jk}^{(\ell)}\rangle}{\mathcal{N}_{jk}^{(\ell)}}$
	$\displaystyle=$	$\displaystyle\sum_{j_{1}=0}^{j}\sum_{k_{1}=0}^{k}\sum_{\ell_{1}=0}^{\ell}\frac{(-1)^{j_{1}+k_{1}+\ell_{1}}(\pi/4)j!k!(4\ell+1)(4\ell-2\ell_{1}-1)!!M_{2(j_{1}+\ell_{1}),2(k_{1}+\ell_{1})}^{(2\ell-2\ell_{1})}}{j_{1}!(j-j_{1})!\Gamma\left(j_{1}+2\ell+\frac{3}{2}\right)k_{1}!(k-k_{1})!\Gamma\left(k_{1}+2\ell+\frac{3}{2}\right)2^{\ell_{1}}\ell_{1}!(2\ell-2\ell_{1})!},$		(26)

where in the second step use has been made of the explicit expressions of the Laguerre and Legendre polynomials.Abramowitz and Stegun (1972) We observe that $a_{jk}^{(\ell)}$ is a linear combination of the moments $M_{pq}^{(r)}$ with $p,q,r=\text{even}$ and degree $2\ell\leq p+q+2r\leq 2(j+k+2\ell)$. By normalization, $a_{00}^{(0)}=1$, $a_{10}^{(0)}=a_{01}^{(0)}=0$. Therefore, the first nontrivial coefficients are those associated with moments of degree four, namely

	$$a_{20}^{(0)}=\frac{4}{15}\langle c^{4}\rangle-1,\quad a_{02}^{(0)}=\frac{4}{15}\langle w^{4}\rangle-1,$$		(27a)
	$$a_{11}^{(0)}=\frac{4}{9}\langle c^{2}w^{2}\rangle-1,\quad{a_{00}^{(1)}=\frac{8}{15}\left[\langle(\mathbf{c}\cdot\mathbf{w})^{2}\rangle-\frac{1}{3}\langle c^{2}w^{2}\rangle\right]}.$$		(27b)

These are the fourth-degree cumulants, which measure the basic deviations of the velocity distribution function $\phi(\mathbf{c},\mathbf{w})$ from the (two-temperaure) Maxwellian $\phi_{M}(c,w)$. While $a_{20}^{(0)}$ and $a_{02}^{(0)}$ measure the kurtosis of the (marginal) translational and rotational velocity distribution functions, respectively, the coefficient $a_{11}^{(0)}$ characterizes the scalar translational-rotational correlations, i.e., $\langle c^{2}w^{2}\rangle\neq\langle c^{2}\rangle\langle w^{2}\rangle$. Finally, $a_{00}^{(1)}$ informs us about the possible existence of orientational correlations, i.e., $\langle c^{2}w^{2}\cos^{2}\vartheta\rangle\neq\langle c^{2}w^{2}\rangle\langle\cos^{2}\vartheta\rangle$ and $\langle\cos^{2}\vartheta\rangle\neq\frac{1}{3}$. To disentangle the latter two effects, let us define the quantities

$${h\equiv\frac{5}{8}\left[\frac{\langle(\mathbf{c}\cdot\mathbf{w})^{2}\rangle}{\langle c^{2}w^{2}\rangle\langle\cos^{2}\vartheta\rangle}-1\right]},\quad b\equiv\frac{10}{3}\left(\langle\cos^{2}\vartheta\rangle-\frac{1}{3}\right).$$

(28)

The parameters $h$, $b$, $a_{11}^{(0)}$, and $a_{00}^{(1)}$ are related by

$${h=\frac{1}{16}\frac{25a_{00}^{(1)}-9b\left(1+a_{11}^{(0)}\right)}{\left(1+\frac{9}{10}b\right)\left(1+a_{11}^{(0)}\right)}.}$$

(29)

In previous works,Brilliantov et al. (2007); Kranz et al. (2009); Rongali and Alam (2014) the angle $\vartheta$ between translational and rotational particle velocities has been used for measuring orientational correlations in a granular gas of rough spheres, by tracking the quantity $\langle\cos^{2}\vartheta\rangle$ or, equivalently, $b$. The use of $b$ has the advantage over $a_{00}^{(1)}$ of being more intuitive (since it directly refers to the angle between translational and rotational particle velocities). However, it has the disadvantage that, since $\cos\vartheta=\mathbf{c}\cdot\mathbf{w}/cw$ is not a polynomial, the parameter $b$ involves an infinite number of coefficients in the expansion (22). More specifically, it is easy to find

$$b=\sum_{j=0}^{\infty}\sum_{k=0}^{\infty}a_{jk}^{(1)}.$$

(30)

Analogously, the quantity $h$ also involves an infinite number of coefficients. For this reason, we prefer instead to use the single coefficient $a_{00}^{(1)}$ for measuring the orientational correlations in a granular gas of rough spheres.Vega Reyes, Santos, and Kremer (2014)

Taking into account the general moment hierarchy (21), the evolution equations for the fourth-degree cumulants are

	$$\frac{\partial\ln\left[1+a_{20}^{(0)}(\tau)\right]}{\partial\tau}=\frac{4}{3}\mu_{20}^{(0)}(\tau)-\frac{4}{3}\gamma(\tau)-\frac{4}{15}\frac{\mu_{40}^{(0)}(\tau)-5\gamma(\tau)}{1+a_{20}^{(0)}(\tau)},$$		(31a)
	$$\frac{\partial\ln\left[1+a_{02}^{(0)}(\tau)\right]}{\partial\tau}=\frac{4}{3}\mu_{02}^{(0)}(\tau)-\frac{4}{15}\frac{\mu_{04}^{(0)}(\tau)}{1+a_{02}^{(0)}(\tau)},$$		(31b)
	$$\frac{\partial\ln\left[1+a_{11}^{(0)}(\tau)\right]}{\partial\tau}=\frac{2}{3}\mu_{20}^{(0)}(\tau)+\frac{2}{3}\mu_{02}^{(0)}(\tau)-\frac{2}{3}\gamma(\tau)-\frac{4}{9}\frac{\mu_{22}^{(0)}(\tau)-\frac{3}{2}\gamma(\tau)}{1+a_{11}^{(0)}(\tau)},$$		(31c)
	$$\frac{\partial\ln\left[1+a_{11}^{(0)}(\tau)+\frac{5}{2}a_{00}^{(1)}(\tau)\right]}{\partial\tau}=\frac{2}{3}\mu_{20}^{(0)}(\tau)+\frac{2}{3}\mu_{02}^{(0)}(\tau)-\frac{2}{3}\gamma(\tau)-\frac{4}{3}\frac{\mu_{00}^{(2)}(\tau)-\frac{1}{2}\gamma(\tau)}{1+a_{11}^{(0)}(\tau)+\frac{5}{2}a_{00}^{(1)}(\tau)}.$$		(31d)

In the steady state ($\partial_{\tau}\to 0$), Eqs. (10) and (31) imply the conditions

	$$\gamma=\mu_{20}^{(0)},\quad\mu_{02}^{(0)}=0,$$		(32a)
	$$\frac{1}{5}\mu_{40}^{(0)}=\frac{2}{3}\mu_{22}^{(0)}=2\mu_{00}^{(2)}=\mu_{20}^{(0)},\quad\mu_{04}^{(0)}=0.$$		(32b)

Before closing this section, and for further use, let us define the marginal distributions $\phi_{c}(c)$, $\phi_{w}(w)$, $\phi_{cw}(c^{2}w^{2})$, $\phi_{\mathbf{c}\cdot\mathbf{w}}((\mathbf{c}\cdot\mathbf{w})^{2})$, and $\phi_{\vartheta}(\cos^{2}\vartheta)$ from the joint velocity distribution function $\phi(\mathbf{c},\mathbf{w})$ as

	$$\phi_{c}(c)=4\pi c^{2}\int d\mathbf{w}\,\phi(\mathbf{c},\mathbf{w}),\quad\phi_{w}(w)=4\pi w^{2}\int d\mathbf{c}\,\phi(\mathbf{c},\mathbf{w}),$$		(33a)
	$$\phi_{cw}(c^{2}w^{2})=\int d\mathbf{c}^{\prime}\int d\mathbf{w}^{\prime}\,\delta(c^{\prime 2}w^{\prime 2}-c^{2}w^{2})\phi(\mathbf{c}^{\prime},\mathbf{w}^{\prime}),$$		(33b)
	$$\phi_{\mathbf{c}\cdot\mathbf{w}}((\mathbf{c}\cdot\mathbf{w})^{2})=\int d\mathbf{c}^{\prime}\int d\mathbf{w}^{\prime}\,\delta((\mathbf{c}^{\prime}\cdot\mathbf{w}^{\prime})^{2}-(\mathbf{c}\cdot\mathbf{w})^{2})\phi(\mathbf{c}^{\prime},\mathbf{w}^{\prime}),$$		(33c)
	$$\phi_{\vartheta}(\cos^{2}\vartheta)=\int d\mathbf{c}^{\prime}\int d\mathbf{w}^{\prime}\,\delta(\cos^{2}\vartheta^{\prime}-\cos^{2}\vartheta)\phi(\mathbf{c}^{\prime},\mathbf{w}^{\prime}).$$		(33d)

The less conventional marginal distribution functions (33b)–(33d) are directly related to the scalar and orientational correlations between the translational and rotational particle velocities. More specifically,

$$\langle c^{2}w^{2}\rangle=\int_{0}^{\infty}dx\,x\phi_{cw}(x),\quad\langle(\mathbf{c}\cdot\mathbf{w})^{2}\rangle=\int_{0}^{\infty}dx\,x\phi_{\mathbf{c}\cdot\mathbf{w}}(x),\quad\langle\cos^{2}\vartheta\rangle=\int_{0}^{1}dx\,x\phi_{\vartheta}(x).$$

(34)

The simplest approximation consists in neglecting all the cumulants $a_{jk}^{(\ell)}$, $j+k+2\ell\geq 2$, in Eq. (22), i.e., $\phi(\mathbf{c},\mathbf{w})\simeq\phi_{M}(c,w)$. This Maxwellian approximation is of course unable to account for the evolution and steady-state values of the cumulants. However, it can capture the main features of the temperature ratio $\theta(\tau)$ and the reduced noise intensity $\gamma(\tau)$. When $\phi(\mathbf{c},\mathbf{w})\simeq\phi_{M}(c,w)$ is inserted into Eq. (17), one obtains

	$$\mu_{20,M}^{(0)}=1-\alpha^{2}+\frac{\kappa(1+\beta)}{(1+\kappa)^{2}}\left[2+\kappa(1-\beta)-\theta(1+\beta)\right],$$		(35a)
	$$\mu_{02,M}^{(0)}=\frac{\kappa(1+\beta)}{(1+\kappa)^{2}}\left[2+\kappa^{-1}(1-\beta)-\theta^{-1}(1+\beta)\right].$$		(35b)

Using these expressions in (19), one gets a closed set of two coupled nonlinear differential equations that can be numerically solved for arbitrary initial values $\theta(0)$ and $\gamma(0)$. Regardless of the values of $\theta(0)$ and $\gamma(0)$, common steady-state values are reached after a few collisions per particle. They are explicitly obtained from Eq. (32a) as

$$\theta_{M}=\frac{1+\beta}{2+\kappa^{-1}(1-\beta)},\quad\gamma_{M}=1-\alpha^{2}+\frac{2(1-\beta^{2})}{2+\kappa^{-1}(1-\beta)}.$$

(36)

It is interesting to note that the Maxwellian prediction for the steady-state value of the temperature ratio is independent of the coefficient of normal restitution $\alpha$. Moreover, $\theta_{M}$ is always smaller than unity (i.e., $T_{r}<T_{t}$), except in the limit of completely rough spheres ($\beta=1$), where $\theta_{M}=1$. In the case of the reduced noise intensity, we can observe that $\gamma_{M}$ monotonically decreases as $\alpha$ increases at fixed $\beta$, as a consequence of a progressively smaller cooling effect. Because of the same reason, $\gamma_{M}$ presents a non-monotonic behavior with respect to $\beta$ at fixed $\alpha$, reaching a maximum value $\gamma_{M,\max}=1-\alpha^{2}+4\kappa(1+2\kappa-2\sqrt{\kappa+\kappa^{2}})$ at $\beta=1-2(\sqrt{\kappa+\kappa^{2}}-\kappa)$. However, as we will see in Sec. IV, $\theta$ actually depends on $\alpha$ but much more weakly than it depends on $\beta$. Also, $\theta$ can reach values slightly larger than unity if $\beta$ is very close to $1$ and $\alpha$ is small enough. In any case, as shown later, the Maxwellian expressions (36) are indeed rather accurate.

The marginal distribution functions defined by (33) adopt simple forms in the Maxwellian approximation,

	$$\phi_{c,M}(c)=\frac{4}{\sqrt{\pi}}c^{2}e^{-c^{2}},\quad\phi_{w,M}(c)=\frac{4}{\sqrt{\pi}}w^{2}e^{-w^{2}},$$		(37a)
	$$\phi_{cw,M}(c^{2}w^{2})=\frac{8cw}{\pi}K_{0}(2cw),\quad\phi_{\mathbf{c}\cdot\mathbf{w},M}((\mathbf{c}\cdot\mathbf{w})^{2})=\frac{4}{\pi}K_{1}(2|\mathbf{c}\cdot\mathbf{w}|),\quad\phi_{\vartheta,M}(\cos^{2}\vartheta)=\frac{1}{2|\cos\vartheta|},$$		(37b)

where $K_{n}(x)$ is a modified Bessel function of the second kind.Abramowitz and Stegun (1972) Upon deriving (37b) use has been made of the mathematical properties $\delta(c^{\prime 2}w^{\prime 2}-x)=\delta(w^{\prime}-\sqrt{x}/c^{\prime})/2c^{\prime 2}w^{\prime}$, $\delta((\mathbf{c}^{\prime}\cdot\mathbf{w}^{\prime})^{2}-x)=\delta(|\cos\vartheta|-\sqrt{x}/cw)/2c^{2}w^{2}|\cos\vartheta|$, and $\delta(\cos^{2}\vartheta-x)=\delta(|\cos\vartheta|-\sqrt{x})/2|\cos\vartheta|$, respectively. Apart from $\langle c^{4}\rangle=\langle w^{4}\rangle=\frac{15}{4}$, it is straightforward to check from Eqs. (34) and (37b) that $\langle c^{2}w^{2}\rangle=\frac{9}{4}$, $\langle(\mathbf{c}\cdot\mathbf{w})^{2}\rangle=\frac{3}{4}$, and $\langle\cos^{2}\vartheta\rangle=\frac{1}{3}$ in the Maxwellian approximation, as expected.

As a much more elaborate approximation that captures the most important non-Maxwellian features of the velocity distribution function, we truncate the exact expansion (22) after $j+k+2\ell=2$, i.e.Vega Reyes, Santos, and Kremer (2014)

	$\displaystyle\frac{\phi(\mathbf{c},\mathbf{w})}{\phi_{M}(c,w)}\simeq$	$\displaystyle{1+}a_{20}^{(0)}\Psi_{20}^{(0)}(\mathbf{c},\mathbf{w})+a_{02}^{(0)}\Psi_{02}^{(0)}(\mathbf{c},\mathbf{w})+a_{11}^{(0)}\Psi_{11}^{(0)}(\mathbf{c},\mathbf{w})+a_{00}^{(1)}\Psi_{00}^{(1)}(\mathbf{c},\mathbf{w})$
	$\displaystyle=$	$\displaystyle 1+{a_{20}^{(0)}}\frac{{15}-20c^{2}+4c^{4}}{8}+{a_{02}^{(0)}}\frac{{15}-20w^{2}+4w^{4}}{8}+{a_{11}^{(0)}}\frac{\left(3-2c^{2}\right)\left(3-2w^{2}\right)}{4}$
		$\displaystyle+{a_{00}^{(1)}}\frac{3(\mathbf{c}\cdot\mathbf{w})^{2}-c^{2}w^{2}}{2}.$		(38)

This approximation differs from that considered in Ref. Santos, Kremer, and dos Santos, 2011 by the addition of the term headed by the coefficient $a_{00}^{(1)}$. The associated marginal distributions are

	$$\frac{\phi_{c}(c)}{\phi_{c,M}(c)}=1+a_{20}^{(0)}\frac{{15}-20c^{2}+4c^{4}}{8},\quad\frac{\phi_{w}(w)}{\phi_{w,M}(w)}=1+a_{02}^{(0)}\frac{{15}-20w^{2}+4w^{4}}{8},$$			(39a)
	$\displaystyle\frac{\phi_{cw}(c^{2}w^{2})}{\phi_{cw,M}(c^{2}w^{2})}=$	$\displaystyle 1+\frac{a_{20}^{(0)}+a_{02}^{(0)}}{8}\left[15+4c^{2}w^{2}-16cw\frac{K_{1}(2cw)}{K_{0}(2cw)}\right]$
		$\displaystyle+\frac{a_{11}^{(0)}}{4}\left[9+4c^{2}w^{2}-12cw\frac{K_{1}(2cw)}{K_{0}(2cw)}\right],$		(39b)
	$\displaystyle\frac{\phi_{\mathbf{c}\cdot\mathbf{w}}((\mathbf{c}\cdot\mathbf{w})^{2})}{\phi_{\mathbf{c}\cdot\mathbf{w},M}((\mathbf{c}\cdot\mathbf{w})^{2})}=$	$\displaystyle 1+\frac{a_{20}^{(0)}+a_{02}^{(0)}}{8}\left[3+4(\mathbf{c}\cdot\mathbf{w})^{2}-12|\mathbf{c}\cdot\mathbf{w}|\frac{K_{0}(2|\mathbf{c}\cdot\mathbf{w}|)}{K_{1}(2|\mathbf{c}\cdot\mathbf{w}|)}\right]$
		$\displaystyle+\frac{a_{11}^{(0)}}{4}\left[1+4(\mathbf{c}\cdot\mathbf{w})^{2}-8|\mathbf{c}\cdot\mathbf{w}|\frac{K_{0}(2|\mathbf{c}\cdot\mathbf{w}|)}{K_{1}(2|\mathbf{c}\cdot\mathbf{w}|)}\right]$
		$\displaystyle-\frac{a_{00}^{(1)}}{2}\left[1-2(\mathbf{c}\cdot\mathbf{w})^{2}+|\mathbf{c}\cdot\mathbf{w}|\frac{K_{0}(2|\mathbf{c}\cdot\mathbf{w}|)}{K_{1}(2|\mathbf{c}\cdot\mathbf{w}|)}\right],$		(39c)
	$$\frac{\phi_{\vartheta}(\cos^{2}\vartheta)}{\phi_{\vartheta,M}(\cos^{2}\vartheta)}=1+\frac{9a_{00}^{(1)}}{8}\left(3\cos^{2}\vartheta-1\right).$$			(39d)

This yields $\langle c^{4}\rangle=\frac{15}{4}(1+a_{20}^{(0)})$, $\langle w^{4}\rangle=\frac{15}{4}(1+a_{02}^{(0)})$, $\langle c^{2}w^{2}\rangle=\frac{9}{4}(1+a_{11}^{(0)})$, $\langle(\mathbf{c}\cdot\mathbf{w})^{2}\rangle=\frac{3}{4}(1+a_{11}^{(0)}+\frac{5}{2}a_{00}^{(1)})$, and $\langle\cos^{2}\vartheta\rangle=\frac{1}{3}(1+\frac{9}{10}a_{00}^{(1)})$, as expected. The latter equality implies $b=a_{00}^{(1)}$, in consistency with the neglect of $a_{jk}^{(\ell)}$ with $j+k+2\ell\geq 3$ in Eq. (30). If, additionally, terms nonlinear in the cumulants are neglected in Eq. (29), we have

$${h\simeq b\simeq a_{00}^{(1)}}$$

(40)

in the Sonine approximation (38). Thus, if $a_{00}^{(1)}<0$ (as will be seen to be the case for most values of $\alpha$ and $\beta$ in the case of uniform spheres), we can expect that $\langle\cos^{2}\vartheta\rangle<\frac{1}{3}$. This implies a tendency of the vectors $\mathbf{c}$ and $\mathbf{w}$ to favor quasinormal mutual orientations (“lifted-tennis-ball” effect).