Kinetic Theory and Memory Effects of Homogeneous Inelastic Granular Gases under Nonlinear Drag

We study a dilute granular gas immersed in a thermal bath made of smaller particles with masses not much smaller than the granular ones in this work. Granular particles are assumed to have inelastic and hard interactions, losing energy in collisions as accounted by a constant coefficient of normal restitution. The interaction with the thermal bath is modeled by a nonlinear drag force plus a white-noise stochastic force. The kinetic theory for this system is described by an Enskog–Fokker–Planck equation for the one-particle velocity distribution function. To get explicit results of the temperature aging and steady states, Maxwellian and first Sonine approximations are developed. The latter takes into account the coupling of the excess kurtosis with the temperature. Theoretical predictions are compared with direct simulation Monte Carlo and event-driven molecular dynamics simulations. While good results for the granular temperature are obtained from the Maxwellian approximation, a much better agreement, especially as inelasticity and drag nonlinearity increase, is found when using the first Sonine approximation. The latter approximation is, additionally, crucial to account for memory effects such as Mpemba and Kovacs-like ones.


Introduction
Since the late 20th century, the study of granular materials has become of great importance in different branches of science, such as physics, engineering, chemistry, and mathematics, motivated by either fundamental or industrial reasons. It is well known that rapid flows in granular gases in the dilute regime are well described by a modified version of the classical Boltzmann's kinetic theory for hard particles. The most widely used model for the granular particles is the inelastic hard-sphere (IHS) one, in which particles are assumed to be hard spheres (or, generally, hard d-spheres) that lose energy due to inelasticity, as parameterized by a constant coefficient of normal restitution.
Theoretical predictions have been tested by different experimental setups in the freely evolving case [1,2]. However, it is rather difficult to experimentally replicate the latter granular gaseous systems due to the fast freezing implied by the dissipative interactions. Then, energy injection is very common in granular experiments [3][4][5][6][7][8][9][10]. In addition, granular systems are never found in a vacuum on Earth. From a quick but attentive glance at our close environment, grains might be found, for example, in the form of dust or pollen suspended in the air, sand, or dirtiness, diving down or browsing through a river, or even forming part of more complex systems such as soils. Therefore, fundamental knowledge about driven granular flows contributes to the understanding of a great variety of phenomena in nature. This is one of the reasons why the study of driven granular flows has become quite important, besides its intrinsic interest at physical and mathematical levels. Consequently, modeling driven granular flows constitutes a solid part of granular matter

The Model
We consider a homogeneous, monodisperse, and dilute granular gas of identical inelastic hard d-spheres of mass m and diameter σ, immersed in a background fluid made of smaller particles. In a coarse-grained description, the interactions between the grains and the fluid particles can be effectively modeled by a drag force plus a stochastic force acting on the grains. If the mass ratio between the fluid and granular particles is not very small, the drag force becomes a nonlinear function of the velocity [24][25][26]. The model, as said in Section 1, has previously been studied in the case of elastic collisions [21][22][23] but not, to our knowledge, in the context of the IHS model. Figure 1 shows an illustration of the system and its modeling. Figure 1. Illustration of the system considered in this paper. A granular gas of hard particles (represented by large yellowish spheres) is coupled to a thermal bath (made of particles represented by the small grayish spheres) via a drag force F drag = −mξ(v)v, where ζ(v) is a velocity-dependent drag coefficient, and a stochastic force F noise = mχ(v)η, where η is a Gaussian white-noise term. In addition, the granular particles are subjected to binary inelastic collisions, represented by the red gleam-like lines.

