Dynamics of an Intruder in Dense Granular Fluids

We investigate the dynamics of an intruder pulled by a constant force in a dense two-dimensional granular fluid by means of event-driven molecular dynamics simulations. In a first step, we show how a propagating momentum front develops and compactifies the system when reflected by the boundaries. To be closer to recent experiments \cite{candelier2010journey,candelier2009creep}, we then add a frictional force acting on each particle, proportional to the particle's velocity. We show how to implement frictional motion in an event-driven simulation. This allows us to carry out extensive numerical simulations aiming at the dependence of the intruder's velocity on packing fraction and pulling force. We identify a linear relation for small and a nonlinear regime for high pulling forces and investigate the dependence of these regimes on granular temperature.


Introduction
The response of an intruder to a pulling force is a rather versatile tool to study the local nonequilibrium dynamics in various complex fluids, such as glasses [3], colloids [4,5] and granular media [1,2,6,7]. The measured force-velocity relations reveal striking nonlinear behaviour close to the glass and/or jamming transition. In the early experiments on colloids evidence was given for a threshold force even in the fluid regime, whereas mode-coupling theory predicts such behaviour only in the glassy phase. In the experiment of Ref. [1,2] two transitions are observed: the first one, called fluidization, separates a regime of continuous to intermittent motion of the intruder. It occurs below the jamming point (the second transition) and depends on the applied pulling force. Experimentally it is observed even for the smallest pulling force with the possibility of a dynamic transition inherent to the vibrated granular fluid for zero applied force. In Ref. [8], the dynamics of an intruder was simulated near the jamming point. One result of this study are velocity-force relations which are linear for packing fraction, η ≤ 0.833 and become nonlinear for η still closer to the jamming point.
In contrast to [8], we discuss a stochastically driven system, describing a fluidized granular state -similar to the experimental setup in [1,2]. A recent mode-coupling theory [9] predicts that such a "thermalized" granular fluid undergoes a glass transition at a packing fraction below the jamming point. This transition is different from both, the jamming transition at zero temperature and the glass transition for either Newtonian or Brownian dynamics.
In this paper we use event driven simulations to analyze the dynamics of an intruder in a two-dimensional system of hard disks close to the glass transition. We compute force-velocity correlations in the linear and nonlinear regime, extract the mobility for the linear regime and discuss scaling for the nonlinear regime. Moreover we discuss the dependence on the granular temperature.

Model
We consider a bidisperse system of N hard disks in an area A. Particles' positions and velocities are denoted by The size ratio R s /R b = 4 : 5 of small to big particles is chosen as in the experiments of Refs. [1,2]. The mass ratio follows from the disk like shape of the particles as m s /m b = 16/25. The intruder with position r 0 and velocity v 0 is chosen twice as large as the small particles, R 0 = 2R s and m 0 = 4m s .
The particles collide inelastically so that in each collision energy is dissipated while momentum is conserved. The simplest model ignores rotational degrees of freedom, allowing for normal restitution only. The collision rules for two colliding particles (say 1 and 2 ) is given in terms of their relative velocity g : where ǫ is the coefficient of restitution and n denotes the unit vector n := (r 1 − r 2 )/|r 1 − r 2 | at contact. In between collisions the particles perform Brownian motion, due to friction with a surrounding medium or with the bottom plate. We model this by a frictional force in the equation of motion: F fr = −γmv. For all finite values of the friction coefficient γ, momentum is lost, whereas for γ = 0, the total momentum is conserved.
Inelastic collisions give rise to a loss of energy ∝ To balance this dissipation the entire system is driven stochastically (like an air-fluidized bed [10,11,12] or a vibrating bottom plate [13,14,15]) by instantaneous kicks. The momentum of the i-th particle p i (t) is changed according to with ξ α i (t)ξ β j (t ′ ) = δ ij δ αβ δ(t−t ′ ) and zero mean ξ α i (t) = 0, and α = x, z denoting the cartesian components. In the frictional case, γ = 0, we kick single particles, whereas for the frictionless case, γ = 0, we kick neighbouring particles with equal and opposite momenta in order to conserve the total momentum locally [16].
We want to investigate the dynamics of a tagged particle under the action of a deterministic pulling force F. This so called intruder with position r 0 and velocity v 0 is subject to systematic kicks mimicking a constant force acting on the intruder. The frequency of these systematic kicks (∆t) −1 is chosen 4 orders of magnitude higher than the frequency of the stochastic driving and the collision frequency (see below). The systematic force does work on the system, injecting momentum and energy. Combining inelastic collisions, stochastic driving, pulling force and friction, we arrive at the following equation of motioṅ We aim to prepare the granular fluid without pulling force to be under stationary conditions. This can be achieved by balancing energy dissipation due to friction and inelastic collisions by the stochastic driving. Without external forcing the mean velocity of the particles vanishes, so that the granular temperature, i , is just the average kinetic energy of the particles. Its time rate of change can be deduced from Eq. 4 by multiplying with p α i (t) yielding 1 2 The explicit calculation of the time derivations on the rhs are beyond the scope of this paper but may be found elsewhere [17,18]. The time rate of change of the granular temperature then reads as follows: where the effective mass is given by m eff = msm b ms+m b . For simplicity, we choose the driving frequency f dr equal to the Enskog frequency ω E [19] and measure length and mass in units, such that R s = 1 and m s = 1. In the stationary state, dT dt = 0, the amplitude of the stochastic driving is given by: We prepare our system in a stationary state with T = 1, yielding a numerical value for p dr depending on γ, ω E and ǫ.

