text:
# Clustering in a One-dimensional Inelastic Lattice Gas

## I Introduction

Dynamics of granular fluids has captured a lot of attention from theoretical physicists for the last few years. In a theorist’s model, the constituent particles of a granular fluid, usually considered to be hard spheres of finite radii, irretrievably dissipate kinetic energy via inelastic binary collisions and interparticle frictional forces. As a result, unlike the microscopic models for the classical kinetic theory of gases, a granular fluid that is not driven by external forces “cools freely”.

Even in the absence of frictional forces, in a stark contrast to hard sphere fluids with elastic inter-particle binary collisions, such a simplified model of granular fluids exhibits complex behaviour at macroscopic scales. In its simplest form, a freely cooling initially homogeneous and isotropic (both in the particles’ position and velocity space) inelastic gas spontaneously forms nontrivial structures in the macroscopic velocity as well as in the macroscopic density field of the gas. A large number of studies, mostly from the point of view of inelastic hydrodynamics, have been carried out to understand the onset for the formation of these structures in two and three dimensions (Refs. young ; young1 ; gold ; ernst1 ; brey1 ; soto ; twan  to cite a few). The (qualitative and quantitative) picture that has emerged from these studies is that in two and three dimensions, the system of inelastic hard spheres suffers from inherent long wavelength linear instabilities. When the system size allows such long wavelengths to be present, these instabilities start to generate inhomogeneities in the macroscopic velocity and the density field of the gas. In the subsequent evolution, these inhomogeneities interact nonlinearly to give rise to macroscopic structures and the entire system evolves into a collection of densely populated clusters that are separated by regions containing particles at very low density gold ; brey2 . At late times, the clusters collide in a very complex manner and merge — a phenomenon that is known as coarsening in the literature poschel ; luding .

The late time evolution of a freely cooling inelastic gas is thus qualitatively completely different from the linear instability mechanisms at early times. However, although a very large number of studies have been devoted to kinetic theory of freely cooling inelastic gases in two and three dimensions, due to the difficulties associated with the nonlinearities in the behaviour of individual clusters and cluster-cluster collisions, a proper theoretical understanding of the long time dynamics of freely cooling inelastic gases has remained elusive. The existing results have only been numerical young1 ; chen ; her ; hill ; poschel ; luding .

At the other extreme, fully analytical solutions have been found for completely inelastic (or “sticky”) granular gas (of point particles) in one dimension Frach1 , and it has been shown that the sticky gas in one dimension is described by the Burgers equation in the inviscid limit Frach1 . In addition, a recent experiment has also observed clustering in a one-dimensional granular gas ko . Due to dimensional reasons, the dynamics of a granular gas in one dimension is qualitatively different from those in two or three dimensions (e.g., vortices cannot form, a strict ordering of particles from left to right is maintained at all times), but structures are still seen to form in the velocity and as well as in the density field of the gas Kad1 ; BN . Ben-Naim et al BN  have studied the formation of these structures, and have conjectured that at the long times, the behaviour of a one-dimensional inelastic gas should be the same as that of the sticky gas.

Our purpose in this paper is to unravel some of these long-time phenomena at a more microscopic level in a simple one-dimensional lattice model that has been introduced in Ref. Bald ; Bald2 . In this model, one considers a system of $`N`$ particles on the integral lattice positions (denoted by $`k`$) of a ring of size $`N`$ with initial velocities chosen randomly from a uniform distribution in $`[1,1]`$. The ordering of the particles is maintained at all times; however, since the particles do not move in this lattice model, the velocities of the particles are not the time derivative of their positions. Instead, $`v_k`$, the velocity of the $`k`$th particle, is simply a scalar quantity associated with the $`k`$th particle, and $`v_k`$ changes only when the $`k`$th particle participates in a collision with one of its neighbours according to the following collision rule. At a collision between $`k`$th and $`(k+1)`$th particles at any time, the post-collisional velocities $`v_{k,(k+1)}^{(+)}`$ are related to their pre-collisional velocities $`v_{k,(k+1)}^{()}`$ by

$`v_{k,(k+1)}^{(+)}=v_{k,(k+1)}^{()}\pm {\displaystyle \frac{1+r}{2}}\left[v_{(k+1)}^{()}v_k^{()}\right],`$ (1)
where $`0r1`$ is the coefficient of restitution. Time is measured by the average number of collisions per particle in this lattice model.

In this paper, we will consider two variants of this lattice model, namely the “random lattice model” and the “systematic lattice model”. Both obey the collision rule (1), but they differ in the way a colliding pair of particles is chosen to update the system. For the random lattice model, at any update, the colliding pair of particles are chosen randomly from all particles that momentarily satisfy the kinetic constraint $`v_{k+1}v_k<0`$. On the other hand, for the systematic lattice model, the colliding pair of particles are the ones that has the momentary global minimum value of $`v_{k+1}v_k`$.

The random lattice model without kinetic constraint, has already been studied as part of a larger class of models, the inelastic Maxwell models Max , in which the collision frequency is chosen to be independent of the incoming velocities of colliding particles. For that model, it has been analytically shown that correlations develop with a diffusively growing correlation length, which consequently affects the temperature decay rate. In addition, the inelastic lattice model in one-dimension with the kinetic constraint, which we study here, was also analyzed in Refs. Bald ; Bald2  in terms of velocity distribution and structure factors remark . It is however important to realize the difference between these existing results and the ones reported in this paper: the existing results mainly concern global quantities, while in this paper, our main thrust is to study the behaviour of the microscopic inhomogeneities arising from the kinetic constraint.

The key feature of both variants of the lattice model with kinetic constraint studied here is that an initial configuration of random velocities of the particles soon develops distinct spatial structures, eventually leading to large positive sharp jumps in the particle velocities intersparsed with relatively smooth variations. We refer to the large positive sharp jumps in the velocity field as shocks and the region between two consecutive shocks, where the velocity variations are smooth, as clusters. These structures have already been observed in Bald , but have not been fully analyzed. Formation of shocks and the subsequent dynamics of the clusters make the lattice model interesting on its own, but its relevance is realized only when the locations of its shocks in $`v(k)`$-profile are compared to those of the one-dimensional inelastic gas Frach1 ; BN  (see Fig. 1). By contrast, there are no shocks in the one dimensional random lattice model without kinetic constraints Max .

This paper is organized in the following manner: in Sec. II, we discuss the generic features of the lattice model and explore its connections with the one-dimensional inelastic gas Frach1 ; BN . In Sec. III, we analyze the formation of shocks at early times and cluster dynamics at late times. We finally end the paper with a short discussion in Sec. IV.

## II Generic features of the lattice model

### II.1 Phenomenology of shock development and the subsequent dynamics for the random lattice model

To start with, in Fig. 2, we show a time sequence of the velocity profile $`v(k)`$ for the random lattice model of $`10000`$ particles with $`r=0.7`$. The initial configuration \[Fig. 2(a)\] is created by choosing $`v_k`$ randomly from a uniform distribution in $`[1,1]`$. As can be seen in Fig. 2(b), velocity correlations set in very rapidly (within $`10`$ collisions per particle). After $`10^4`$ collisions per particle, shocks and clusters can be clearly identified \[Fig. 2(c)\]. According to the kinetic constraint, neighbouring particles cannot collide across a shock, and as a result, each cluster evolves independently of the others, until two neighbouring clusters collide and coalesce to form a new bigger cluster \[this mechanism is illustrated by the evolution of velocity profile of Fig. 2(c) to that of 2(d); Fig. 2(d) corresponds to a state after $`10^6`$ collisions per particle\].

From the above phenomenological description, one can identify two separate regimes: the initial homogeneous regime, and the subsequent clustered regime at late times. To make the evolution of the system from one regime towards the other quantitatively more precise, we notice that it corresponds to a symmetry breaking of the probability distribution $`P(\mathrm{\Delta }v_k,t)`$ of relative velocities $`\mathrm{\Delta }v_k=v_{k+1}v_k`$. Initially, for each $`k`$, the distribution of $`v_k`$ is uniform in $`[1,1]`$ (i.e., a “box” distribution). Naturally, $`P(\mathrm{\Delta }v_k,t=0)`$ is a piecewise linear distribution on $`[2,2]`$ (a “triangular” distribution), symmetric around $`\mathrm{\Delta }v_k=0`$. Under the effect of the dynamics, this symmetry is not preserved. This is easily understood from the fact that only the neighbouring particle pairs with $`\mathrm{\Delta }v_k<0`$ can collide, and after a collision between the $`k`$th and $`(k+1)`$th particle, $`\mathrm{\Delta }v_{k1}`$, $`\mathrm{\Delta }v_k`$ and $`\mathrm{\Delta }v_{k+1}`$ change in a way that, in general, does not preserve the symmetry of $`P(\mathrm{\Delta }v_k,t)`$. As this process continues, from phenomenological considerations, we know that a few shocks remain at late times. Now, since shocks correspond to large values of $`\mathrm{\Delta }v_k`$, this means that in the clustered regime, the probability to have a large positive value of $`\mathrm{\Delta }v_k`$ is expected to be greater than that of a correspondingly large negative value.