Enskog-Fokker-Planck Equation
The full dynamics of the system can be studied from the inelastic homogeneous Enskog-Fokker-Planck equation (EFPE), where f is the one-particle VDF, so that n = dv f (v; t) is the number density, and J[v| f , f ] is the usual Enskog-Boltzmann collision operator defined by Here, α is the coefficient of normal restitution (see below), v 12 = v 1 − v 2 is the relative velocity, σ = (r 1 − r 2 )/σ is the intercenter unit vector at contact, g c = lim r→σ + g(r) is the contact value of the pair correlation function g(r), + d σ ≡ d σ Θ(v 12 · σ), Θ being the Heaviside step-function and v i refers to the precollisional velocity of the particle i. Within the IHS model, the collisional rules are expressed by [18,30] From Equation (3), one gets (v 12 · σ) = −α(v 12 · σ); this relation defines the coefficient of normal restitution, which is assumed to be constant. The second term on the left-hand side of Equation (1) represents the action of a net force F = F drag + F noise describing the interaction with the particles of the background fluid. The deterministic nonlinear drag force is F drag = −mξ(v)v, where the drag coefficient ξ(v) depends on the velocity. In turn, F noise = mχ 2 (v)η is a stochastic force, where χ 2 (v) measures its intensity, and η is a stochastic vector with the properties of a zero-mean Gaussian white noise with a unit covariance matrix, i.e., where i and j are particle indices, and I is the d × d unit matrix so that different Cartesian components of η i (t) are uncorrelated. The functions ξ(v) and χ 2 (v) are constrained to follow the fluctuation-dissipation theorem as /m being the thermal velocity associated with the background temperature T b . The drag coefficient ξ is commonly assumed to be independent of the velocity. However, a dependence on v cannot be ignored if the mass of a fluid particle is not much smaller than that of grain [24][25][26]. The first correction to ξ = const is a quadratic term [21][22][23], namely where ξ 0 is the drag coefficient in the zero-velocity limit and γ controls the degree of nonlinearity of the drag force.

Dynamics
It is well known that, in the case of driven granular gases [11,12,14,[17][18][19]31,32], there exists a competition between the loss and gain of energy due to inelasticity and the action of the thermal bath, respectively. This eventually leads the granular gas to a steady state, in contrast to the freely cooling case [18].
The basic macroscopic quantity characterizing the time evolution of the system is the granular temperature, defined analogously to the standard temperature in kinetic theory as While in the case of elastic collisions, the asymptotic steady state is that of equilibrium at temperature T b , i.e., lim t→∞ T(t) = T b , in the IHS model, the steady state is a nonequilibrium one and, moreover, lim t→∞ T(t) = T st < T b . From the EFPE, one can derive the evolution equation of the granular temperature, which is given by where is the cooling rate and is the excess kurtosis (or fourth cumulant) of the time-dependent VDF. The coupling of T(t) to a 2 (t) is a direct consequence of the quadratic term in the drag coefficient. As for the cooling rate ζ(t), it is a consequence of inelasticity and, therefore, vanishes in the elastic case (conservation of energy). Insertion of Equation (2) into Equation (9) yields [18] Here, v th (t) = 2T(t)/m is the time-dependent thermal velocity and is the time-dependent collision frequency. Let us rewrite Equation (8) in dimensionless form. First, we introduce the reduced quantities (13) where ν b = g c K d nσ d−1 v b is the collision frequency associated with the background temperature T b . Note that the control parameter ξ * 0 measures the ratio between the characteristic times associated with collisions and drag. In the molecular case, ξ * 0 depends on the bath-to-grain density, size, and mass ratios, but otherwise, it is independent of T b [21,26]. In terms of the quantities defined in Equation (13), Equation (8) becomeṡ where henceforth, a dot over a quantity denotes a derivative with respect to t * , and we have taken into account that ζ(t)/ν(t) = 2µ 2 (t * )/d and ν(t)/ν b = θ 1/2 (t * ). Equation (14) is not a closed equation since it is coupled to the full VDF through a 2 and µ 2 . More generally, taking velocity moments on the EFPE, an infinite hierarchy of moment equations can be derived. In dimensionless form, it readṡ where (1 + a 2 ), and M 6 = d(d+2)(d+4)

8
(1 + 3a 2 − a 3 ), a 3 being the sixth cumulant. Equation (15) is trivial for = 0 and = 2. The choice = 4 yieldṡ Equations (14)- (16) are formally exact in the context of the EFPE, Equation (1). Nevertheless, they cannot be solved because of the infinite nature of the hierarchy (15) and the highly nonlinear dependence of the collisional moments µ on the velocity moments of the VDF. This forces us to devise tractable approximations in order to extract information about the dynamics and steady state of the system.

