# Hydrodynamics of a Granular Gas in a Heterogeneous Environment

## Abstract

We analyze the transport properties of a low density ensemble of identical macroscopic particles immersed in an active fluid. The particles are modeled as inelastic hard spheres (granular gas). The non-homogeneous active fluid is modeled by means of a non-uniform stochastic thermostat. The theoretical results are validated with a numerical solution of the corresponding the kinetic equation (direct simulation Monte Carlo method). We show that a steady flow in the system that is accurately described by Navier-Stokes (NS) hydrodynamics, even for high inelasticity. Surprisingly, we find that the deviations from NS hydrodynamics for this flow are stronger as the inelasticity decreases. The active fluid action is modeled here with a non-uniform fluctuating volume force. This is a relevant result given that hydrodynamics of particles in complex environments, such as biological crowded environments, is still a question under intense debate.

1 \articlenumberx \doinum10.3390/—— \pubvolumexx 2017 \copyrightyear2017 \externaleditorAcademic Editor: name \historyReceived: date; Accepted: date; Published: date \TitleHydrodynamics of a Granular Gas in a Heterogeneous Environment \AuthorFrancisco Vega Reyes and Antonio Lasanta \AuthorNamesFrancisco Vega Reyes and Antonio Lasanta \corresCorrespondence: alasanta@ing.uc3m.es

## 1 Introduction

Biological Systems, and more specifically, systems with active particles, are in general out of equilibrium, which implies that for these systems the results of equilibrium statistical mechanics do not in general apply. The term active matter refers to a system with an ensemble of particles that are able to transform energy - internal or environmental- into movement [Marchetti and et al.(2013)]. A correct description of transport phenomena in this type of systems can be achieved via a specific kinetic theory, that takes into account the peculiar energy processing of these particles, see for example [Bertin et al.(2009)Bertin, Droz, and Grégoire, Bertin et al.(2006)Bertin, Droz, and Grégoire, Chou et al.(2012)Chou, Wolfe, and Ihle, Ihle(2016)] for the well known Vicsek-like models [Vicsek et al.(1995)Vicsek, CzirÃ³k, Ben-Jacob, Cohen, and O.Shochet].

In recent years, the study of active particles in crowded environments has become the subject of intense study [Bechinger et al.(2016)Bechinger, Leonardo, Lowen, Reichhardt, and Volpe]. One of the motivations is the large number of promising applications derived from the use of active particles as micro-machines [Bechinger et al.(2016)Bechinger, Leonardo, Lowen, Reichhardt, and Volpe]. In addition, the behaviour of inert particles inside complex active environments remains an open question. For instance, there are experimental studies on inert Brownian particles trapped by a harmonic potential in a fluid composed by bacteria [Argun et al.(2016)Argun, Moradi, PinÃ§e, Bagci, Imparato, and Volpe]. Also, Sándor et al. [Sándor et al.(2017a)Sándor, LibÃ l, Reichhardt, and Reichhardt] study active Brownian particles on a travelling-wave substrate. Other works study particles which can be trapped by holes [Sándor et al.(2017b)Sándor, LibÃ l, Reichhardt, and Reichhardt, Sándor et al.(2017c)Sándor, LibÃ l, Reichhardt, and Reichhardt]. More realistic active systems such as biological tissues, share important analogies with granular matter, namely glass transition, cage effect, fluidization, solidification and the appearance of giant density fluctuations [Angelini et al.(2011)Angelini, Hanezzo, Trepat, Marquez, Fredberg, and Weitz, Malinverno et al.(2017)Malinverno, Corallino, and et al., Garcia et al.(2015)Garcia, Hannezo, Elgeti, Joanny, Silberzan, and Gov, Bi et al.(2015)Bi, Lopez, Schwarz, and Manning]. There, active complex non-homogeneous forces appear. Also the non-Gaussian behaviour of diffusing non-active particles in heterogeneous media has been studied recently for a spatially varying diffusion coefficient [Malgaretti et al.(2016)Malgaretti, Pagonabarraga, and Rubi].

On the other hand, the transport properties of granular matter are of interest at a fundamental level and also for applications, since granular matter is present in a variety of industry and technology sectors [Dufty(2001), Brilliantov and Pöschel(2004), Goldhirsch(2003), de Bruyn(2011)]. Granular transport theories usually make a connection with traditional continuum mechanics theories, due to fundamental similarities observed with classical fluids [Dufty(2001), Goldhirsch(2003)]. Granular particles loose a fraction of their kinetic energy after collision and for this reason we need to continuously excite the system, in order to sustain the dynamics. This excitation usually comes from the system boundaries, which results in inhomogeneous flows that may eventually become steady [Batchelor(1953), Dufty(2001), Goldhirsch(2003)].

We propose in this work a study that tries to model a system of macroscopic inert particles that are immersed in an active fluid with inhomogeneous temperature. For this, we model the system as a granular fluid [Haff(1983)] subject to a nonuniform stochastic volume force (in the form of a white noise [Puglisi et al.(2012)Puglisi, Gnoli, Gradenigo, Sarracino, and Villamaina]). The inelastic cooling term, inherent in a granular fluid [Goldhirsch(2003)], and a stochastic force allow us in this work to model the energy sink and source terms that are characteristic to active particle systems [Marchetti and et al.(2013), Bechinger et al.(2016)Bechinger, Leonardo, Lowen, Reichhardt, and Volpe]. More specifically, we focus on a very low density granular system (i.e., a granular gas) surrounded by a low density (inhomogeneous) fluid that can be considered as a coarse grained version of an active matter system. Since the system has very low density, particle velocity correlations can be neglected and thus the the inelastic Boltzmann equation applies [Goldhirsch(2003), Puglisi(2014)]). Moreover, systems of solid particles suspended in a low density interstitial fluid usually fall into the grain-grain collision dominated regime [Jaeger and Nagel(1992)], as opposed to the interstitial fluid viscosity-dominated regime, that applies for most cases of interstitial liquids [Geminard et al.(1999)Geminard, Losert, and Gollub].