To detect the locations of the shocks by following the evolution of $`P(\mathrm{\Delta }v_k,t)`$ as a function of time, one can identify the location $`k`$ of a shock by requiring $`\mathrm{\Delta }v_k>C|\text{min}_j\mathrm{\Delta }v_j|`$, where $`C>0`$ is a constant. The choice of the numerical value of $`C`$ is arbitrary, and this arbitrariness can be used to tune the minimum value that $`\mathrm{\Delta }v_k`$ must have in order to qualify for a shock. Furthermore, if one defines the instant $`t_p`$ when the system passes from one regime to the other by the minimal time where at least one shock becomes visible, then the above requirement can also be used to characterize $`t_p`$. Of course, the precise value of $`t_p`$ depends on the chosen value of $`C`$.

In the clustered regime, due to inelastic collisions, the amplitude of the velocity profile within each cluster decreases, and the velocities of the particles approach that of the center of mass of the cluster itself. This process continues for a while until two neighbouring clusters, with center of mass velocities $`V_j`$ and $`V_{j+1}`$ collide and coalesce and the shock separating them disappears \[at time $`t`$, if we number the clusters $`j=1,2,\mathrm{},j_1(t)`$ from the left to the right, then such a cluster-cluster collision takes place only if $`V_j>V_{j+1}`$\]. Precisely this mechanism is responsible for making the velocity profile of Fig. 2(c) evolve to that of Fig. 2(d).

At very late times, due to periodic boundary condition, only a single cluster survives. Its amplitude also decreases in time, and the velocities of all particles asymptotically converge to the center of mass velocity of the whole system.

### II.2 Dependence on the coefficient of restitution for the random lattice model

One of the first questions that one can ask for the lattice model is what is the effect of inelasticity on the clustering properties of the system. Clearly, for $`r=1`$, the collisions are elastic, and upon collision, the velocities of two neighbouring particles are simply interchanged. Such a dynamics cannot lead to clustering or any structure formation.

On the other hand, our simulations show that the lattice model exhibits clustering for any value of $`r`$ smaller then $`1`$, and it is only $`t_p`$ defined in Sec. II.1 that becomes an increasing function of $`r`$. The most remarkable observation, however, is that for a given initial velocity configuration of the particles, not only the locations of the shocks persisting at late times are the same for all $`r`$, but also the macroscopic velocity profiles are nearly independent of $`r`$. If we consider the macroscopic velocity profile of the lattice model with a coefficient of restitution $`r_1`$ at a time $`t_1`$ well into the clustered regime, then for any $`r_2>r_1`$, (originating from the same initial velocity configuration of the particles) there exists a time $`t_2<t_1`$, for which the two macroscopic velocity profiles are almost the same (see Fig. 3). Effectively, this implies that for all $`r`$, the dynamical behaviour of the system is identical to that for $`r=0`$ at late times. We also observe that at intermediate times, the maximum number of shocks (and correspondingly, clusters) in the system are observed for $`r=0`$ (see Sec. III.1 in this regard).

### II.3 Dependence on the order of collisions and on the initial velocity configuration of the particles

The next important question one can ask is if not $`r`$, exactly what influences the positions of the shocks? A priori, one expects that the locations of the shocks and the long time dynamics of the system depend on the initial velocity configuration of the particles, as well as on the sequence in which the collisions are performed.

Surprisingly however, the locations of the shocks depend only on the initial configuration of the particles’ velocities, e.g, two different random collision sequences result in the same locations of the shocks. On the other hand, the variant of the lattice model (random or systematic) determines the macroscopic shape of the velocity profile inside the clusters. To demonstrate this phenomenon, starting from the same initial velocities of the particles, we simulated both the random and the systematic lattice models and compared the locations of the shocks (see Fig. 1).

### II.4 Relation between the lattice model and one-dimensional inelastic gas

The empirical observations in the preceding sections suggest that so long as one is only interested in the velocity profile $`v(k)`$, a close relation exists between the one-dimensional inelastic gas and the lattice model from the following consideration: the dynamics of a one-dimensional inelastic gas for a given initial configuration is completely deterministic — after any collision, the next colliding pair is automatically determined by the instantaneous minimal value of $`(x_{k+1}x_k)/(v_kv_{k+1})`$, for $`v_kv_{k+1}>0`$. In this sense, one can view the one-dimensional inelastic gas as a variant of the lattice model too, where the colliding particles are chosen in a very complicated manner. Section II.3 then indicates that for a given initial configuration of the particles’ velocities as a function of the particle index, the locations of the shocks in the one-dimensional inelastic gas and those in the lattice model are the same. That this is indeed the case is numerically verified in Fig. 1.

In our simulations of the one-dimensional inelastic gas, we followed the procedure outlined in Ref. BN , namely that to avoid the inelastic collapse, the collisions were taken to be elastic when the relative velocities between the neighbouring particles are smaller than a given cutoff. Our results confirm the existence of “Burgers” shocks in the $`v(x)`$ profile BN . More importantly, we note that for the one-dimensional inelastic gas, there is an inverse correspondence of shocks and clusters between this $`v(x)`$ and $`v(k)`$ profile — each shock in the $`v(x)`$ profile \[i.e., a large particle density $`n(x)`$\] corresponds to a cluster in the $`v(k)`$ profile and the shocks in the $`v(k)`$ profile correspond to regions with finite gradients in $`v(x)`$.

## III Analysis

### III.1 Development of shocks

The lesson that we learnt from Sec. II is that the locations of the shocks in this lattice model are essentially determined from the initial configuration of the particles’ velocities. This immediately gives rise to the following question: how can one predict the locations of the shocks seen at early times from the initial velocity configuration of the particles?

As it turns out, most of the shocks observed in the clustered regime already exist at $`t=0`$. These are the locations marked by large positive jumps in the $`v(k)`$ profile between two neighbouring lattice sites at $`t=0`$, surrounded by relatively small velocity variations. To illustrate this point, we present a schematic diagram for such a velocity profile for a few particle in Fig. 4. In the subsequent dynamics at short times, the particles on the right of B and on the left of A collide with each other and “thermalize” (i.e., the velocity variations reduce due to inelasticity), but this thermalization may not necessarily convert a relatively large $`v_Bv_A>0`$ to $`v_Bv_A<0`$ quickly enough to make the A and B collide. In other words, such a large positive jump remains preserved and may eventually give rise to a shock.

To check if such a scenario is correct or not, one needs to identify the locations of such large positive jumps in the initial $`v(k)`$ profile, and contrast them with the observed locations of the shocks at early times. A convenient way to single out such jumps \[say, in Fig. 2(a)\] is to define a $`(2n+1)`$-coarse grained velocity profile $`\overline{v}(k)={\displaystyle \underset{i=kn}{\overset{k+n}{}}}v_i/(2n+1)`$ of the particles, which, roughly speaking, can single out one such large jump over the surrounding small velocity variations in a window of $`(2n+1)`$ lattice sites in the $`\overline{v}(k)`$ profile. Having fine-tuned $`n`$, a surprisingly close match between the observed shock locations at early times and the initial configuration of the particles’ velocities can be found (see Fig. 5).

Notice, however, in Fig. 5 that not all the large jumps in the initial configuration of the particles’ velocities have turned out to be shocks at a later time. Indeed, in general, whether such an initial large positive jump at a particular location develops into a shock or not really depends on the magnitude of the jump itself in relation to particle velocities in its immediate vicinity, and as well as on the coefficient of restitution $`r`$. In this regard, the dependence on the coefficient of restitution is not difficult to understand, as $`r`$ quantifies the “transport of velocity” between two colliding particles. Intuitively speaking, an inelastic collision between two neighbouring particles can be viewed as a combination of dissipation (thermalization) and transport of velocities. For $`r=1`$ only transport of velocities can take place. At the other limit, for $`r=0`$, there is no transport of velocities but only dissipation. As a result, for $`r=0`$, small velocity variations around any large jump thermalize immediately, leading to an early appearance of the clustered regime, and a large number of initial jumps end up becoming shocks. With increasing $`r`$, some of the initial large positive jumps are eliminated by transport of velocities, and clustering regime, with a smaller number of selected shocks (and clusters), appears later.

The appearance of the shocks at the locations where relatively large positive jumps exist is therefore caused by the thermalization of the smaller velocity variations around these jumps — due to thermalization, the smaller velocity variations become even smaller, but the large positive jumps remain preserved under the dynamics. In relation to their surroundings, the magnitude of these large jumps thus starts to grow in time. Such a scenario applies independently of the details of the model, and gives a qualitative explanation for the observation that starting from a given initial configuration, shocks appear at the same locations in the $`v(k)`$ profile for both variants of the lattice model and as well as for the one-dimensional inelastic gas Frach1 ; BN . Moreover, it also clearly demonstrates that the shock formation “instability” in the $`v(k)`$ profile in all these models is a “relative instability” and not an “absolute” one. Such behaviour has also been found in the (long-wavelength) linear instability of the macroscopic flow field for two- and three-dimensional inelastic gases (see Refs. gold ; ernst1 ; twan  for example).

### III.2 Cluster dynamics

#### III.2.1 Dynamics within a single cluster

Once the shocks develop, our simulations show that for a given variant of the lattice model, the velocity profile within each cluster has the same characteristic macroscopic shape. For the deterministic lattice model, this shape is simply linear with a negative slope. Although we cannot provide an analytical derivation of this shape, if we interpret the dynamics along the lines that at each updating step of a cluster, two neighbouring particles collide at the location where the $`v_{k+1}v_k`$ is minimum (which is a mechanism that tries to align all $`v_{k+1}v_k`$ values to its maximum value within each cluster), then the linear profile appears to be intuitively reasonable.