Implementation as event-driven simulation
In order to apply an event-driven algorithm to the dynamics, as described by Eq. (4), we need to generalize the code to account for damping (γ = 0). As usual, events include collisions of particles, wall collisions, subbox-wall collisions and (discrete) driving events (kicks). Standard event driven algorithms calculate the time of an upcoming event and advance all particles to that time under the assumption that the particle motion is ballistic. Hence, in between events, the velocities do not change. The condition for a collision is , the difference in the particles' trajectories at the collision time t • must be equal to the sum of their radii. Plugging in the trajectories subject to ballistic motion one finds a quadratic equation, that must have a positive solution for (t • − t 0 ) if and only if a collision will occur [20,21]. In our case, the velocities decrease as Inserting this into the condition for a collision, we find a similar condition for a collision except that (t Hence, collisions will occur at the same place and in the same order but at a different time and with different particles' velocities, provided the place of the collision is within reach of the particle. In detail, the maximum distance a particle is able to travel is limited because of the friction slowing the particle down. The maximal distance is given by r max that must be larger than the distance to the place at which the event will take place. Hence, (r i (t • ) − r i (t 0 )) < v i (t 0 )/γ or (t • − t 0 ) < 1/γ. This is consistent with Eq. (9) which requires log (1 − γ · (t • − t 0 )) to be negative for a positive collision time. Advancing a particle to a different type of event ( e.g. kicks) is done in the same way.
Changing an existing event-driven simulation in the way discussed enables us to simulate systems of hard disks and spheres subject to friction almost as fast as without friction. The systematic force on the intruder as described in Eq. 3 is implemented as frequent kicks on the intruder. In order to avoid an inelastic collapse, i.e., a infinite number of collisions in a finite time interval, we use the same method as described in [16] to circumvent it.

Equilibration and data aquisition
For packing fractions up to η = 0.8 we used the same data sets from [22] as initial conditions. For even denser systems, we used a compactified system acquired as described in section 4.1.
The data shown in section 4.2 were obtained by first equilibrating 100 different configurations and then switching on the force on the intruder. The final velocity of the intruder was measured after stationarity was attained. The time to reach a stationary state depends on the packing fraction, the restitution coefficient and the applied force on the intruder. A typical trajectory and its corresponding fluctuating velocity for the largest packing fraction η = 0.8 and force F = 100000 is shown the inset of Fig. 4. The intruder moves approximately 3R 0 in x-direction before fluctuating around its stationary average velocity.

Results
The time rate of change of the total momentum P = 1 N i m i v i follows from Eq. (4): which is solved by Here we have used that the collisions as well as the random driving conserve momentum. In the frictional case with γ = 0, the total momentum goes to a constant, P x = F/(N γ), whereas for the frictionless case with γ = 0, the momentum grows linearly with time, P x = F t/N . In the following we shall mainly discuss the first case, because the frictional model is closer to experiment. However, the frictionless case allows us to study the propagation of momentum in an inelastic fluid [23,24].

Frictionless state ( γ = 0 )
To that end, we consider a system with aspect ratio A = 5 : 1, N = 5000 and fixed walls, which reflect the particles elastically. Momentum is conserved on average except for the pulling force which constantly feeds (small) momenta into the system, which are propagated by collisions away from the source into all of the available area. Below the jamming transition momentum transport is given by ballistic motion of particles as well as by collisions. If the intruder starts at the left hand side of the system and is pulled by the external force, then only particles in the neighbourhood of the intruder will feel the local momentum input for short times. The momentum given to the intruder is distributed among the particles in front (i.e. in the direction of the pulling  Fig. 1 Momentum transport with γ = 0 and hard walls. A vector denoting the particle's velocity is assigned to each particle's position. (T = 0.1, ǫ = 0.9, F = 500, η = 0.8, t = 2.6, 6. 3, 8.5) force) of the intruder. Again, these particles distribute their momenta by collisions as well as by ballistic transport. A front of particles carrying the momentum fed into the system by the intruder is formed (see Fig. 1 top), propagates through the system (see Fig. 1 center) and ultimately collides with the hard wall (see figure 1 bottom). At this instant the momentum is reflected by the wall and many particles end up in a highly compactified state with a packing fraction η ≈ 0.839.
To analyze the momentum wave quantitatively, we plot P x averaged over z in Fig. 2. The propagation front is well defined, its center of mass moves with velocity V cm = 16.51 and it broadens with time such that it's width increases linearly with time, ∆ = V br t (with V br = 10.89). Both velocities increase with pulling force, while their ratio is approximately constant. Hence the observed propagation front cannot be identified with linear sound but is presumably a shockwave in agreement with the observations of Ref. [24].