In previous theoretical works, granular gas fluidization has been studied only in the simplistic case of a homogeneous interstitial fluid [Williams and MacKintosh(1996), Cafiero et al.(2000)Cafiero, Luding, and Herrmann, Visco et al.(2006)Visco, Puglisi, Barrat, Trizac, and van Wijland, Villamaina et al.(2008)Villamaina, Puglisi, and Vulpiani, Fiege et al.(2009)Fiege, Aspelmeier, and Zippelius, Gradenigo et al.(2011)Gradenigo, Sarracino, Villamaina, and Puglisi]. And yet this theoretical homogeneous state is not quite a realistic situation [Losert et al.(2000)Losert, Bocquet, Lubensky, and Gollub]. Neat examples of inhomogeneous granular gas suspensions are the large sand plumes that can be observed in the Earth’s atmosphere for weeks [Yu and et al.(2015)]. A homogeneous energy source is obviously not a realistic situation either in active fluids, at a biological level [Marchetti and et al.(2013)]. For this reason we focus in this work on the more elaborate case of a non-uniform interstitial fluid.

However, the intent of this work is not experimental validation but to extract a more clear picture of the relevant physics of this kind of systems. As a result, we show there are granular gas flows immersed in an inhomogeneous active fluid hat can be accurately described with Navier-Stokes (NS) hydrodynamics, even at high inelasticities. Furthermore, we surprisingly found cases where NS hydrodynamics is a better approximation for the granular gas than for a perfectly elastic gas subjected to the same boundary forcing conditions. In particular, we will show that, once the granular gas is fluidized, the intensity of the spatial gradients in response to external excitations (those coming from the boundaries) can actually be weaker for granular gases than for elastic gases.

In a previous experimental work [Ojha et al.(2004)Ojha, Lemieux, Dixon, Liu, and Durian], local balance between the total energy input (fluctuating force plus boundary heating) and the energy sink was found. We use now this kind of balance condition. As we will show below, this local balance results, in the specific set-up considered in the present work, in a steady flow with uniform heat flux throughout the system. This balance condition is in fact analogous to the one occurring in well known non-Newtonian granular flows like the uniform shear flow [Campbell(1990), Vega Reyes et al.(2010)Vega Reyes, Santos, and Garzó].

The structure of the paper is as follows: in Section 2 we present a description of the system and the corresponding steady state equations. In section 3 we describe the flows where there is a balance between inelastic cooling coming from particle collisions and volume energy input (and with no shear). Section 4 is devoted to a brief discussion on the steady flows resulting from adding a weak shear Section 5 presents a discussion on the transport coefficients and rheology properties of the studied flows and the paper conclusions are drawn in Section 6. The simulation method and transport coefficient equations are discussed in the Appendices .1 and .2 respectively.

## 2 System and Steady base state equations

### 2.1 Description of the system

Our system consists of a set of identical inert smooth hard spheres (or disks) immersed in an active fluid. The system is limited by two infinite parallel hard walls, located at planes respectively. The walls act like two distinct kinetic energy sources, characterized with temperatures . They may have relative movement (wall velocities respectively), eventually inducing the particles to continuously flow along the channel between them.

Collisions between the inert particles do not preserve energy (i.e., particle collisions are inelastic). The particles have a diameter (radius) that we denote as and their mass is . For inelastic smooth hard spheres/disks, inelasticity may be characterized accurately, in a range of experimental situations, by a constant coefficient of normal restitution [Foerster et al.(1994)Foerster, Louge, Chang, and Alla, Grasselli et al.(2015)Grasselli, Bossis, and Morini]. This coefficient ranges from 1 for purely elastic collisions to 0 for purely inelastic collisions. In addition, we model the interstitial active fluid injection of thermal energy onto the granular gas particles as a stochastic force . The equation of motion for a particle can be written

(1) |

where is the force due to inelastic collisions, is the force exerted by the heterogeneous medium and the mass of particles (set to 1 for simplicity).

We can model this interaction by means of a fluctuating volume force [Van Kampen(1992)], that we denote as and fulfills the conditions [Marini-Bettolo-Marconi et al.(2007)Marini-Bettolo-Marconi, Tarazona, and Cecconi, Gradenigo et al.(2011)Gradenigo, Sarracino, Villamaina, and Puglisi, Van Kampen(1992), Ojha et al.(2004)Ojha, Lemieux, Dixon, Liu, and Durian]

(2) |

being the fluctuating force intensity [Marini-Bettolo-Marconi et al.(2007)Marini-Bettolo-Marconi, Tarazona, and Cecconi, Gradenigo et al.(2011)Gradenigo, Sarracino, Villamaina, and Puglisi] (that has dimensions of velocity squared times time), the unit matrix in dimensional space, is the Kronecker delta function and indicates average of magnitude over a sufficiently long time interval (long compared with the characteristic microscopic time scale). It has been shown in the case of the large mass limit and homogeous energy injection, that equation 1 for a particle can be written as a Langevin equation corresponding to a granular brownian particle [Sarracino et al.(2010)Sarracino, Villamaina, Constantini, and Puglisi].

On the other hand, collisions between two grains follow the inelastic smooth hard sphere collisional model

(3) | |||

(4) |

with and . The primes denote pre-collisional velocities.

Notice that in (2) we have (a space-dependent noise intensity). The noisy term appears in the inelastic Boltzmann equation as a Fokker-Planck-like operator [Garzó et al.(2013)Garzó, Chamorro, and Vega Reyes, Marconi et al.(2011)Marconi, Puglisi, Rondonic, and Vulpiani]