For the random lattice model, however, the slope of the macroscopic velocity profile is smoothly varying, as can be seen in, e.g., Fig. 2(d) \[occasionally, flat profiles can be observed too, but for the time being, we leave them for Sec. III.2.2\]. It turns out that the functional form of the smooth macroscopic velocity profile can be obtained through a mean-field approach by averaging over an ensemble of random collision sequence realizations within a cluster, as we describe below.

The idea behind this mean-field approach is the following: we consider a cluster of $`M`$ particles, and denote the velocity of the $`k`$th particle at time $`t`$, averaged over random collision sequences, by $`v_k(t)`$. From the very definition of a cluster in the introduction, we assume that $`v_{k+1}(t)v_k(t)<0`$ at all times. Thereafter, as we observe from direct simulation measurements that the probability of collisions between any two neighbouring particle pairs are equally likely, in this mean field approach, $`v_k(t)`$ is easily seen to satisfy the equation

$`v_k(t+1)=(12p)v_k(t)`$ (2)
$`+p[v_k(t)\epsilon \{v_k(t)v_{k1}(t)\}]`$
$`+p[v_k(t)+\epsilon \{v_{k+1}(t)v_k(t)\}].`$
Here, $`p`$ is an effective probability of a collision between any two neighbouring particles in this mean-field theory and $`\epsilon =(1+r)/2`$. The three terms on the r.h.s. of Eq. (5) respectively originate from the events when the $`k`$th particle is not involved in a collision, when there is a collision between the $`k`$th and the $`(k1)`$th particle and when there is a collision between the $`k`$th and the $`(k+1)`$th particle. As Eq. (2) very simply reduces to

$`v_k(t+1)=(12p\epsilon )v_k(t)+p\epsilon v_{k1}(t)`$ (3)
$`+p\epsilon v_{k+1}(t),`$
the interesting point to note is that Eq. (3) is the discrete (both in space and time) form of the diffusion equation

$`{\displaystyle \frac{v(k,t)}{t}}=D(r){\displaystyle \frac{^2v(k,t)}{k^2}},`$ (4)
where mean-field theory predicts $`D(r)=p\epsilon `$. One has to keep in mind however that Eq. (4) holds only for monotonically decreasing $`v(k,t)`$ as a function of $`k`$. A similar equation has also been found in the one-dimensional lattice model without the kinetic constraint Max .

To solve Eq. (4), the macroscopic velocity profile of a cluster at time $`t_0`$ can be generally expanded in a Fourier series. Since particle velocities are not transported across the boundaries of a cluster, $`{\displaystyle \frac{v_k(t)}{k}}`$ must vanish at the boundaries of a cluster, the Fourier series contains only terms $`\mathrm{cos}[\pi jk]`$. The amplitude of the $`j`$th such term decreases as $`e^{\pi ^2j^2D(r)(tt_0)}`$, so that for large values of $`(tt_0)`$, only the slowest decay mode $`j=1`$ survives. Hence at long times we expect $`v(k,t)`$ to be given by

$`v(k,t)=A(t)\mathrm{cos}\left[{\displaystyle \frac{\pi k}{M}}\right]`$ (5)
$`=A(t_0)e^{\pi ^2D(r)(tt_0)/M^2}\mathrm{cos}\left[{\displaystyle \frac{\pi k}{M}}\right].`$
The mean-field result (5) allows us to compare the half-cosine form of $`v(k,t)`$ with the macroscopic velocity profile of the particles within a cluster, as observed in the simulation for the random lattice model. This comparison is shown in Fig. 6(a) — the fact that we cannot distinguish the simulation data from the half-cosine shape of $`v(k,t)`$ as predicted by Eq. (5) is an indication of how well the mean-field approach works to describe the average shape of an isolated cluster. In addition, by following the dynamics from the simulation, we can also obtain an empirical functional form of $`D(r)`$ \[shown in Fig. 6(b)\].

There are two features of the $`D(r)`$ vs. $`r`$ curve that require further elaboration. The first one of them is that $`D(r)`$ is an increasing function of $`r`$, indicating that the amplitude of the cluster $`A(t)`$ decreases faster as $`r`$ increases. This observation seems to contradict the fact that the dissipation decreases with increasing $`r`$. The point to notice however is that the mean-field approach really describes the time evolution within a cluster in which $`{\displaystyle \frac{v_k(t)}{k}}<0`$. From that point of view, the amplitude of the cluster decreases not only from the dissipation, but also from transport of velocities. This is easily seen from the fact that after a collision, the ordering of the magnitudes of the velocities of the colliding particles are simply exchanged, i.e., the particle on the right has a higher velocity, which results in a local positive slope in the $`v(k)`$ profile. For increasing $`r`$, transport of velocities become more effective, and thus the negative slope of $`{\displaystyle \frac{v_k(t)}{k}}`$ within a cluster decays to zero faster. Only if $`r=1`$, the velocities of the particles within an isolated cluster order themselves in an increasing order of magnitude. However, $`D(r)`$ has a singularity at $`r=1`$ and the $`r=1`$ case cannot be treated within the scope of our mean-field theory. Secondly, from the way we presented our mean-field approach, it may seem that the $`r`$-dependence of $`D(r)`$ can be obtained through the dependence of $`p`$ and $`\epsilon `$ on $`r`$. This is actually not true. The mean-field approach (2) neglects fluctuations in the particles’ velocities around $`v(k,t)`$, but in an actual simulation, these fluctuations are very important as they decide the sign of $`\mathrm{\Delta }v_k(t)=v_{k+1}(t)v_k(t)`$ for the particles within the cluster and thereby control which collisions are possible and which are not.

#### III.2.2 Interacting clusters for the random lattice model

With the background of Sec. III.2.1, it is now clear that as time progresses, the amplitude of each cluster present in the random lattice model at late times effectively decays exponentially as $`\mathrm{exp}[\pi ^2D(r)(tt_0)/M^2]`$. Such a decay brings the velocities of the particles within a cluster closer and closer to the centre-of-mass velocity of the cluster itself. Let us assume that at time $`t_01`$, there are $`j_1`$ number of clusters present in the whole system of the random lattice model. If we number these clusters by $`j`$ such that $`j=1,2,\mathrm{},j_1`$ and denote the number of particles within the clusters, the amplitudes of the clusters and their velocities of their centre of masses respectively by $`M_j`$, $`A_j`$ and $`V_j`$, then with increasing time, we observe the following dynamics: (i) when $`V_{j1}<V_j<V_{j+1}`$, the shocks on the left and on the right of the $`j`$th cluster cannot respectively be smaller than $`V_jV_{j1}>0`$ and $`V_{j+1}V_j>0`$. In that case, the amplitude $`A_j`$ simply decreases to $`0`$, forming a “flat” velocity profile (see e.g., Fig. 3). (ii) On the other hand, when $`V_j>V_{j+1}`$, the magnitude of the shock between the $`j`$th and the $`(j+1)`$th clusters decreases to zero, and the two clusters coalesce together to form a new bigger cluster in a finite time. For a given configuration of clusters at time $`t`$, the mean-field theory of Sec. III.2.1 provides us with a way to measure the time $`t^{}=t+\mathrm{\Delta }t_{min}`$ of the first coalescence of two clusters in the system. In fact, $`\mathrm{\Delta }t_{min}=\text{min}_j\mathrm{\Delta }t_j`$, where for each pair of neighbouring clusters $`(j,j+1)`$ with $`V_jV_{j+1}>0`$, $`\mathrm{\Delta }t_j`$ is obtained by solving the equation

$`A_{j+1}(t_0)\mathrm{exp}[\pi ^2D(r)\mathrm{\Delta }t_j/M_{j+1}^2]`$ (6)
$`+A_j(t_0)\mathrm{exp}[\pi ^2D(r)\mathrm{\Delta }t_j/M_j^2]=V_jV_{j+1}.`$
Due to the fact that Eq. (6) is transcendental in nature, a closed form analytical solution for $`\mathrm{\Delta }t_j`$ is impossible to obtain, let alone the value of $`\mathrm{\Delta }t_{min}`$. Nevertheless, Eq. (6) provides us with a glimpse of how complicated it is to theoretically study the cluster-cluster collisions and coalescence in the random lattice model, and in general, in inelastic gases. If two clusters $`l`$ and $`l+1`$ are the first ones to coalesce (at time $`t_0+\mathrm{\Delta }t_1`$) in the random lattice model, then a new cluster with $`M_l+M_{l+1}`$ particles and center of mass velocity $`(M_lV_l+M_{l+1}V_{l+1})/(M_l+M_{l+1})`$ is formed. At time $`t_0+\mathrm{\Delta }t_1`$, the shape of the new cluster is different from a half-cosine, but as the mean-field theory of Sec. III.2.1 suggests, very soon the shape of the new cluster converges to a half-cosine, unless the newly formed cluster collides with another one in the meantime.

#### III.2.3 A “hydrodynamic” description of the random lattice model

So far, we have analyzed the system-wide properties of shocks and clusters. We have also seen that there exists an effective dynamics in terms of diffusion equations (4) and (5) within a single cluster. Based on this collected wisdom on the random lattice model so far, one can naturally ask if it is possible to express the system-wide properties of shocks and clusters in terms of an effective (“hydrodynamic”) equation.

It turns out that indeed such an equation can be constructed for the random lattice model:

$`{\displaystyle \frac{v}{t}}=D(r)[\mathrm{\Theta }(^+v)^+v\mathrm{\Theta }(^{}v)^{}v],`$ (7)
where $`\mathrm{\Theta }`$ denotes a unit step function, $`(^+v)_k=v_{k+1}v_k`$ and $`(^{}v)_k=v_kv_{k1}`$ are the discrete gradients operating to the right and left respectively and $`{\displaystyle \frac{}{t}}`$ the discrete time gradient. It is important to stress that in a numerical implementation of Eq. (7), all pairs of particles with negative relative velocities collide simultaneously in a unit timestep. The main difference between the numerical implementation of Eq. (7) and the random lattice model lies in the order of the collisions, but in the light of the previous results, we expect to find the shocks at the same locations. Moreover, inside a given cluster, (7) reduces to (4), so that the cluster velocity profiles should be the same as in the random lattice model.

Indeed, we find that the macroscopic velocity profile obtained numerically from Eq. (7) compares very well with the actual computer simulation results of the random lattice model (see Fig. 7), as expected.

## IV Discussion

In this paper, we have extensively studied structure formation and their dynamics in a one-dimensional inelastic lattice model for granular gases. These structures appear in the form of clusters separated by shocks. The locations of the shocks in this lattice model at early times are decided from the initial configurations of the particles’ velocities in a very robust manner. In order to predict the locations of the shocks, we have demonstrated a procedure to process the initial velocity configurations data of the particles, and this procedure works equally well also for the one-dimensional inelastic gas Frach1 ; BN . The coefficient of restitution does play a role to decide which of the large jumps in the particles’ velocities at early times yield shocks, but at late times, the macroscopic velocity profile for a given model with a given initial configuration of the particles’ velocities is independent of the coefficient of restitution.

Thus, we observe that in terms of detailed structure formation, systems with any $`r<1`$ “flow” towards the sticky limit $`r=0`$. Such “universality” has been found in Ref. BN  in terms of global quantities, while our results suggest that a broader “universality” holds for microscopic quantities, such as the locations of shocks (when the same initial velocity configuration is considered).

In addition, for the random lattice model, we have also studied the dynamics of an isolated cluster and cluster-cluster collisions in detail and obtained an effective “hydrodynamic” equation. We hope that the analyses presented here can be successfully used to study coarsening problems in realistic granular gases (i.e., in two or three dimensions).

S.O. and D.P. are financially supported by the Dutch research organization FOM (Fundamenteel Onderzoek der Materie).


format: text

text:
# Fluid-like behavior of a one-dimensional granular gas

## I Introduction

Scientists and engineers have been studying granular materials for nearly two centuries for their relevance both in natural processes (landslides, dunes, Saturn rings) and in industry (handling of cereals and minerals, fabrication of pharmaceuticals etc.). general1 ; general2 ; general3 ; general4  The understanding of the “granular state” still represents an open challenge and one of the most active research topics in non-equilibrium statistical mechanics and fluid dynamics. For instance, a way to attack the problem consists in fluidizing the grains by shaking them so that the system behaves as a non-ideal gas, a problem relatively easier to study. The difficulty, but also the beauty, of the dynamics of granular gases, meant as rarefied assemblies of macroscopic particles, stems from the inelastic nature of their collisions which leads to a variety of very peculiar phenomena. Several theoretical methods have been employed to deal with granular gases ranging from hydrodynamic equations, kinetic theories to molecular dynamics. Engineers often prefer the strategy of the continuum description because it gives a better grasp of real life phenomena, while natural scientists tend to opt for a microscopic approach, to better control each step of the modelization. The latter, as far as the interaction between particles is concerned, regards granular systems as peculiar fluids, and treats them through the same methods which have been successfully applied to ordinary fluids. hansen  This allows not only, to employ concepts already developed by physicists and chemists, but also to stress analogies and substantial differences.

The purpose of the present paper is to establish such a connection for a system of stochastically driven inelastic hard-rods constrained to move on a ring. The elastic version of this system has a long tradition lord ; tonks ; takahashi ; frenkel ; gursey ; salsburg ; percus  and is particularly suitable to test approximations and theories since many of its equilibrium properties can be derived in a closed analytical form. Even though, the one dimensional geometry introduces some peculiarities not shared by real fluids, we shall show that the model provides many useful information and a very rich phenomenology which closely recalls the behavior of microscopic particles confined in tubules or cylindrical pores with little interconnection. A second reason to investigate such a model is to show how the inelasticity of interactions influences not only the average global properties of a system, but also its microscopic local structure.

A basic requirement to a theoretical description of a granular gas is to provide an equation of state linking the relevant control parameters and possibly to relate it to the microscopic structure of the system. This connection is well known for classical fluids, where thermodynamic and transport properties are linked to the microscopic level via the correlation function formalism.

One dimensional models have been employed by several authors as simple models of granular gases. mcnamara ; sela ; grossman ; zhou ; zhou2 ; kadanoff ; mackintosh ; puglisi ; puglisi2 ; bennaim ; cordero ; baldassarri  The differences between the various models stem chiefly from the choice of the thermalizing device. In fact, granular gases would come to rest unless supplying energy compensating the losses due to the inelastic collisions. We call, by analogy, “heat bath” the external driving mechanism maintaining the system in a statistically steady state.

For the history, the first 1D models, which were proposed, had no periodic boundary conditions and the energy was injected by a vibrating wall (stochastic or not). This kind of external driving however, was not able to keep the system homogeneous, because only the first and last particle had a direct interaction with the wall. kadanoff  As an alternative, a uniform heating mechanism, namely a Gaussian white noise acting on each particle, was introduced. mackintosh  Later Puglisi et al.puglisi  added a second ingredient, consisting of a friction term that prevents the kinetic energy from diverging. With such a modification the system reaches a steady regime and time averages can be safely computed. In the present paper we shall focus on this last model characterizing its steady state properties.

The layout is the following: in section II we introduce the model, in section III we obtain numerically and by approximate analytical arguments equations for the average kinetic energy and pressure. Section IV is devoted to fluctuations of the system observables around their average values. In section V we study the diffusion properties of the system. Finally, in section VI we present the conclusions.

## II The model

Inelastic hard sphere models are perhaps the simplest models able to capture the two salient features of granular fluids, namely the hard core repulsion between grains and the dissipation of kinetic energy due to the inelastic collisions. Since many of the equilibrium properties of the 1D elastic hard rods are known in closed analytical form, such a system represents an excellent reference model even for the inelastic case. Let us consider $`N`$ identical impenetrable rods, of coordinates $`x_i(t)`$, mass $`m`$ and size $`\sigma `$, constrained to move along a line of length $`L`$. Periodic boundary conditions are assumed. The hard-core character of the repulsive forces among particles reduces the interactions to single binary, instantaneous collision events occurring whenever two consecutive rods reach a distance $`d_i(t)=x_i(t)x_{i1}(t)`$ equal to their length $`\sigma `$. When two inelastic hard-rods collide, their post-collisional velocities (primed symbols) are related to pre-collisional velocities (unprimed symbols) through the rule:

$$v_i^{}=v_i\frac{1+r}{2}(v_iv_j)$$
(1)
where $`r`$ indicates the coefficient of restitution. The interaction of each particle with the heat-bath is represented by the combination of a viscous force proportional to the velocity and a stochastic force. Then each particle follows the so called Kramers dynamics

$`{\displaystyle \frac{dx_i}{dt}}`$ $`=`$ $`v_i`$ (2)
$`m{\displaystyle \frac{dv_i}{dt}}`$ $`=`$ $`m\gamma v_i+\xi _i(t)`$ (3)
where $`\gamma `$ is the viscous friction coefficient, $`\xi _i(t)`$ is a Gaussian white noise with zero average and correlation

$$\xi _i(t)\xi _j(s)=2\gamma mT_b\delta _{ij}\delta (ts),$$
(4)
$`T_b`$ is the “heat-bath temperature” and $``$ indicates the average over a statistical ensemble of noise realizations.

We have developed a numerical simulation code for hard rods interacting through momentum conserving but energy dissipating collisions. In our simulations the motion between two consecutive collisions is governed by the dynamics (2,3). Thus, we determine the instant when the first collision among the $`N`$ particles occurs and change their velocities and positions according to the equations of motion. The effect of the collision is taken into account by updating the velocities after each collision according to the rule (1).

We tested our code on the elastic case ($`r=1`$) and checked that our simulations faithfully reproduced the well known properties of the equilibrium hard rod system.

## III Average properties

We begin by considering the steady state properties of the model. The aim is to derive relations connecting the microscopic parameters to the “thermodynamic” observables such as temperature and pressure and eventually to obtain an “equation of state” relating these two quantities. In order to achieve this goal we assume that the system is homogeneous, so that its density $`\rho `$ is constant.

### III.1 Kinetic temperature

Collisions and the Kramers’ dynamics entail that the time derivative of the average kinetic energy per particle is

$$\frac{d}{dt}\frac{1}{2}mv(t)^2=\gamma (T_bT_g)w(t),$$
(5)
where $`T_g=m_{i=1}^Nv_i^2/N`$ is the granular temperature and $`w(t)`$ is the average power dissipated by collisions, given by $`w=\frac{1r^2}{8}m\delta v^2/\tau _c`$, where $`\delta v`$ is the difference between the pre-collisional velocities of the colliding pair. The average collision time $`\tau _c`$ is estimated by assuming a mean free path $`\lambda =(L\sigma N)/N`$, where $`L\sigma N`$ is the free volume. We obtain, in terms of the system density $`\rho =N/L`$,

$$\tau _c=\frac{\lambda }{v}=\frac{1\sigma \rho }{\rho }\sqrt{\frac{m}{T_g}}.$$
(6)
Thus the average power dissipated per grain reads:

$$w=\frac{1r^2}{8}\delta v^2\frac{\rho \sqrt{T_gm}}{1\sigma \rho }.$$
(7)
In order to estimate $`T_g`$ we assume that $`m\delta v^24T_g`$ since the velocities of the colliding pairs are strongly correlated. Thus imposing the solution of Eq. (5) to be stationary we obtain for $`T_g`$ the following expression:

$$T_g=\frac{T_b}{1+\frac{1r^2}{2\gamma }\frac{\rho }{1\rho \sigma }\sqrt{\frac{T_g}{m}}}.$$
(8)
In figure 1 we compare formula (8) with the results of numerical simulations at various densities. In spite of the simplicity of the argument used to derive Eq. (8), the agreement between $`T_g`$ extracted by simulations and its theoretical estimate is fairly good.

### III.2 Kinetic pressure

In a granular system the total pressure, $`P`$, can be obtained via its mechanical or kinetic definition, i.e. as the impulse transferred across a surface in the unit of time.allen ; barrat  The pressure contains both the ideal gas and the collisional contribution ($`P_{id}`$ and $`P_{exc}`$ respectively)

$$P=P_{id}+P_{exc}=\rho T_g(\rho )+\frac{\sigma }{Lt_{ob}}\underset{k=1}{\overset{M_c}{}}\delta p_k$$
(9)
where the second equality stems from the virial theorem. erpenbeck  Here $`t_{ob}`$ is the observation time, the sum runs over the $`M_c`$ collisions and $`\delta p_k=m\delta v_k`$ represents the impulse variation due to the $`k`$-th collision.

An approximate formula for $`P_{exc}`$ can be derived as follows. The average collision frequency per particle can be estimated as $`\tau _c^1=(M_c/t_{ob})/N`$. By replacing in Eq. (9) $`t_{ob}`$ with $`(M_c/N)\tau _c`$ and using $`\tau _c`$ given by Eq. (6), we obtain, for the excess part of the pressure, the expression

$$P_{exc}=\frac{N\sigma }{\tau _cL}m\delta v_c=\frac{\rho ^2\sigma }{1\sigma \rho }T_g(\rho ).$$
(10)
Collecting pieces together we arrive at

$$P(\rho )=T_g(\rho )\left[\rho +\frac{\sigma \rho ^2}{1\sigma \rho }\right]=T_g(\rho )\frac{\rho }{1\sigma \rho }$$
(11)
which reproduces the well known Tonks formula tonks  in the case of elastic particles and constitutes the sought equation of state for the inelastic system. Let us recall that in the elastic case, Eq. (11) can be written in the virial form

$$P(\rho )=T_g(\rho )\rho [1+\rho \sigma g(\sigma )]$$
(12)
showing the connection between the macroscopic and the microscopic level, since $`g(\sigma )=1/(1\rho \sigma )`$ is the equilibrium pair correlation at contact.

We see, from figure (2), that the presence of the prefactor $`T_g(\rho )`$, which is decreasing function of the density, makes $`P(\rho )`$ to increase more slowly than the corresponding pressure of the elastic system in the same physical conditions (i.e. same density and contact with the same heat bath).

Equations (8) and (10) for temperature and pressure coincide, in the limit $`\gamma 0`$, $`\sigma 0`$, $`\gamma T_b=\mathrm{\Omega }=\text{const}`$, with those derived by Williams and MacKintosh mackintosh .

Following the standard approach to fluids, we define, even for the inelastic system, the response of the density to a uniform change of the pressure for a fixed value of the heat bath temperature:

$$\chi _T=\frac{1}{\rho }\frac{\rho }{P}$$
(13)
which is plotted in figure 3. We observe that the response of the inelastic system to a compression is much larger than the corresponding elastic system at the same density, due to the tendency to cluster.

## IV Fluctuations

So far we have considered only the global average properties of the granular gas. It is well known, on the other hand, that these system may exhibit strong spontaneous deviations from their uniform state. In this section we shall study fluctuations of the main observables in order to understand the qualitative effect of inelasticity on such a peculiar fluid.

### IV.1 Velocity distributions

One of the signatures of the inelasticity of the collisions is represented by the shape of the velocity distribution function (VDF), $`P(v)`$. Non-Gaussian VDF’s, displaying low velocity and high velocity overpopulated regions, have been measured experimentally losert ; menon ; Blair ; urbach ; aranson ; kudrolli  and in numerical simulations.puglisi ; breyvel  In figure 4 we show two VDF corresponding to two different values of $`r`$.

Theoretical, numerical and experimental studies have shown that the VDF for inelastic ($`r<1`$) gases usually displays overpopulated tails. The literature seems to indicate the lack of a universal VDF: in $`d>1`$ the solution of the homogeneous Boltzmann equation with inelastic collisions (with a stochastic driving similar to ours but without viscosity) has overpopulated tails of the kind $`\mathrm{exp}(v^{3/2})`$.ernst

### IV.2 Energy fluctuations

Interestingly, the energy fluctuations of our system, $`Em_iv_i^2`$/2, display a scaling with respect to the number of particles. We are interested in the quantity $`(E^2E^2)/E^2`$ as a function of $`N`$ (at fixed density $`\rho `$). Defining $`v^n=v^nP(v)𝑑v`$, we have

$$E^2E^2=\frac{m^2}{4}\underset{i,j}{}\left(v_i^2v_j^2v_i^2v_i^2\right)$$
(14)
Under the hypothesis that the variables $`v_i`$ are independently distributed, we get

$$E^2E^2=\frac{m^2}{4}[Nv^4Nv^2^2]$$
(15)
Since for Gaussian variables $`v^4=3v^2^2`$, we find

$$\frac{E^2E^2}{E^2}=\frac{2}{N}$$
(16)
which is a well known formula for equilibrium systems huang . This scaling is fairly well verified in figure 5.

In the case of a granular fluid, $`P(v_i)`$ is no more Gaussian, exhibiting fatter tails, so one observes $`v^43v^2^2`$ (for example in the case $`r=0.6`$ we have $`v^4/v^23.3`$). This leads to the conclusion that the scaling $`1/N`$ of formula (16) still holds, but with a coefficient larger than $`2`$. Simulation runs for $`r=0.6`$ confirm this prediction (Fig. 5). The renormalization of the multiplicative constant occurring in the inelastic system could be interpreted also as an “effective reduction” of the number of degrees of freedom. Indeed, the inelastic system has the tendency to cluster, as it will be shown, and therefore the effective number of independent “particles” appears smaller. Another appealing interpretation is that the inelastic systems possesses an effective “specific heat” larger than that of elastic systems.

### IV.3 Velocity correlations

A universal signature of the inelasticity is the presence of correlations between the velocities of the particles. We measured the structure function of the velocities $`v_i`$, $`S_v(k)=\stackrel{~}{v}(k)\stackrel{~}{v}(k)`$, where $`\stackrel{~}{v}(k)`$ is the Fourier transform of $`v_i`$. In figure 6, we show three $`S_v(k)`$ corresponding to the elastic ($`r=1`$) and inelastic system with $`r=0.8`$ and $`r=0.6`$.

As mentioned above, the elastic systems is characterized by uncorrelated velocities and this reflects on a constant structure function. A certain degree of correlation is instead evident in the inelastic system. In fact, the inelasticity reduces by a factor $`r`$ the relative velocity of two colliding particles and this leads to an increasing correlation among velocities. However, the noise induced by the bath competes with these correlations, making the structure function not very steep. More specifically, $`S_v(k)`$ can be fitted, in the middle range of $`k`$ values, by an inverse power $`k^{0.5}`$, while at high $`k`$ values it reaches a constant plateau. This is the fingerprint of a persistent internal noise (velocity fluctuations are not completely frozen by inelastic collisions).Ernstprogram

### IV.4 Distribution of interparticle spacing and of collision times

The probability distribution, $`P(\delta x)`$, of distances between nearest neighbor particles $`\delta x=x_ix_{i1}`$, shown in figure 7, provides information about the spatial arrangement of the system. In the elastic case one easily finds

$$P(\delta x)=\frac{1}{\lambda }\mathrm{exp}[(\delta x\sigma )/\lambda ]$$
(17)
for $`\delta x\sigma `$ and $`0`$ for $`\delta x<\sigma `$ with $`\lambda =(1\rho \sigma )/\rho `$. The presence of inelasticity modifies such a simple exponential law in the way shown in figure 7. In this case, the probability of finding two particles at small separation increases together with that of finding large voids. Such a picture is consistent with the idea of the clustering phenomenon:zanetti  two particles, after the inelastic collision, have a smaller relative velocity and therefore reach smaller distances, eventually producing dense clusters and leaving larger empty regions with respect to the elastic case.

On the contrary, the probability distribution of collision times, shown in figure 8, appears to always follow the theoretical (elastic) form $`P(t)=1/\tau _c\mathrm{exp}(t/\tau _c)`$. Apart from a trivial rescaling due to the change of the thermal velocity with $`r`$, it seems not to depend appreciably on the coefficient of restitution.

Such a finding is in contrast with the situation observed in 2D vibrated granular systems Paolotti ; Blair . A possible explanation for this discrepancy is the following: in the inelastic 1D system, there is correlation between the relative velocities and the free-paths (or free-times), otherwise the distribution $`P(t)`$ and $`P(\delta x)`$ would have had the same shape due to the trivial relation $`x=vt`$, $`v`$ being the the average velocity of the rods. In particular, the fact that the peak of $`P(\delta x)`$ in $`\delta x=0`$ does not yield a corresponding peak in the $`t=0`$ region of $`P(t)`$ suggests that the shorter the distance between particles the smaller their relative velocity.