Maxwellian Approximation
The simplest approximation consists of assuming that the VDF remains very close to a Maxwellian during its time evolution so that the excess kurtosis a 2 can be neglected in Equation (14), and the reduced cooling rate µ 2 can be approximated by [11,12,17,18,28,33,34] In this Maxwellian approximation (MA), Equation (14) becomeṡ This is a closed equation for the temperature ratio θ(t * ) that can be solved numerically for any initial temperature. The steady-state value θ st in the MA is obtained by equating to zero the right-hand side of Equation (18), which results in a fourth-degree algebraic equation.

First Sonine Approximation
As we will see later, the MA given by Equation (18) provides a simple and, in general, rather accurate estimate of θ(t * ) and θ st . However, since the evolution of temperature is governed by its initial value only, the MA is unable to capture memory phenomena, such as Mpemba-or Kovacs-like effects, which are observed even in the case of elastic particles [21][22][23]. This is a consequence of the absence of any coupling of θ with some other dynamical variable(s).
The next simplest approximation beyond the MA consists of incorporating a 2 into the description but assuming it is small enough as to neglect nonlinear terms involving this quantity, as well as higher-order cumulants, i.e., a k 2 → 0 for k ≥ 2 and a → 0 for ≥ 3. This represents the so-called first Sonine approximation (FSA), according to which Equations (14) and (16) becomė where we have used [11,12,17,18,28,33,34] with Equations (19) make a set of two coupled differential equations. In contrast to the MA, now the evolution of θ(t * ) is governed by the initial values of both θ and a 2 . This latter fact implies that the evolution of temperature depends on the initial preparation of the whole VDF, this being a determinant condition for the emergence of memory effects, which will be explored later in Section 4.1.

Steady-State Values
The steady-state values θ st and a st 2 in the FSA are obtained by equating to zero the right-hand sides of Equations (19), i.e., Eliminating a st 2 in Equation (22) Figure 2 compares the MA and FSA predictions of θ st for three-and two-dimensional granular gases with ξ * 0 = 1. We observe that the breakdown of equipartition (as measured by 1 − θ st ) is stronger in 2D than 3D and increases with increasing inelasticity but decreases as the nonlinearity of the drag force grows. Apart from that, the deviations of the MA values with respect to the FSA ones increase with increasing nonlinearity and inelasticity, the MA values tending to be larger (i.e., closer to equipartition) than the FSA ones. The FSA predictions of a st 2 are displayed in Figure 3. First, it is quite apparent that the departure from the Maxwellian VDF (as measured by the magnitude of a st 2 ) is higher in 2D than 3D. It is also noteworthy that a st 2 starts growing with increasing γ, reaches a maximum at a certain value γ = γ max (α, ξ * 0 ), and then it decreases as γ increases beyond γ max (α, ξ * 0 ); this effect is more pronounced for small α.  Another interesting feature is that a st 2 takes negative values (in the domain of small inelasticity) only if γ is smaller than a certain value γ c . Of course, a st 2 (α, γ) α=1 = 0 for any γ (since the steady state with α = 1 is that of equilibrium), but ∂ α a st Thus, the critical value γ c is determined by the condition ∂ α a st 2 (α, γ c ) α=1 = 0. Interestingly, the result obtained from the FSA, Equation (24), is quite simple, namely which is independent of ξ * 0 .

Special Limits Absence of Drag
Let us first define a noise temperature T n as T n = T b ξ * 2/3 0 ∝ (ξ 0 T b ) 2/3 , so that θ 3/2 /ξ * 0 = (T/T n ) 3/2 . Now we take the limit of zero drag, ξ 0 → 0, with finite noise temperature T n . This implies T b → ∞, and thus, the natural temperature scale of the problem is no longer T b but T n , i.e., θ st → 0 but T st /T n = finite. From Equations (23) we see that F 0 (0) = d, F 1 (0) = 0, G 0 (0) = 0, and G 1 (0) = −d. Therefore, Equations (22) reduce tȯ where, for the sake of generality, we have undone the linearizations with respect to a st 2 . By the elimination of T n /T st 3/2 , one simply gets (d + 2)µ st 2 = µ st 4 , from which one can then obtain a st 2 upon linearization [11,12]. The steady-state temperature is given by T st /T n = (d/µ st 2 ) 2/3 .