(5) |

where is the collisional integral for inelastic hard spheres and whose expression may be found elsewhere [Goldhirsch(2003), Brey et al.(1998)Brey, Dufty, Kim, and Santos], is particle velocity, and its modulus.

Taking into account the system geometries, we consider only , that is the simplest space-dependence our geometry configuration allows for an inhomogeneous active fluid. The steady base states may be deduced from the reduction of the kinetic equation (5)

(6) |

We have numerically solved the kinetic equation (6) by means of the Direct Simulation Monte Carlo method [Brilliantov and Pöschel(2004), Bird(1994)], applying in each case the corresponding boundary conditions and heterogeneous fluctuating force properties. For more details on the simulation algorithm, see the corresponding section in the Appendix .1.

### 2.2 Steady base state equations

Multiplying by velocity momenta the kinetic equation (6) and performing velocity integrals, we may easily obtain from (6) the mass, momentum and energy balance equations [Vega Reyes and Urbach(2009), Garzó et al.(2013)Garzó, Chamorro, and Vega Reyes]

(7) | |||

(8) |

where is the flow velocity field (with its y-component), is the particle density, and and are momentum flux tensor (also called stress tensor) and heat flux vector respectively, and is the cooling rate, a magnitude that emerges from the energy balance due to inelasticity in the collisions [Goldhirsch(2003), Brilliantov and Pöschel(2004)]. Equations (7) and (8) contain no approximations; i. e. they are exact for this type of geometry, even if the steady state were not hydrodynamic. Notice also that the stress tensor element is homogeneous and is not a function of time. This condition breaks if the flow is not laminar, but we will only consider laminar flows in this work.

We take now into account the characteristics of NS hydrodynamics in (7) and (8). In (7) we use that all diagonal elements of the stress tensor are equal and we also use that the horizontal heat flux is null (they both raise from fluxes terms that are of second order in the gradients [Burnett(1935)]) and thus, from the first equation in (7), we obtain that . We also use the fluxes expressions at NS order [Brey et al.(1998)Brey, Dufty, Kim, and Santos]: , , where is the viscosity and are the heat flux transport coefficients. With this, the steady state forms of (7) and (8) are

(9) | |||

(10) |

Next, we take into account that the NS transport coefficients and cooling rate scale in the following forms with temperature [Brey et al.(1998)Brey, Dufty, Kim, and Santos]: , , , . Here, the coefficients , are related to the elastic gas [Chapman and Cowling(1970)] NS viscosity and thermal conductivity , and thus

(11) | |||

(12) |

where is the gamma function [Abramowitz and Stegun(1965)] and for disks and for spheres. Thus, the steady state balance equations (9), (10) yield

(13) | |||

(14) |

where , and we call it effective thermal conductivity and is the reference unit of temperature, that we choose as (temperature at the lower, and colder, wall). Here, is the dimensionless cooling rate [Brilliantov and Pöschel(2004), Vega Reyes et al.(2013)Vega Reyes, Santos, and Garzó] and the expressions of the reduced transport coefficients viscosity , thermal conductivity , and that we used are the ones obtained by Garzó and Montanero [Garzó and Montanero(2002)] for a granular gas heated by a white noise. Their expressions are displayed in the transport coefficients section in the Appendix .2.

We use a scaled space variable that allows us to find an analytical solution for the flow velocity (13) and the temperature (14) profiles. This variable is defined by the relation . When expressed in terms of this scaled variable, equations (9)-(10) read

(15) | |||

(16) |

where and are dimensionless magnitudes in the system and we call them local shear rate and temperature curvature, respectively. In (15), (16) and henceforth we use dimensionless variables indicated with a hat, where we have chosen as set of units: for temperature (already defined), for mass, for pressure (the steady state hydrostatic pressure, that as we said is a global constant for the system). As a unit for length we use

(17) |

that is the reference mean free path [Chapman and Cowling(1970)] and for time we use

(18) |

(inverse of reference collision frequency). These two definitions imply that the unit for velocity is . We indicate the spatial coordinate in dimensionless form as , and analogously for all dimensionless magnitudes.

The reader may infer the physical meaning of and with the aid of (7), (8): (reduced shear rate) measures how rapidly varies the flow velocity field, in units of the mean free path; measures the degree of imbalance between energy volume sinks and sources per unit time (as we explain above in Eq. (18), this unit being the inverse of the collision frequency). The differential equation (16) has been observed for the granular temperature in similar experimental configurations [Losert et al.(2000)Losert, Bocquet, Lubensky, and Gollub].

With the presence of the stochastic volume force in the energy balance equation (8), it is straightforward to deduce from (14)- (16) and the definitions above that has the form

(19) | |||||

that, as it can be noticed, depends in general on position through .

## 3 Steady Base States with energy balance and no shear

We now want to look up for steady states in a granular gas fluidized by an heterogeneous source that may be accurately described by hydrodynamics at NS order. We will denote them as ’reference states’. As in rarefied gases, a close to unity Knudsen number () signals hydrodynamics failure. Moreover, hydrodynamics at first order in the gradients (NS equations) are usually only valid if takes low values [Goldhirsch(2003), Bird(1994)]. The Knudsen number is defined as a representative microscopic time or length scale over a macroscopic time or length scale. However, the choice of reference scales is not unique and picking an appropriate Knudsen number for each specific flow problem is a subtle question [Boyd et al.(1995)Boyd, Chen, and Candler].