### IV.5 Density fluctuations

We turn, now, to the study of the structural properties of the inelastic hard-rod gas, by considering the pair correlation function or the structure factor. As mentioned in section II, the virial equation (12) relates the pressure of an elastic hard-rod system to its microscopic structure. In the presence of inelasticity, however, we expect that the tendency to cluster is mirrored by a change in the structural properties of the fluid. Therefore we considered the behavior of the static (truly speaking steady state) structure factor

$$S(k)=1+\rho 𝑑x[g(x)1]e^{ikx}$$
(18)
for different values of $`\rho \sigma `$ and inelasticity.

The spatial structure of the system is determined, as in ordinary fluids, by the strong repulsive forces. Their role is seen in the oscillating structure of $`g(x)`$. The inelastic nature of the collisions provides a correction to $`g(x)`$, which can be better appreciated by studying the small wavelength behavior of $`S(k)`$ which develops a peak at small $`k`$ recalling the Ornstein-Zernike behavior

$$S(k)\frac{1}{S^1(0)+c_2k^2}.$$
(19)
The coefficient $`c_2`$ is negative for hard rods, whereas it is positive for the inelastic system. For hard rods, $`S(k)`$ is known percus  and reads:

$$S(k)=\frac{1}{1+2b\rho \sigma \frac{\mathrm{sin}(k\sigma )}{k\sigma }+(b\rho \sigma )^2\frac{\mathrm{sin}^2(k\sigma /2)}{(k\sigma /2)^2}}$$
(20)
with $`b=\frac{1}{1\rho \sigma }`$.

Figure 9 shows the typical behaviors of $`S(k)`$ for elastic and inelastic systems. The numerically computed structure factor of the elastic system agrees rather well with equation (20). The inelastic system, instead, displays a peak in the small $`k`$ region reflecting the tendency of the fluid to cluster. The peak increases with inelasticity, demonstrating that the energy dissipation in collisions is responsible for these long range correlations. Incidentally, we comment that such a behavior of $`S(k)`$ could be attributed to the presence of a long range attractive effective potential between the rods, as a result of dynamical correlations.ernst2

MackIntosh and Williams mackintosh  found that, in the case of randomly kicked rods in the absence of viscosity, the pair correlation function decays as an inverse power law, $`g(x)x^\eta `$ with $`\eta 0`$ for $`r1`$ and $`\eta 1/2`$ for $`r0`$. Correspondingly, one expects $`S(k)`$ to diverge as $`k^{1/2}`$ as $`k0`$, for very inelastic systems. In other words, the inelasticity leads to long range spatial correlations which are revealed by the peak at small $`k`$ of $`S(k)`$. We remark that, in spite of the apparent similarity between the equations of state for elastic and inelastic system, their structural properties are radically different. Such a phenomenon is the result of the coupling of the long-wavelength modes of the velocity field with the stochastic non-conserved driving force. In fact, due to the inelastic collisions, the velocities of the particles tend to align, thus reducing the energy dissipation. On the other hand, these modes adsorb energy from the heat bath and grow in amplitude, and only the presence of friction prevents these excitations from becoming unstable. The density field, which is coupled to the velocity field by the continuity equation, also develops long range correlations, and the structure factor displays a peak at small wave-vectors.

## V Transport properties

One dimensional hard-core fluids exhibit an interesting connection between the microscopic structural properties and diffusive ones. In the present section, we present numerical results for the collective diffusion and for the diffusion of a tagged particle and show how these are connected to the structure.

### V.1 Collective Diffusion and Self-Diffusion

Let us turn to analyze the perhaps simplest transport property of the hard rods system, namely the self-diffusion, i.e. the dynamics of a grain in the presence of $`N1`$ partners. The problem is highly non-trivial since the single grain degrees of freedom are coupled to those of the remaining grains. Such a single-filing diffusion is also relevant in the study of transport of particles in narrow pores. kollmann

The diffusing particles can never pass each other. The excluded volume effect represents a severe hindrance for the particles to diffuse. In fact, a given particle in order to move must wait for a collective rearrangement of the entire system. Only when the cage of a particle expands, the tagged particle is free to diffuse further. This is a peculiar form of the so called cage effect which is enhanced by the one-dimensional geometry. In addition, the cage effect produces a negative region and a slow tail in the velocity autocorrelation function.

As an appropriate measure of the self-diffusion, we consider the average square displacement of each particle from its position at a certain time, that we assume to be $`t=0`$ without loss of generality

$$R(t)=\frac{_{i=1}^N[x_i(t)x_i(0)]^2}{N}$$
(21)
At an early stage, the self-diffusion is expected to display ballistic behavior, $`R(t)|v|^2t^2`$, with $`|v|^2=T_g/m`$, before any perturbation (heat bath and collisions) change the free motion of particles, i.e. when $`t1/\gamma `$ and $`t\tau _c`$.

A system of non-interacting (i.e. non-colliding) particles subjected to Kramers’ dynamics (2,3) displays, after the ballistic transient, normal self-diffusion of the form $`R(t)2D_0t`$ with $`D_0=T_g/\gamma `$. This is well verified in figure 10 (circles).

Lebowitz and Percus percuslebow  studied the tagged particle diffusion problem for systems governed by non-dissipative dynamics without heat bath and found a diffusive behavior described by

$$R(t)=2D_{coll}t=\lambda |v|t$$
(22)
with $`\lambda =(1\rho \sigma )/\rho `$. Since $`|v|=\sqrt{T/2\pi m}`$ one obtains

$$D_{coll}=\frac{1\rho \sigma }{\rho }\sqrt{\frac{T}{2\pi m}}$$
(23)
On the other hand, almost in the same years, Harris harris  studied the behavior of $`R(t)`$ in the case of $`N`$ identical Brownian particles with hard-core interactions , i.e. obeying a single-filing condition, and obtained a sub-diffusive behavior increasing as

$$R(t)=2\lambda \left(\frac{D_0t}{\pi }\right)^{1/2}$$
(24)
where $`D_0`$ is the single (noninteracting) particle diffusion coefficient.

In figure 10, we study the self-diffusion $`R(t)`$ for elastic and inelastic systems in the presence of heat bath and viscosity, obtaining two different regimes separated by a typical time $`\tau _c`$. In the first transient regime $`t<\tau _c`$ we we observe the ballistic motion. In the second stage, instead, we expect the sub-diffusive behavior, $`R(t)t^{1/2}`$, predicted by Harris and other authors. harris ; pincus ; olandesi . The inelastic system displays the same sub-diffusive behavior, but with a multiplicative constant larger than $`1`$, i.e. at equal times the granular (inelastic) fluid has a larger absolute value of $`R(t)`$.

It is interesting to analyze the connection between this transport property and the compressibility of the system, as remarked by Kollmann. kollmann

### V.2 Connection between self-diffusion and structure

We follow Alexander and Pincus pincus  argument in order to derive a formula for the self-diffusion and show the connection with the compressibility of the system. Let us consider the two time correlator:

$$\rho _k(t)\rho _k(0)=\underset{ij}{}e^{ikx_i(t)}e^{ikx_j(0)}$$
(25)
we set $`x_i(t)=X_i+(x_i(t)X_i)=X_i+u_i`$, where $`X_i`$ are the nodes of the 1D lattice $`X_i=ia`$, ($`a`$ being the lattice spacing). Expanding around $`X_i`$

$$\rho _k(t)\rho _k(0)k^2\underset{ij}{}\text{e}^{ik(X_iX_j)}u_i(0)u_j(t)$$
(26)
hence

$$\widehat{u}_k(t)\widehat{u}_k(0)=\frac{\rho _k(t)\rho _k(0)}{Nk^2}$$
(27)
We assume, now, that the density correlator varies as

$$\rho _k(t)\rho _k(0)=\rho _k(0)\rho _k(0)e^{Dk^2t}$$
(28)
where $`D=D_0/S(0)`$ (this is demonstrated in Appendix A) is the collective diffusion coefficient, with $`D_0=T_b/\gamma `$. Now, the mean square displacement per particle (21) can be written as $`R(t)=_l[u_l(t)u_l(0)]^2/N`$ and, in Fourier components, reads

$$R(t)=\frac{2}{N}\underset{k}{}\widehat{u}_k(0)\widehat{u}_k(0)\widehat{u}_k(t)\widehat{u}_k(0)$$
(29)
employing Eq. 27 we find