Homogeneous Cooling State
If, in addition to ξ 0 → 0, we take the limit T n → 0, the asymptotic state becomes the homogeneous cooling state. In that case, T does not reach a true stationary value, but a 2 does. As a consequence, Equation (26a) is not applicable, but Equation (26b), with T n = 0, can still be used to get (d + 2)µ st 2 (1 + a st 2 ) = µ st 4 , as expected [11,12,17].

Collisionless Gas
If the collision frequency ν b is much smaller than the zero-velocity drag coefficient ξ 0 , the granular dynamics is dominated by the interaction with the background fluid and the grain-grain collisions can be neglected; therefore, the grains behave as Brownian particles. In that case, the relevant dimensionless time is no longer t * = ν b t but τ = ξ 0 t = ξ * 0 t * and the evolution equations (19) become It is straightforward to check that the steady-state solution is θ st = 1 and a st 2 = 0, regardless of the value of γ, as expected.

Comparison with Computer Simulations
We have carried out DSMC and EDMD computer simulations to validate the theoretical predictions. The DSMC method is based on the acceptance-rejection Monte Carlo Metropolis decision method [35] but adapted to solve the Enskog-Boltzmann equation [36,37], and the algorithm is, consequently, adjusted to agree with the inelastic collisional model [12,17] and reflect the interaction with the bath [23]. On the other hand, the EDMD algorithm is based on the one exposed in Ref. [23], but is adequated to the IHS collisional model. The main difference between DSMC and EDMD is that the latter does not follow any statistical rule to solve the Boltzmann equation but solves the equations of motion of the hard particles. Simulation details about the characteristics of the schemes and numerical particularities can be found in Appendix A.
In Figure 4, results from simulations are compared with the theoretical predictions of θ st (from MA and FSA) and of a st 2 (from FSA) in a three-dimensional (d = 3) IHS system with ξ * 0 = 1. It can be observed that both the DSMC and EDMD results agree with each other. From Figure 4a, one can conclude that, as expected, FSA works in the prediction of θ st much better than MA for values of γ close to γ max (α, ξ * 0 ) (which corresponds to the maximum magnitude of a st 2 ). Moreover, FSA gives reasonably good estimates for the values of a st 2 , although they get worse for increasing inelasticity, i.e., decreasing α. One might also think that the increase in γ produces a poorer approach; however, according to the theory, the performance of FSA improves if γ > γ max (α, ξ * 0 ), which corresponds to a decrease in |a st 2 |. Of course, nonlinear terms or higher-order cumulants might play a role that is not accounted for within FSA. Apart from the steady-state values, we have studied the temporal evolution of θ and a 2 , starting from a Maxwellian VDF at temperature T b , i.e., θ(0) ≡ θ 0 = 1 and a 2 (0) ≡ a 0 2 = 0. Note that this state is that of equilibrium in the case of elastic collisions (α = 1), regardless of the value of the nonlinearity parameter γ. The theoretical and simulation results are displayed in Figure 5 (a-d)) and FSA predictions, respectively. All states are initially prepared with a Maxwellian VDF at the bath temperature, i.e., θ 0 = 1 and a 0 2 = 0.
We observe that the relaxation of θ is accurately predicted by MA, except for the later stage with small α and/or large γ, in accordance with the discussion of Figure 4. This is remedied by FSA, which exhibits an excellent agreement with simulation results in the case of θ and a fair agreement in the case of a 2 , again in accordance with the discussion of Figure 4. It is also worth mentioning the good mutual agreement between DSMC and EDMD data, even though fluctuations are much higher in a 2 than in θ because of the rather small values of |a 2 |.

Memory Effects
Whereas the temperature relaxation from Maxwellian initial states is generally accurate from MA, it misses the explicit dependence of the temperature evolution on the fourth cumulant (see Equation (14)), which, however, is captured by FSA (see Equation (19a)). This coupling of θ to a 2 is a signal of preparation dependence of the system, hence, a signal of memory effects, as occurs in the elastic case reported in Refs. [21][22][23].