It was shown in a previous work [Vega Reyes and Urbach(2009)] that Couette-Fourier flows (i.e., those occurring in our system, under the appropriate circumstances) are characterized by three representative and independent microscopic vs. macroscopic relative scales (or Knudsen numbers) that happen to be constant throughout the system: , (look at eqs. (15) and (16)) and , with , where measures the size of the temperature gradient imposed by the boundaries, this gradient being referred to mean free path units. These numbers arise from the specific form of the relevant differential equations for steady states, and the boundary geometry. Since the geometry of the differential equations (15), (16) and boundary conditions are the same in this work, we may safely use the same set of reference Knudsen numbers (see previous work [Vega Reyes and Urbach(2009)] for more details on this issue). Furthermore, another previous analysis [Vega Reyes et al.(2013)Vega Reyes, Santos, and Garzó] for granular gases without an interstitial fluid showed that has a much weaker effect on the emergence non-Newtonian behavior.

Thus, let us analyze the class of solutions for steady flows where this is the only active source of spatial gradients; i.e., and

(20) | |||

(21) |

Moreover, in the limit this state corresponds to the classic Fourier flow of a molecular gas. In order to obtain (with ), and taking into account (19) we find

(22) |

Equation (20) implies that flow velocity is constant respect to , which, from the definition of leads to also constant (and to all physical effects, flow velocity may be considered as null). Moreover, equation (21) implies that temperature is linear in , which (since ) yields and thus ; i.e., the stochastic force generates a Fourier-like temperature profile [Vega Reyes et al.(2010)Vega Reyes, Santos, and Garzó]. At a physical level, this makes sense if we consider that the granular gas is fluidized by a low viscosity interstitial fluid that is heated from the boundaries. It is easy to deduce from (8), (10) and (19) that implies that inelastic cooling is balanced by volume energy input and thus, that heat flux is uniform throughout the system.

Therefore our base steady states are (not previously reported) granular flows with uniform heat flux. Summarizing, these NS steady states have the following properties: i) , , ii) constant hydrostatic pressure , iii) null normal stress differences () and null [Burnett(1935)], iv) constant , and v) linear profiles, or equivalently, [Vega Reyes et al.(2010)Vega Reyes, Santos, and Garzó]. Properties i) to iii) are fulfilled by all Couette-Fourier steady flows at NS order [Vega Reyes and Urbach(2009)] and conditions iv) and v) are fulfilled by all Couette-Fourier steady flows with uniform heat flux [Vega Reyes et al.(2010)Vega Reyes, Santos, and Garzó].

We need to remark that the relation in (22) between fluctuating force intensity and temperature is not a hypothesis but a consequence of the experimental fact that there is local energy balance. On the other hand, there is an essential difference of this new flow respect to previously studied granular flows with uniform heat flux (the simple shear flow included [Campbell(1989), Campbell(1990)]).

Very surprisingly, we can show that contrary to other steady granular flows, the kind of state defined by (22) can be accurately described in the context of Navier-Stokes equations. In order to prove this feature, we have simulated the system described by equations (20), (21) and (22) by means of the DSMC method (more details on simulations are given in the Appendix .1). We have checked that the temperature profiles obtained under steady state in effect follow the behaviour , inherent to gas steady flows with uniform heat flux [Vega Reyes et al.(2010)Vega Reyes, Santos, and Garzó]. Results can be seen in Figure 1, for which the agreement with the NS theory prediction is excellent, independently of the inelasticity value.

## 4 Weakly sheared steady states

We discuss here weakly sheared states, in order to capture the gradual the eventual appearance of non-Newtonian behavior. Thus, starting out of the reference base states in the previous subsection; i.e., still with condition (22) being fulfilled, we analyze now sheared states (), for which we have from (19),

(23) |

For this type of flow we have also performed DSMC simulations. Relation (23) helps us understand that an increase in shear rate (one of the relevant Knudsen numbers) implies an increase in (one of the other Knudsen numbers). Thus, for increasing shear we should gradually depart from NS hydrodynamics as both and increase. However, relation (23) is valid only as long as NS hydrodynamics applies. For non-Newtonian regime, the relation between the temperature curvature , local shear rate and coefficient of restitution would be more intricate and in general it must be solved numerically [Vega Reyes et al.(2013)Vega Reyes, Santos, and Garzó]. A non-linear theory is beyond the scope of the present work, but, fortunately, we may use results from DSMC to analyze deviations from NS hydrodynamics.

The results we obtained are quite surprising (see Table 1). In effect, we can see that similar values of local shear rate tend to yield greater values of the other Knudsen number for more elastic gases than for more inelastic gases. This means that, contrary to what would be expected, the Knudsen numbers in this type of flow tend globally to be smaller for higher inelasticities, which could indicate that departures from NS hydrodynamics are more important, for similar shear rates, for more elastic gases. But we will confirm this point in the following section 5, by measuring from DSMC simulations the set of transport coefficients and rheological properties and comparing with NS theory predictions.

0.99 | 0.173 | 0.058 | 0.115 | 0.036 |

0.9 | 0.191 | 0.034 | 0.120 | 0.018 |

0.8 | 0.193 | 0.027 | 0.119 | 0.017 |

0.7 | 0.193 | 0.024 | 0.119 | |

0.6 | 0.188 | 0.115 | ||

0.5 | 0.187 | 0.113 |

We also performed an additional series with varying (not shown in figures) and checked that the measured transport coefficients do not depend on its value, analogously to the result in a previous work on granular flows with uniform heat flux [Vega Reyes et al.(2010)Vega Reyes, Santos, and Garzó].

## 5 Transport coefficients and rheology

Results in the preceding section suggest that: a) our steady states are well described by NS hydrodynamics, even for strong inelasticities, as long as shearing from the boundaries is not too strong; b) deviations from NS hydrodynamics are stronger for more elastic gases. Both observations would be contrary to what occurs in steady granular flows (like the simple shear flow [Campbell(1990)], just to put one example).