Force-velocity relation ( γ = 0 )
In the following, we consider a system with N = 20000 grains, subject to friction (γ = 1) with hard walls in the z-direction and periodic boundary conditions in the xdirection, allowing for a stationary current. We analyze the steady state motion of the intruder for small and large forces as a function of packing fraction for T = 1 and subsequently compare these results to the same system with T = 0.04.

Linear Regime: Mobility of the intruder
For a small pulling force we expect a linear relation between the velocity of the intruder and the force: defining the mobility µ of the intruder. In Fig. (3) we plot the velocity of the intruder versus force and indeed do observe a linear regime for small pulling force. From the slope we extract the mobility which is shown in Fig.  (4). The breakdown of the mobility when the glass transition is approached is clearly visible for both investigated ǫ. The more inelastic system with ǫ = 0.7 shows increased mobilities as compared to ǫ = 0.9. Since the particles are more "sticky", they tend to stay closer after a collision than in the elastic limit, i.e., stronger density fluctuations are inherent to systems with stronger inelasticity. This enlarges the accessible space for the intruder compared to the system with ǫ = 0.9, increasing its average velocity. Moreover, the momentum of the intruder in the direction of the pulling force is reduced due to collisions with other particles in front.
In the center of mass frame, the intruder is reflected backwards. For the more inelastic system this kind of backscattering is less effective (see Eq. 1), so that the velocity in the direction of the force and hence the mobility is larger for the more inelastic system.

Nonlinear regime
As the pulling force is increased, deviations from linear behaviour are expected and observed. To investigate these systematically we have applied forces in the range 1 ≤ F ≤ 10 5 and show our data in a scaling plot in Fig.  (5). Scaling velocities by v dr = p dr /m eff and forces by v 2 dr · η 1/2 collapses the data for large forces. The velocity of the intruder scales algebraically with the pulling force, according to v I ∝ F β , β = 0.55.
The crossover between the linear and the algebraic regime is accompanied by range of forces in which the intruder velocity increases superlinearly with the applied force. This is observed at high packing fractions only, see Fig. 3. We conjecture that this shear thinning is due to the formation of vortices which can only be formed at high packing fractions, where momentum is conveyed almost instantaneously as the disks are very close to each other. For low packing fractions, the dis- tance of neighbouring particles is too large for the formation of vortices because the damping γ dissipates most of the momentum.

Temperature dependence
The crossover from linear to nonlinear response depends on the thermal velocity, v th ≃ 2T /m 0 ≈ 0.7. For small packing fractions, this crossover actually occurs at v th as can be seen e.g. in Fig. 3 for the smallest packing fractions. Decreasing the temperature is expected to shrink the linear regime also for higher densities. Hence we try to explore the force velocity relation in the same system but with smaller thermal velocity. Here we choose T = 0.04. For low pulling forces, we expect the mobility to be larger than in the case T = 1 since the decreased thermal motion does not disturb the intruder travelling through the almost resting surrounding disks. For high forces, we expect the intruder velocity to not depend on the temperature, since in this case, the intruder's velocity is at least one order of magnitude higher than the thermal velocity (see Fig. 5).
These expectations are indeed born out by the data. We plot the force velocity relation for packing fractions η = 0.6 and 0.775 for both temperatures in Fig. 6. For high pulling forces F ≥ 10 3 , the intruder velocities for low and high temperature collapse for both packing fractions as expected. For η = 0.6 and T = 0.04, the crossover from the linear regime to the nonlinear regime is shifted to smaller forces. For η = 0.775 and T = 0.04, the linear regime is not even detectable. The crossover is shifted by roughly two decades. To explore the linear regime would require forces as small as F = 10 −2 , at which the average intruder velocity would be too small to be detectable.

Conclusion and Outlook
We have generalised an event-driven algorithm to include friction. As a first application, we have studied the dynamics of an intruder pulled by an external force. In contrast to previous simulations [8], we consider a fluidized granular medium, which is expected to undergo a glass transition at a packing fraction below random close packing [9]. We do indeed find a dramatic decrease of the mobility around η = 0.8, in agreement with previous simulations without pulling force and consistent with a glass transition. For large pulling force, the data can be collapsed by scaling, following a power law dependence v I ∝ F β with β = 0.55.
In the frictionless case the pulling force generates a momentum wave propagating through the sample and thereby compactifying it. We plan to investigate momentum transport in more detail in the future. Furthermore the generalised event driven algorithm will be useful more generally in the context of frictional granular matter fluidized by air or water flow. Work along these lines is in progress.