Mpemba Effect
We start the study of memory effects with the Mpemba effect [38][39][40][41][42]. This counterintuitive phenomenon refers to situations in which an initially hotter sample (A) of a fluid-or, more generally, a statistical-mechanical system-cools down sooner than an initially colder one (B) in a cooling experiment. We will refer to this as the direct Mepmba effect (DME). Analogously, the inverse Mpemba effect (IME) occurs in heating experiments if the initially colder sample (B) heats up more rapidly than the initially hotter one (A) [21,23,40,41,43]. In the special case of a molecular gas (i.e., α = 1), an extensive study of both DME and IME has recently been carried out [21,23]. Figure 6a,b present an example of DME and IME, respectively. As expected, FSA describes the evolution and crossing for temperatures of samples A and B very well. On the contrary, MA does not predict this memory effect. In addition, from Figure 6c,d we can conclude that FSA captures the relaxation of a 2 toward a st 2 = 0 quite well.

Kovacs Effect
Next, we turn to another interesting memory effect: the Kovacs effect [44,45]. In contrast to the Mpemba effect, the Kovacs effect has a well-defined two-stage protocol and does not involve a comparison between two samples. In the context of our system, the protocol proceeds as follows. First, the granular gas is put in contact with a bath at temperature T b1 and initialized at a temperature T 0 > T st 1 , T st 1 = θ st T b1 being the corresponding steady-state temperature (note that θ st is independent of T b1 at fixed ξ * 0 ). The system is allowed to relax to the steady state during a time window 0 < t < t K , but then, at t = t K , the bath temperature is suddenly modified to a new value T b , such that T(t K ) = T st , T st = θ st T b being the new steady-state value. If the system did not retain a memory of its previous history, one would have T(t) = T st for t > t K , and this is, in fact, the result given by the MA. However, the temperature exhibits a hump for t > t K , before relaxing to T st . This hump is a consequence of the dependence of ∂ t T on the additional variables of the system. According to Equation (14), and maintained in the FSA, Equation (19a), the first relevant quantity to be responsible for a possible hump is the excess kurtosis of the VDF, as occurs in the elastic limit [22]. In fact, at time t * = t * K , such that θ(t * K ) = θ st , the slope of the temperature according to FSA, Equation (19a), readṡ Thus, a nonzero difference a st 2 − a 2 (t * K ) implies the existence of a Kovacs-like hump, its sign being determined by that of this difference; that is, we will obtain an upward hump if a 2 (t * K ) < a st 2 or a downward hump if a 2 (t * K ) > a st 2 . For simplicity, in our study of the Kovacs-like effect, we replace the first stage of the protocol (0 < t * < t * K ) by just generating the state at t * = t * K with θ(t * k ) = θ st and a 2 (t * K ) = a st 2 (see Appendix A). The effect is illustrated in Figure 7 for the same system as in Figure 6 with the choices a 2 (t * K ) = −0.35 < a st 2 and a 2 (t * K ) = 0.4 > a st 2 . Again, the DSMC and EDMD results agree with each other and with the theoretical predictions. However, in the case a 2 (t * K ) = −0.35 (upward hump), Figure 7a, we observe that the theoretical curve lies below the simulation results. This might be caused by a nonnegligible value of the sixth cumulant a 3 (t * K ) = −0.375, as reported in Ref. [23] in the elastic case. Apart from this small discrepancy, FSA captures the magnitude and sign of the humps, as well as the relaxation of the fourth cumulant, very well.