Let us analyze in more detail to what extent the reference (no-shear) steady states described in section 3 and the sheared states described in section 4 are strictly NS flows or not, by analyzing the relevant transport properties. This analysis will be done as a function of the degree of inelasticity in grains collisions and intensity of shearing. In order to analyze non-Newtonian behavior, we need to define the cross thermal conductivity transport coefficient and the reduced normal stresses, defined by

(24) |

where is the thermal conductivity for a gas with elastic particle collisions [Chapman and Cowling(1970)]. For hydrodynamics at NS order, and . The degree of divergence from these numerical values quantifies eventual departures from NS hydrodynamics. The appearance of and occurs with the emergence in the fluxes of terms of second order in the gradients (Burnett order). For our geometry, the only remaining second order terms are of form , in the heat flux and of the form , , in the stress tensor (moment flux) respectively (the rest of the gradient second order terms in the fluxes are null for the geometry of our system) [Burnett(1935), Chapman and Cowling(1970)].

In Fig. 2 we present the theory-simulation comparison for reduced viscosity and effective thermal conductivity , for different shearing levels. For sufficiently low shearing, the agreement between theory and simulation is very good for the range of inelasticities represented here. In the limit of no shear, the agreement for the thermal conductivity is excellent from down to . The viscosity, however, cannot be measured on the reference (null shear) steady states. Of course, this limitation is not intrinsic of the granular gas since the same happens for classic fluids. Nevertheless, a relatively good viscosity measurement can be obtained in the series with lowest shear used in the simulations (), for which the agreement at high inelasticities between theory and simulation is actually slightly better than for the thermal conductivity measured at . For both coefficients, and at high inelasticities (), the agreement is improved when the zeroth order distribution function in the Chapman-Enskog method [Chapman and Cowling(1970)] is not the Maxwellian but an approximation to the solution of the homogeneous state for the heated granular gas (slightly different from the Maxwellian [Garzó et al.(2007)Garzó, Santos, and Montanero]). This approximation improve (represented with dashed lines in Fig. 2) consists in taking into account the first term in an expansion in Sonine polynomials around the Maxwellian (for more reference on the Sonine polynomials expansion method in kinetic theory, see for instance [Burnett(1935), Chapman and Cowling(1970), Truesdell(1973)]).

We report in Fig. 3 (a) the reduced normal stresses measured in the simulations, as a function of the coefficient of restitution and for several series with different shearing intensities. As we see, deviations from the NS behavior () are small (less than ) even at moderately high values of shear (close to ). In this case, the coefficient of restitution seems not to have an important impact on the behavior of both and . In Fig. 3 (b) we represent the value of the cross conductivity coefficient . As we can see, the deviations from linear behavior (i.e., from ) are more significant here: up to a maximum of approximately for the highest shear rate series represented here. However, as expected, the NS behavior is recovered for small shear rate (always less than for the lowest non-null shear rate series). It is remarkable that, again, departures from NS behavior are stronger in the quasi-elastic limit than for more inelastic gases (i.e., more inelastic flows show less significantly non-Newtonian behavior here). This is consistent with the non-trivial observation in the previous section that increases as we approach the quasi-elastic limit. As we already noticed, a higher implies higher Knudsen number, and thus, our observations are self-consistent and not an artifact of the theory. The NS behavior (, ) is strictly recovered for the reference states (), for all values of inelasticity.

## 6 Conclusion

In summary, our results on the transport properties of our system strongly indicate that NS hydrodynamics may apply for steady flows of inert particles in a heterogeneous medium, as long as the relevant Knudsen number scales imposed by the boundary conditions in the problem do keep small. This is in no way different to the condition that is required to a molecular gas in order to obey hydrodynamics at NS order [Boyd et al.(1995)Boyd, Chen, and Candler].

More generally, the fact that this result is independent of the degree of inelasticity supports the hypothesis that inelasticity by itself does not cause necessarily the breakdown of scale separation in granular gases such as it appears in the Homogeneous Cooling State. A recent work on the aging time to hydrodynamics for the homogeneous cooling state of the granular gas also supports this hypothesis [Vega Reyes et al.(2014)Vega Reyes, Santos, and Kremer]. The dynamics is actually a bit more complex and the energy source necessary to fluidize the system may cancel out the gradients inherently produced by inelasticity, bringing the granular gas back to the realm of NS hydrodynamics. In fact, theory-simulation comparisons seem to indicate that a part of the disagreement between theory and simulation can be ruled out when we consider an improve in the approximation to the reference distribution function in the Chapman-Enskog method (when a Sonine polynomial expansion to first order [Burnett(1935), Truesdell(1973)] for the homogeneous cooling state [van Noije and Ernst(1998)] is used instead of the Maxwellian [Garzó et al.(2007)Garzó, Santos, and Montanero]) and thus, the relatively small disagreement would not be due to an eventual failure of hydrodynamics but to the inherent weakness of the perturbative solution used in the Chapman-Enskog procedure [Chapman and Cowling(1970), Brilliantov and Pöschel(2004)] for calculation of the transport coefficients.

Therefore, there is no a priori reason for concern respect to the validity of hydrodynamics for granular gases. This would depend essentially on the specific properties and geometry of the flow of interest, not just the inelasticity.

Considering a granular gas fluidized by an active fluid should be regarded, not as an artifact, but as a natural situation since, in the first place, rapid granular flows are possible just because an energy input is able to compensate for inelastic cooling (otherwise, clustering instabilities and total freezing of the dynamics occur [Goldhirsch(2003)]). However, an interstitial fluid is a very frequent configuration for granular dynamics problems present in nature [Bagnold(1954)]. Applications to biology problems are therefore promising. Moreover, it has been experimentally proved that this type of fluctuating force like the one we use to model the active fluid reaches local equilibrium with a granular particle [Ojha et al.(2004)Ojha, Lemieux, Dixon, Liu, and Durian]. Based on these observations, our model for energy input from the active fluid should be closer, in most cases, to experimental situations [Mognetti et al.(2013)Mognetti, Šarić, Angioletti-Uberti, Cacciuto, Valeriani, and Frenkel] where, interestingly, particle clustering instabilities also may occur.