$$R(t)=2\underset{k}{}\rho _k(0)\rho _k(0)\frac{1\text{e}^{Dk^2t}}{N^2k^2}$$
(30)
Approximating the sum with an integral ($`_kL/(2\pi )_{\pi /a}^{\pi /a}𝑑k`$ and recalling that

$$NS(k)=\rho _k(0)\rho _k(0)$$
(31)
we obtain

$$R(t)=\frac{2L}{N}_0^{\pi /a}\frac{dk}{2\pi }S(k)\frac{1\text{e}^{Dk^2t}}{k^2}$$
(32)
and therefore

$$R(t)\frac{2}{\rho \pi }S(0)\sqrt{\pi Dt}=\frac{2}{\rho }\sqrt{\frac{D_0S(0)t}{\pi }}$$
(33)
Notice that such a formula in the case of hard-rods is identical to formula (24).

We see that the tagged particle diffusion depends on the structure of the fluid. In the granular fluid the $`k0`$ part of the spectrum is enhanced and thus we expect a stronger tagged particle diffusion. This is what we observe. Physically there are larger voids and particles can move more freely. Let us notice that as far as the collective diffusion is involved the spread of a group of particles is faster in the presence of repulsive interactions than without. tarazona

## VI Conclusions

In this paper we have studied a one-dimensional system of inelastic hard-rods coupled to a stochastic heat-bath with the idea that it can represent a reference system in the area of granular gases to test theories and approximations. Due to the relative simplicity of the one dimensional geometry we have shown that it is possible to obtain relations between the macroscopic control parameters such as kinetic temperature, pressure and density. We tested these analytical predictions against the numerical measurements and found a fairly good agreement. It also appears that many properties of the heated one-dimensional inelastic hard rod system are similar to those of ordinary fluids. However, when we have considered how various physical observables fluctuate about their equilibrium values, many relevant differences have emerged. These range from the non-Gaussian behavior of the velocity distribution, the peculiar form of the distribution of distances between particles and of the energy fluctuations to the shape of the structure factor at small wave-vectors. Finally, we have found that the diffusive properties of the system are affected by the inelasticity and in particular the self-diffusion is enhanced.

To conclude, in spite of the similarity between ordinary fluids and granular fluids, which has been recognized for many years and has made possible to formulate hydrodynamical equations for granular media in rapid, dilute flow, the presence of anomalous fluctuations in the inelastic case indicates the necessity of a treatment which incorporates in a proper way both the local effects such as the excluded volume constraint and the long ranged velocity and density correlations. Such a program has been partially carried out by Ernst and coworkers Ernstprogram , but needs to be completed regarding the description of the fluid structure.

## VII Appendix A

In the case of over-damped dynamics (i.e. large values of $`\gamma `$) one finds that the collective diffusion is given by tarazona

$$\frac{\rho (x,t)}{t}=\frac{1}{\gamma }\frac{}{x}\left\{\rho (x,t)\frac{\mu (\rho (x))}{x}\right\}$$
(34)
where $`\mu (\rho (x))`$ is the local chemical potential. Expanding $`\mu `$ about its average value $`\rho _0`$ we obtain:

$$\frac{\mu (\rho (x))}{x}=\left[\frac{\delta \mu }{\delta \rho }\right]_{\rho _0}\frac{\rho (x)}{x}$$
(35)
Substituting into Eq. (34) we find

$$\frac{\rho (x,t)}{t}=\frac{1}{\gamma }\rho _0\left[\frac{\delta \mu }{\delta \rho }\right]_{\rho _0}\frac{^2\rho (x)}{x^2}$$
(36)
and using

$$S(0)=\frac{K_BT_b}{\rho _0}\left[\frac{\rho _0}{\mu }\right]_T$$
(37)
in the case of elastic hard rods we obtain:

$$\frac{\rho (x,t)}{t}=\frac{1}{\gamma }\frac{K_BT_b}{S(0)}\frac{^2\rho (x)}{x^2}=\frac{D_0}{S(0)}\frac{^2\rho (x)}{x^2}$$
(38)
Thus the renormalized diffusion coefficient is $`D=D_0/S(0)`$.


text:
# Spin-filter magnetoresistance in magnetic barrier junctions

## I Introduction

Since the observation of tunnel magnetoresistance (TMR) in magnetic tunnel junctions , there has been a lot of interest in structures consisting of two ferromagnetic (FM) electrodes separated by one insulator or semiconductor (FM/I(S)/FM). This is due to the possible applications as magnetic field sensors and memory cells in magnetic random access memories . The TMR in the magnetic tunnel junctions is limited by the spin polarizations of the FM electrodes. Recently, large magnetoresistance ratios as high as 30% were reported in Fe/Al2O3/Fe  and CoFe/Al2O3/Co .

The TMR of such FM/I(S)/FM tunnel junctions can be understood in terms of a two-band model in which the energy bands of the FM electrodes are split into spin-up and spin-down bands with different density of states at the Fermi energy. When the magnetization of the electrodes is parallel, the spin-up (spin-down) electrons tunnel from a majority (minority) spin to a majority (minority) spin band, whereas in the antiparallel alignment, the spin-up (spin-down) electrons are forced to tunnel from a majority (minority) spin to a minority (majority) spin band. This gives rise to a change in the conductance when the magnetizations are switched. During the last 10 years, a lot of theoretical papers were published on this topic (see reviews in Refs. ).

Using magnetic insulating layers, instead of having FM electrodes, one can obtain very high values for the TMR . Tunneling experiments using ferromagnetic semiconductors (FMSs), in particular the Eu chalcogenides, as the magnetic tunnel barriers , have shown spin polarization in the tunnel current which in favorable cases exceeds 99% . When a magnetic barrier, which acts as spin filter, is used in a tunnel junction, due to the spin splitting of its conduction band below $`T_C`$ (Curie temperature), tunneling electrons see a spin dependent barrier height. In this case the probability of tunneling for one spin channel will be much larger than the other, and a highly spin polarized current may result. More recently, LeClair et al.  using a spin-filter barrier and an FM electrode, obtained a large TMR in an Al/EuS/Gd tunnel junction which is a new method for injecting spins into semiconductors. The spin filtering and the magnetoresistance have also been investigated theoretically in semimagnetic semiconductors .

In this paper, using the transfer matrix method, we study theoretically the effect of the thickness of the magnetic barrier and the applied bias on the TMR for tunneling through a FM/FMS/NM tunnel junction with an asymmetrical barrier. We assume that the electron wave vector parallel to the interfaces and the spin direction of the electron are conserved in the tunneling process through the whole system.

In section 2, we describe the model and present a general formula for the tunneling current through the magnetic tunnel junction. Numerical results for a typical tunnel junction are presented in section 3. A brief summary is given in section 4.

## II Model and formalism

Consider a FM/FMS/NM sandwiched structure in the presence of DC bias $`V_a`$ as shown in Fig. 1, where $`\varphi _L`$ and $`\varphi _R`$ are the barrier heights in the left- and right-hand side of the FMS layer above $`T_C`$, respectively. For a tunnel junction with different electrode materials, the difference in barrier height at both metal-insulator interfaces makes the barrier oblique, resulting in an asymmetric current-voltage behavior. In this case, an asymmetry parameter $`\mathrm{\Delta }\varphi =\varphi _L\varphi _R`$ is introduced at zero bias to account for the tilted barrier potential . Here we consider the case that the tunneling electron with energy $`E_x`$ is incident from the left and transmits to the right along the $`x`$ direction. In a free-electron approximation of the spin-polarized conduction electrons, the longitudinal part of the effective one-electron Hamiltonian may be written as

$$H_x=\frac{\mathrm{}^2}{2m_j^{}}\frac{d^2}{dx^2}+U_j(x)+V^\sigma 𝐡(x)\sigma ,$$
(1)
where $`m_j^{}`$ ($`j`$=1-3) are the electron effective masses in three regions labeled in Fig. 1, and

$$U_j(x)=\{\begin{array}{cc}0,& x<0,\\ E_{FL}+\varphi _L(\mathrm{\Delta }\varphi +eV_a)x/d,& 0x<d,\\ eV_a,& xd,\end{array}$$
(2)
where $`E_{FL}`$ is the Fermi energy of the FM electrode and $`d`$ is the barrier width. $`V^\sigma `$ which is a spin-dependent potential, denotes the $`sf`$ exchange coupling between the spin of tunneling electrons and the localized $`f`$ spins in the magnetic barrier. This term, within the mean field approximation, is proportional to the thermal average of the $`f`$ spins, $`S_z`$ (a 7/2 Brillouin function), and can be written as $`V^\sigma =I\sigma S_z`$. Here, $`\sigma =\pm 1`$ which corresponds to $`\sigma =,`$, respectively and $`I`$ is the $`sf`$ exchange constant in the FMS layer. $`𝐡(x)\sigma `$ is the internal exchange energy where $`𝐡(x)`$ is the molecular field in the FM electrode and $`|𝐡|=h_0`$. Although the transverse momentum $`\mathrm{}𝐤_{}`$ is omitted from the above notations, the summation over $`𝐤_{}`$ is carried out in our calculations.

The Schrödinger equation for a biased barrier layer can be simplified by a coordinate transformation whose solution is the linear combination of the Airy function Ai\[$`\rho (x)`$\] and its complement Bi\[$`\rho (x)`$\] . Considering all three regions of the FM/FMS/NM junction shown in Fig. 1, the eigenfunctions of the Hamiltonian (1) with eigenvalue $`E_x`$ have the following forms

$$\psi _{j\sigma }(x)=\{\begin{array}{cc}A_{1\sigma }e^{ik_{1\sigma }x}+B_{1\sigma }e^{ik_{1\sigma }x},& x<0,\\ A_{2\sigma }\text{Ai}[\rho _\sigma (x)]+B_{2\sigma }\text{Bi}[\rho _\sigma (x)],& 0x<d,\\ A_{3\sigma }e^{ik_3x}+B_{3\sigma }e^{ik_3x},& xd,\end{array}$$
(3)
where

$$k_{1\sigma }=\sqrt{2m_1^{}(E_x+h_0\sigma )}/\mathrm{},$$
(4)
and

$$k_3=\sqrt{2m_3^{}(E_x+eV_a+E_{FR}E_{FL})}/\mathrm{},$$
(5)
are the electron wave vectors along the $`x`$ axis. $`A_{j\sigma }`$ and $`B_{j\sigma }`$ are constants to be determined from the boundary conditions, while

$$\rho _\sigma (x)=\frac{x}{\lambda _0}+\beta _\sigma ,$$
(6)
with

$$\lambda _0=\left[\frac{\mathrm{}^2d}{2m_2^{}(\mathrm{\Delta }\varphi +eV_a)}\right]^{1/3},$$
(7)
$$\beta _\sigma =\frac{[E_xE_{FL}\varphi _LV^\sigma ]d}{(\mathrm{\Delta }\varphi +eV_a)\lambda _0}.$$
(8)
Upon applying the boundary conditions such that the wave functions and their first derivatives are matched at each interface point $`x_j`$, i.e., $`\psi _{j,\sigma }(x_j)=\psi _{j+1,\sigma }(x_j)`$ and $`(m_j^{})^1[d\psi _{j,\sigma }(x_j)/dx]=(m_{j+1}^{})^1[d\psi _{j+1,\sigma }(x_j)/dx]`$, we obtain a matrix formula that connects the coefficients $`A_{1\sigma }`$ and $`B_{1\sigma }`$ with the coefficients $`A_{3\sigma }`$ and $`B_{3\sigma }`$ as follows:

$`\left[\begin{array}{c}A_{1\sigma }\\ B_{1\sigma }\end{array}\right]=M_{total}\left[\begin{array}{c}A_{3\sigma }\\ B_{3\sigma }\end{array}\right],`$ (13)
where

$`M_{total}`$ $`=`$ $`{\displaystyle \frac{k_3}{k_{1\sigma }}}\left[\begin{array}{cc}ik_{1\sigma }& \frac{1}{\lambda _0}\frac{m_1^{}}{m_2^{}}\\ ik_{1\sigma }& \frac{1}{\lambda _0}\frac{m_1^{}}{m_2^{}}\end{array}\right]\left[\begin{array}{cc}\text{Ai}[\rho _\sigma (x=0)]& \text{Bi}[\rho _\sigma (x=0)]\\ \text{Ai}^{}[\rho _\sigma (x=0)]& \text{Bi}^{}[\rho _\sigma (x=0)]\end{array}\right]`$ (26)
$`\times \left[\begin{array}{cc}\text{Ai}[\rho _\sigma (x=d)]& \text{Bi}[\rho _\sigma (x=d)]\\ \text{Ai}^{}[\rho _\sigma (x=d)]& \text{Bi}^{}[\rho _\sigma (x=d)]\end{array}\right]^1\left[\begin{array}{cc}ik_3& \frac{1}{\lambda _0}\frac{m_3^{}}{m_2^{}}\\ ik_3& \frac{1}{\lambda _0}\frac{m_3^{}}{m_2^{}}\end{array}\right]^1`$
$`\times \left[\begin{array}{cc}e^{ik_3d}& 0\\ 0& e^{ik_3d}\end{array}\right]^1.`$
Since there is no reflection in region 3, the coefficient $`B_{3\sigma }`$ in Eq. (3) is zero and the transmission coefficient of the spin $`\sigma `$ electron which is defined as the ratio of the transmitted flux to the incident flux can be written as

$$T_\sigma (E_x,V_a)=\frac{k_3m_1^{}}{k_{1\sigma }m_3^{}}\left|\frac{1}{M_{total}^{11}}\right|^2,$$
(27)
where $`M_{total}^{11}`$ is the left-upper element of the matrix $`M_{total}`$ defined in Eq. (26).

At $`T=0`$ K, the spin-dependent current density for the magnetic tunnel junction in the free-electron model is given by the formula

$$J_\sigma =\frac{em_1^{}}{4\pi ^2\mathrm{}^3}\left[eV_a_{E_0^\sigma }^{E_FeV_a}T_\sigma (E_x,V_a)𝑑E_x+_{E_FeV_a}^{E_F}(E_FE_x)T_\sigma (E_x,V_a)𝑑E_x\right],$$
(28)
where $`E_0^\sigma `$ is the lowest possible energy that will allow transmission and is given by $`E_0^{}=`$max$`\{h_0,(eV_a+E_{FR}E_{FL})\}`$ for spin-up and $`E_0^{}=h_0`$ for spin-down electrons.

The tunnel conductance per unit area is given by $`G=_\sigma J_\sigma /V_a`$. In this case, the TMR can be described quantitatively by the relative conductance change as

$$\text{TMR}=\frac{G_{}G_{}}{G_{}},$$
(29)
where $`G_{}`$ and $`G_{}`$ correspond to the conductances in the parallel and antiparallel alignments of the magnetizations, respectively.

## III Numerical results

Taking the Fe/EuS/Al tunnel junction as an example, we calculate the tunnel current and TMR of the junction according to the Eqs. (27)-(29). We have chosen Fe and EuS because they have cubic structures and the lattice mismatch is only 4% . The appropriate parameters for EuS which have been used in this article are: $`S`$=7/2, $`I`$=0.1 eV . The parameters $`E_{FL}`$ and $`h_0`$ for Fe layer are taken corresponding to $`k_F`$=1.09 Å-1 and $`k_F`$=0.42 Å-1 (for itinerant $`d`$ electrons) . In the Al layer $`E_{FR}=11.7`$ eV . The barrier heights at the interfaces, which can be derived from the work functions  of Fe and Al and the electron affinity of EuS , are taken as $`\varphi _L`$=1.94 eV and $`\varphi _R`$=1.7 eV. In practice, the effective masses of electrons may differ from that of free electron, but here for simplicity, we assume all electrons have the same mass, $`m`$ as free electrons. We show the numerical results at $`T`$=0 K. Thus, in the case of parallel (antiparallel) alignment $`S_z=S`$ ($`S_z=S`$).

In Fig. 2, the TMR is shown as a function of the width of the magnetic barrier $`d`$, when the bias voltage $`V_a`$=0.5 V is applied to the junction. At zero temperature and for the nearly normal incidence, the electrons with $`E_x`$ near $`E_F`$ carry most of the current. Furthermore, the Fermi energy and the barrier heights at fixed temperature, are constant. Therefore, with increasing the thickness of the EuS layer from zero, when the thickness approaches half the wavelength of the electron wave in the magnetic barrier, the tunnel conductance with the parallel alignment increases faster than with the antiparallel alignment. In this case, we can expect a peak in the TMR. By increasing the thickness further, the tunnel conductance and hence the TMR decreases exponentially.

Figure 3 shows the spin-dependent current densities for spin-up and spin-down electrons as a function of applied bias. In the parallel alignment, the tunnel current for spin-up electrons is much higher than the spin-down ones, while in the antiparallel alignment, the discrepancy between the tunnel currents of the two spin channels increases very slowly. The origin of this effect is that the tunnel current for each spin channel depends on the density of states in the Fe layer and the heights of the tunnel barrier. In the parallel (antiparallel) alignment, the majority spin electrons tunnel through a barrier with low (high) height, whereas the minority spin electrons tunnel through a barrier with high (low) height, thus we can expect a high (low) spin-polarized current.

It is clearly observed that for both the parallel and antiparallel alignments at low voltages the tunnel current curves vary linearly. With increasing the bias voltage, the effective width of the barrier becomes narrower and hence the slope of the curves increases. This effect is very strong for spin-up electrons in the parallel alignment, and causes the TMR increases with increasing the bias voltage.

In order to reveal another aspect of the spin filtering phenomenon in the system, we have displayed the voltage dependence of the TMR in Fig. 4. As the figure shows for $`d`$=7 and 10 Å, with increasing the bias voltage the TMR increases linearly and decreases at high voltages. This decrease which is not identical for different barrier thicknesses is a consequence of the narrowing of the barrier width. At low voltages, the tunneling probability for one spin channel is higher than the other, whereas at high voltages, by narrowing the barrier, this tunneling probability increases for both channels; thus, the tunnel conductance ratio, $`G_{}/G_{}`$, tends toward unity, and hence the TMR decreases. This dependence of the TMR on bias voltage is not similar to the results of the FM/I(S)/FM tunnel junctions , because in the present study, we have a spin-dependent barrier. In FM/I(S)/FM junctions, both the majority and minority electrons see identical barrier heights. In this case, as the bias voltage increases, the imbalance between the number of majority and minority tunneling states diminishes and, therefore, the TMR decreases.

From these results one can see that, for the FM/FMS/NM tunnel junctions, the spin transport and hence the TMR can be controlled by the applied bias and the thickness of the FMS layer.

## IV Concluding remarks

In summary, we have shown how using a magnetic barrier and a FM electrode one leads to a new magnetoresistive device. The obtained results for spin-filter magnetoresistance show that the spin currents depends on the relative magnetization orientation of the FM electrode and the spin-filter barrier. It was also found that the low spin-polarized current will be amplified through the magnetic barrier; thus, a very large spin-polarized current from the FM electrode is injected into the NM electrode and a high value for the TMR is then obtained.

Our analysis of the TMR can be potentially useful to achieve larger TMR by optimally adjusting the material parameters. The two materials systems, ferrites and garnets  are both insulating and ferrimagnetic above room temperature. Therefore, they are well suited for obtaining large TMR in FM/FMS/NM. However, for successful spintronics applications, future efforts will have to concentrate on fabricating FM semiconductors in which ferromagnetism will persist at higher temperatures.

## Acknowledgments

I would like to express sincere thanks to Dr. Patrick R. LeClair for helpful discussion.