Conclusions
In this work, we have looked into the dynamics of a dilute granular gas immersed in a thermal bath (at temperature T b ) made of smaller particles but with masses comparable to those of the grains. To mathematically characterize this system, we have worked under the assumptions of Boltzmann's kinetic theory, describing the system by the one-particle VDF, whose evolution is monitored by the EFPE, Equation (1), for the IHS model of hard d-spheres. The action of the bath on the dynamics of the granular gas is modeled by a nonlinear drag force and an associated stochastic force. At a given dimensionality d, the control parameters of the problem are the coefficient of normal restitution (α), the (reduced) drag coefficient at zero velocity (ξ * 0 ), and the nonlinearity parameter (γ). After a general presentation of the kinetic theory description in Section 2, we obtained the evolution equation of the reduced temperature θ(t * ) ≡ T(t)/T b (Equation (14)), which is coupled explicitly with the excess kurtosis, a 2 , and depends on every velocity moment through the second collisional moment µ 2 (which is nonzero due to inelasticity). Therefore, the whole dynamics in the context of the EFPE is formally described by Equation (14) and the infinite hierarchy of moment equations given by Equation (15). In order to give predictions, we proposed two approximations. The first one is MA, which consists of assuming a Maxwellian form for the one-particle VDF, whereas the second one, FSA consists of truncating the Sonine expansion of the VDF up to the first nontrivial cumulant a 2 . Their evolution equations are given by Equations (18) and (19), respectively. The predictions for the steady-state values are exposed in Figures 2 and 3, which show some small discrepancies in θ st between MA and FSA as we increase the inelasticity (decreasing α). Moreover, we observed that, for fixed α and ξ * 0 , a st 2 gets its maximum value when the nonlinearity parameter is γ = γ max (α, ξ * 0 ). Another interesting feature is the existence of a critical value γ c , such that for γ > γ c , the values of a st 2 are always positive for every value of α, while for γ < γ c , we find a st 2 < 0 for inelasticities small enough. Interestingly, the value of γ c given by Equation (25) is found to be independent of ξ * 0 . In addition, some already known limits are recovered in Section 3.2.2.
Furthermore, in order to check the predictions from MA and FSA equations, we carried out DSMC and EDMD simulations for hard spheres (d = 3) with fixed ξ * 0 = 1 (which corresponds to comparable time scales associated with drag and collisions). First, from Figure 4a, we can conclude that, whereas MA provides good predictions of θ st , except for large inelasticities and values of γ close to γ max , FSA is much more accurate because it takes into account the influence of a st 2 . The latter approach is generally reliable for a st 2 , as observed in Figure 4b, although, not unexpectedly, it slightly worsens as |a st 2 | grows. Relaxation curves starting from a Maxwellian initial state in Figure 5 show that FSA agrees very well with both DSMC and EDMD; however, MA exhibits good agreement during the first stage of the evolution but becomes less reliable as the steady state is approached.
A relevant feature of these systems, as already studied in the elastic case [21][22][23], is the emergence of memory effects, which are not contemplated by MA. FSA predicts the emergence of the Mpemba effect very well for both DME and IME, as can be seen in Figure 6. Analogously, Kovacs-like humps, both upward and downward, are correctly described by FSA, as observed in Figure 7, although the FSA humps are slightly less pronounced (especially the upward one) than the simulation ones. This is presumably due to the role played by a 3 and higher-order cumulants, as occurs in the elastic limit reported in Ref. [23].
To conclude, we expect that this work will motivate research about this type of system and the emergence of memory effects. For instance, one can extend the study to other collisional models (such as that of rough spheres), to nonhomogeneous states, or to a more detailed description of the memory effects observed. Institutional Review Board Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available in the online repository: https://github.com/amegiasf/GranularNonlinearDrag.

Acknowledgments:
The authors are grateful to the computing facilities of the Instituto de Computación Científica Avanzada of the University of Extremadura (ICCAEx), where the simulations were run.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: In the free-streaming stage, each particle velocity is updated according to an Euler numerical algorithm of a Langevin-like equation derived from an Itô interpretation of the Fokker-Planck part of the EFPE (see Ref. [23] where Y i is a random vector drawn from a Gaussian probability distribution with unit variance, P(Y) = (2π) −d/2 e −Y 2 /2 . In the implementations of the DSMC algorithm, we used N = 10 4 hard spheres (d = 3) and a time step ∆t = 10 −2 λ/v b , λ = ( √ 2πnσ 2 ) −1 being the mean free path.

Appendix A.2. Event-Driven Molecular Dynamics
EDMD methods compute the evolution of particles driven by events: particle-particle collisions, boundary effects, or other more complex interactions. Analogously to the splitting described in the DSMC algorithm, free streaming of particles occurs between two consecutive events. Here, we need to consider the influence of the stochastic and drag forces not only in the velocities but also in the positions of the N granular particles, {r i } N i=1 . In order to account for this, we followed the approximate Green Function algorithm proposed in Ref. [46]. Whereas the velocities are updated according to Equation (A1), the positions follow where W i = Y i + √ 5/3Y i , Y i being another random vector drawn from P(Y) = (2π) −d/2 e −Y 2 /2 .