In addition, the results in the present work have helped us to produce further results, since we have already applied NS hydrodynamics with a non-uniform stochastic force also to problems and applications like granular segregation [Vega Reyes et al.(2014)Vega Reyes, Garzó, and Khalil], where excellent agreement is found between granular impurity segregation criteria resulting from NS theory and computer simulations. Finally, we plan to extend this work by developing the corresponding non-Newtonian hydrodynamic theory and also by calculating the linear stability criteria of theses base flows. On the other hand, this work could be relevant for the study of transport in crowded active enviroments [Bechinger et al.(2016)Bechinger, Leonardo, Lowen, Reichhardt, and Volpe] where the time scales of the active fluids is faster than the inert particles relaxation time. The existence of a well defined NS hydrodynamics in heterogeneous active bath [Argun et al.(2016)Argun, Moradi, PinÃ§e, Bagci, Imparato, and Volpe], is also an important result in order to explain biological systems. Another possible future work is to study the effect of sorting of active swimmers with inert particles in heterogeneous active media. Also the study of the behaviour of particles, mean square displacement, mobility, diffusion coefficient depending on the media. A study of the hydrodynamics of more complex forms of the volume force, for example with time dependency, such as for a traveling wave [Sándor et al.(2017a)Sándor, LibÃ l, Reichhardt, and Reichhardt], is also a promising line of research.

### .1 Computer simulations

As we said, we used for this work data from the direct simulation Monte Carlo method (DSMC) [Pöschel and Schwager(2005)] of the kinetic equation (6). The standard DSMC algorithm consists of two basic steps: collisions (particle-particle collisions and, in this case, particle-boundaries collisions) and free drift [Bird(1994)]. The boundaries are modeled as hard walls, similarly to previous works. In the presence of a volume force, there is an additional step that accounts for the action of the fluctuating volume force. Implementation of this step, is described for instance in [Montanero and Santos(2000), Gradenigo et al.(2011)Gradenigo, Sarracino, Villamaina, and Puglisi]. However, for this work, we consider a more general situation where the noise intensity is space-dependent. More specifically, the fluctuating force intensity has a profile obeying equation (22); i.e., of the form . For the solution to be self-consistent, the temperature profile should be obtained, as we explained, from the differential equation (21), that is the same differential equation that the ’LTu’ profiles obey [Vega Reyes et al.(2010)Vega Reyes, Santos, and Garzó]. For this reason, we have extracted simulation temperature profiles from LTu states obtained in a previous work [Vega Reyes et al.(2010)Vega Reyes, Santos, and Garzó]. We introduced them in condition (22), which defines a force intensity profile that is maintained constant in time during the simulation: . The initial temperature of the granular gas was however always set to the lower wall temperature (constant): . Once the simulation starts, and independently of the value of inelasticity, we always obtained after a transient a steady temperature profile with the same form of a Fourier flow, as expected (see Fig. 1). We have considered in all simulation figures presented in this work a system with and (, being the average density in the system). We also performed additional series with different (not shown in figures), and checked that with this variation we obtained the same results.

### .2 Navier-Stokes Transport coefficients

The Navier-Stokes (NS) transport coefficients that we use in this work have been calculated by Garzó and Montanero in a previous work [Garzó and Montanero(2002)]. Their expressions need to be inserted in Eqs. (19), (22) in order to complete the analytical solution of the hydrodynamic profiles, that result from the differential equations (15), (16) for sheared states and (20), (21) for reference states (no shear). We write these transport coefficients in this appendix too, for the sake of completion. The reduced viscosity , and reduced heat flux coefficients are [Garzó and Montanero(2002)]

(25) |

Here, , , are the collisional frequencies associated to the transport coefficients and their expressions were calculated by Brey and Cubero, together with the expression for the reduced cooling rate [Brey and Cubero(2001)]

(26) | |||

(27) | |||

(28) |

The coefficient in (25)-(28) was calculated by van Noije and Ernst [van Noije and Ernst(1998)] and has the expression

(29) |

As we already mentioned, the additional (dashed) theoretical curves in Fig. 2 correspond to the transport coefficients calculated with an improved approximation to the reference distribution function in the Chapman-Enskog method. In effect, in this improved procedure the reference distribution function an approximate form of the homogeneous state. This approximate form consists in an expansion in Sonine polynomials, where we only retain the first non-Maxwellian contribution. For more reference on this alternative procedure, please refer to [Garzó et al.(2007)Garzó, Santos, and Montanero] or [Garzó et al.(2009)Garzó, Vega Reyes, and Montanero].

###### Acknowledgements.

Support of the Spanish Government through Grant FIS2010-12587 (F. V. R.), MTM2014-56948-C2- 2-P (A. L.) is acknowledged. Partial support from FEDER funds and by Junta de Extremadura through Grant No. GRU10158 is also acknowledged (F. V. R.). A. L. thanks the hospitality of UEX where during a stay this work was done.### References

- Marchetti, M.C.; et al.. Hydrodynamics of soft active matter. Rev. Mod. Phys. 2013, 85, 1143.
- Bertin, E.; Droz, M.; Grégoire, G. Hydrodynamic equations for self-propelled particles: microscopic derivation and stability analysis. J. Phys. A 2009, 42, 445001.
- Bertin, E.; Droz, M.; Grégoire, G. Boltzmann and hydrodynamic description for self-propelled particles. Phys. Rev. E 2006, 74, 022101.
- Chou, Y.; Wolfe, R.; Ihle, T. Kinetic theory for systems of self-propelled particles with metric-free interactions. Phys. Rev. E 2012, 86, 021120.
- Ihle, T. ChapmanâEnskog expansion for the Vicsek model of self-propelled particles. J. Stat. Mech. 2016, p. 083205.
- Vicsek, T.; CzirÃ³k, A.; Ben-Jacob, E.; Cohen, I.; O.Shochet. Novel Type of Phase Transition in a System of Self-Driven Particles. Phys. Rev. Lett. 1995, 75, 1226.
- Bechinger, C.; Leonardo, R.D.; Lowen, H.; Reichhardt, C.; Volpe, G. Active particles in complex and crowded enviroments. Rev. Mod. Phys. 2016, 88, 045006.
- Argun, A.; Moradi, A.R.; PinÃ§e, E.; Bagci, G.B.; Imparato, A.; Volpe, G. Non-Boltzmann stationary distribution and nonequilibrium relations in active baths. Phys. Rev. E 2016, 94, 062150.
- Sándor, C.; LibÃ l, A.; Reichhardt, C.; Reichhardt, C.J.O. Collective transport for active matter run-and-tumble disk systems on a travellin-wave substrate. Phys. Rev. E 2017, 95, 012607.
- Sándor, C.; LibÃ l, A.; Reichhardt, C.; Reichhardt, C.J.O. Dynamic phases of active matter systems with quenched disorder. Phys. Rev. E 2017, 95, 032606.
- Sándor, C.; LibÃ l, A.; Reichhardt, C.; Reichhardt, C.J.O. Dewetting and spreading transitions for active matter on random pinning substrate. J. Chem. Phys. 2017, 146, 204903.
- Angelini, T.E.; Hanezzo, E.; Trepat, X.; Marquez, M.; Fredberg, J.J.; Weitz, D.A. Glass-like dynamics of collective cell migration. Proc. Nat. Acad. Sci. 2011, 108, 4714.
- Malinverno, C.; Corallino, S.; et al., F.G. Endocytic reawakening of motility in jammed epithelia. Nature Materials 2017, 16, 587.
- Garcia, S.; Hannezo, E.; Elgeti, J.; Joanny, J.F.; Silberzan, P.; Gov, N.S. Physics of active jamming during collective cellular motion in a monolayer. PNAS 2015, 14, 1040.
- Bi, D.; Lopez, J.H.; Schwarz, J.M.; Manning, M.L. A density-independent glass transition in biological tissues. Nature Physics 2015, 11, 1074.
- Malgaretti, P.; Pagonabarraga, I.; Rubi, M.J. Rectification and Non-Gaussian Diffusion in Heterogeneous Media. Entropy 2016, 18, 349.
- Dufty, J.W. Kinetic theory and hydrodynamics for a low density granular gas. Adv. Complex Sys. 2001, 4, 397.
- Brilliantov, N.V.; Pöschel, T. Kinetic Theory of Granular Gases; Oxford University Press, Oxford, 2004.
- Goldhirsch, I. Rapid granular flows. Annu. Rev. Fluid Mech. 2003, 35, 267–293.
- de Bruyn, J. Physics 2011, 4, 86.
- Batchelor, G.K. The theory of homogeneous turbulence; Cambridge University Press, Cambridge, UK, 1953.
- Haff, P.K. Grain flow as a fluid-mechanical phenomenon. J. Fluid Mech. 1983, 134, 401–430.
- Puglisi, A.; Gnoli, A.; Gradenigo, G.; Sarracino, A.; Villamaina, D. Structure factors in granular experiments with homogeneous fluidization. J. Chem. Phys. 2012, 136.
- Puglisi, A. Transport and Fluctuations in Granular Fluids; Springer International Publishing, 2014.
- Jaeger, H.M.; Nagel, S.R. Science 1992, 255, 1523–1531.
- Geminard, J.C.; Losert, W.; Gollub, J.P. Frictional mechanics of wet granular material. Phys. Rev. E 1999, 59, 5881.
- Williams, D.R.M.; MacKintosh, F.C. Driven granular media in one dimension: Correlations and equation of state. Phys. Rev. E 1996, 54, R9.
- Cafiero, R.; Luding, S.; Herrmann, H.J. Two-dimensional granular gas of inelastic spheres with multiplicative driving. Phys. Rev. Lett. 2000, 84, 6014.
- Visco, P.; Puglisi, A.; Barrat, A.; Trizac, E.; van Wijland, F. Fluctuations of Power Injection in Randomly Driven Granular Gases. J. Stat. Phys. 2006, 125, 533.
- Villamaina, D.; Puglisi, A.; Vulpiani, A. The fluctuation–dissipation relation in sub-diffusive systems: the case of granular single-file diffusion. J. Stat. Mech. 2008, p. L10001.
- Fiege, A.; Aspelmeier, T.; Zippelius, A. Long-Time Tails and Cage Effect in Driven Granular Fluids. Phys. Rev. Lett. 2009, 102, 098001.
- Gradenigo, G.; Sarracino, A.; Villamaina, D.; Puglisi, A. Fluctuating hydrodynamics and correlation lengths in a driven granular fluid. J. Stat. Mech. 2011, p. P08017.
- Losert, W.; Bocquet, L.; Lubensky, T.C.; Gollub, J.P. Particle Dynamics in Sheared Granular Matter. Phys. Rev. Lett. 2000, 85, 1428.
- Yu, H.; et al.. The fertilizing role of African dust in the Amazon rainforest: A first multiyear assessment based on data from Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations. Geophy. Res. Lett. 2015, 42, 1984–1991.
- Ojha, R.P.; Lemieux, P.A.; Dixon, P.K.; Liu, A.J.; Durian, D.J. Statistical mechanics of a gas-fluidized particle. Nature 2004, 427, 521.
- D’Anna, G.; Mayor, P.; Barrat, A.; Loreto, V.; Nori, F. Observing brownian motion in vibration-fluidized granular matter. Nature 2003, 424, 909–912.
- Garzó, V.; Tenetti, S.; Subramaniam, S.; Hrenya, C.M. Enskog kinetic theory for monodisperse gasâsolid flows. J. Fluid Mech. 2012, 712, 129–168.
- Lun, C.K.K.; Savage, S.B., Granular Gas Dynamics; Springer, 2003; chapter Kinetic Theory for Inertia Flows of Dilute Turbulent Gas-Solids Mixtures, p. 263. ed. T. Pöschel and N. Brilliantov.
- Campbell, C.S. Rapid Granular Flows. Annu. Rev. Fluid Mech. 1990, 22, 57–92.
- Vega Reyes, F.; Santos, A.; Garzó, V. Non-Newtonian granular hydrodynamics: what do the inelastic simple shear flow and the elastic Fourier flow have in common? Phys. Rev. Lett. 2010, 104, 028001.
- Foerster, S.F.; Louge, M.Y.; Chang, H.; Alla, K. Measurements of the collision properties of Small Spheres. Phys. Fluids 1994, 6, 1108–1115.
- Grasselli, Y.; Bossis, G.; Morini, R. Measurements of the collision properties of Small Spheres. Eur. Phys. J. E 2015, 38, 8.
- Van Kampen, N.G. Stochastic Processes in Physics and Chemistry; Elsevier, Amsterdam, 1992.
- Marini-Bettolo-Marconi, U.; Tarazona, P.; Cecconi, F. Theory of thermostatted inhomogeneous granular fluids: A self-consistent density functional description. J. Chem. Phys. 2007, 126, 164904.
- Sarracino, A.; Villamaina, D.; Constantini, G.; Puglisi, A. Granular Brownian motion. J. Stat. Mech. 2010, p. P04013.
- Garzó, V.; Chamorro, M.G.; Vega Reyes, F. Transport properties for driven granular fluids in situations close to homogeneous steady states. Phys. Rev. E 2013, 87, 032201.
- Marconi, U.M.B.; Puglisi, A.; Rondonic, L.; Vulpiani, A. Fluctuation–dissipation: Response theory in statistical physics. Phys. Rep. 2011, 461, 111.
- Brey, J.J.; Dufty, J.W.; Kim, C.S.; Santos, A. Hydrodynamics for granular flow at low density. Phys. Rev. E 1998, 58, 4638–4653.
- Bird, G.I. Molecular Gas Dynamics and the Direct Simulation of Gas Flows; Oxford University Press, Oxford, 1994.
- Vega Reyes, F.; Urbach, J.S. Steady base states for Navier–Stokes granular hydrodynamics with boundary heating and shear. J. Fluid Mech. 2009, 636, 279.
- Burnett, D. Velocity Distribution in a non-uniform gas. Proc. Lond. Math. 1935, 39, 385.
- Chapman, C.; Cowling, T.G. The Mathematical Theory of Non-Uniform Gases, 3 ed.; Cambridge University Press, Cambridge, 1970.
- Abramowitz, M.; Stegun, I.A. Handbook of Mathematical Functions; Dover Publications, Inc., New York, 1965.
- Vega Reyes, F.; Santos, A.; Garzó, V. Steady base states for non-Newtonian granular hydrodynamics. J. Fluid Mech. 2013, 719, 431–464.
- Garzó, V.; Montanero, J.M. Transport coefficients of a heated granular gas. Physica A 2002, 313, 336.
- Boyd, I.D.; Chen, G.; Candler, G.V. Phys. Fluids 1995, 7, 210.
- Campbell, C.S. The stress tensor for simple shear flows of a granular material. J. Fluid Mech. 1989, 449, 58.
- Garzó, V.; Santos, A.; Montanero, J.M. Modified Sonine approximation for the Navier– Stokes transport coefficients of a granular gas. Physica A 2007, 376, 94.
- Truesdell, C. Mathematical aspects of the kinetic theory of gases; Instituto de Matematica, Universidade Federal do Rio de Janeiro, 1973.
- Vega Reyes, F.; Santos, A.; Kremer, G.M. Role of roughness on the hydrodynamic homogeneous base stateof inelastic spheres. Phys. Rev. E 2014, 89, 020202(R).
- van Noije, T.; Ernst, M. Velocity distributions in homogeneous granular fluids: Velocity distributions in homogeneous granular fluids: the free and the heated case. Gran. Matt. 1998, 1, 57.
- Bagnold, R.A. The physics of Blown Sand and Desert Dunes; Dover Publications Inc., Mineola, New York, 1954.
- Mognetti, B.M.; Šarić, A.; Angioletti-Uberti, S.; Cacciuto, A.; Valeriani, C.; Frenkel, D. Living Clusters and Crystals from Low-Density Suspensions of Active Colloids. Phys. Rev. Lett. 2013, 111, 245702.
- Vega Reyes, F.; Garzó, V.; Khalil, N. Hydrodynamic granular segregation induced by boundary heating and shear. Phys. Rev. E 2014, 89, 052206.
- Pöschel, T.; Schwager, T. Computational granular dynamics: models and algorithms; Springer-Verlag, Berlin, 2005.
- Montanero, J.M.; Santos, A. Computer simulation of uniformly heated granular fluids. Gran. Matt. 2000, 2, 53–64.
- Brey, J.J.; Cubero, D. In Granular Gases; Pöschel, T.; Luding, S., Eds.; Lecture Notes in Physics, Springer, Berlin, 2001; p. 59.
- Garzó, V.; Vega Reyes, F.; Montanero, J.M. Modified Sonine approximation for granular binary mixtures. J. Fluid Mech. 2009, 623, 387–411.