text:
# Phys. Rev. E, Vol. 53, 5382 (1996) A model for collisions in granular gases

## I Introduction

A rich variety of systems one encounters in nature may be considered as “granular gas” Campbell . The most important difference between a “gas” of granular particles and a regular gas is the inelastic nature of the inter–particle collisions. The steady removal of kinetic energy from the granular gas due to dissipative collisions causes a variety of non–equilibrium processes, which have been subjects of experimental (e.g. exp1 ; exp2 ; exp3 ; exp4 ; exp5 ; exp6 ; exp7 ; exp8 ; exp9 ) and theoretical (e.g. theo1 ; theo2 ; theo3 ; theo4 ; theo5 ) interest. Particularly in recent time many of the experimental results have been reproduced and investigated using various techniques such as Cellular Automata (e.g. CA1 ; CA2 ; CA3 ), Monte–Carlo Methods MonteCarlo , Lattice–Gas models LatticeGas , Molecular Dynamics in two twoDMD1 ; twoDMD2 ; twoDMD3 ; twoDMD4  and three threeDMD1 ; threeDMD2 ; threeDMD3  dimensions and hybrid methods hybrid1 ; hybrid2 ; hybrid3 ; hybrid4 .

The loss of kinetic energy of a pair of inelastically colliding grains can be described using the restitution coefficients for the normal and tangential components of the relative motion $`ϵ^N`$ and $`ϵ^T`$

$`\left(\stackrel{}{g}^N\right)^{}`$ $`=`$ $`ϵ^N\stackrel{}{g}^N\left(0ϵ^N1\right)`$ (1)
$`\left(\stackrel{}{g}^T\right)^{}`$ $`=`$ $`ϵ^T\stackrel{}{g}^T\left(1ϵ^T1\right),`$ (2)
where $`\stackrel{}{g}^N`$ and $`\stackrel{}{g}^T`$ are the relative velocities of the particles in normal and tangential direction before the collision and $`\left(\stackrel{}{g}^N\right)^{}`$, $`\left(\stackrel{}{g}^T\right)^{}`$ after the collision.

Recently, the collision properties of small spheres have been investigated experimentally Foerster . These investigations have shown that the type of the collision (sliding or sticking) depends on the ratio of $`g^N`$ and $`g^T`$. The results were explained with different models for each type, and the coefficients in these models were independent from the velocity. On the other hand, laboratory experiments with ice balls Bridges  as well as with spheres of other materials (for an overview see impact ) have shown that the normal restitution coefficient $`ϵ^N`$ depends significantly on the impact velocity. As already seen, the tangential restitution coefficient depends on the impact parameters as well.

The behaviour of the sheared granular material may be significantly different if the restitution coefficients depend on the impact velocity. This dependence should be taken into account in order to get an adequate model of the stress distribution LunSavage . It is also known that the the parameters $`ϵ^N`$ and $`ϵ^T`$ crucially influence the global dynamics of granular systems (e.g. ArakiTremaine ; WisdomTremaine ).

In the present study we investigate how the restitution coefficients depend on the relative impact velocity. For the normal component of the relative motion we derive an expression for the normal force acting between the colliding particles, which accounts for the dissipation in the bulk of material. One particular application of the results presented here is the explanation of experiments with ice balls Bridges  which are of importance for the investigation of the dynamics of planetary rings Planetenbuch . A static model for the inelastic impact of metal bodies was presented in Johnson  which based on the assumption of fully plastic indentation and constant mean contact pressure and led analytically to a proportionality $`ϵ^N\left(g^N\right)^{1/4}`$ for arbitrary material constants. In contrary, our quasistatic approach does not request other additional assumptions and can be adapted to different experimental results by changing the coefficients in the differential equation which describes the time dependence of the deformation. From these coefficients, material coefficients can be estimated Caputh .

Our result contains the Hertz theory of elastic impact Hertz  and the theory of the viscoelastic impact by Pao pao  as special cases. For the tangential component of the relative motion we consider a microscopic model of the contact of colliding particles. We derive a mean–field expression for the tangential inter–particle force. The result contains the model of the tangential force of colliding particles by Haff and Werner CundallStrack ; HaffWerner  as a special case, and we are able to treat different tangential collisional behaviors within the framework of one single model.

In section II we formulate the collision model and derive the equations for the normal and tangential components of the relative motion of the colliding grains. In section III we present the results for the restitution coefficients for the proposed model and discuss the dependence of the coefficients on the components of the impact velocity. A model for the dynamics of granular materials is proposed. In the last section IV we summarize the results. Details of the derivations are given in the appendixes A and B.

## II The collision model

We consider the inelastic collision between two spherical particles $`i`$ and $`j`$. The values $`\stackrel{}{r}_i`$, $`R_i`$, $`\dot{\stackrel{}{r}}_i`$, $`\dot{\stackrel{}{\omega }}_i`$, $`m_i`$, and $`J_i`$ are the position of the center of sphere $`i`$, its radius, velocity, angular velocity, mass and momentum of inertia, respectively. The relative velocity of the surfaces of the colliding particles at the point of contact is given by (e.g. CundallStrack ; impact )

$`\stackrel{}{g}_{ij}`$ $`=`$ $`\left(\dot{\stackrel{}{r}}_i\stackrel{}{\omega }_i\times R_i\stackrel{}{n}\right)\left(\dot{\stackrel{}{r}}_j+\stackrel{}{\omega }_j\times R_j\stackrel{}{n}\right)`$ (3)
$`=`$ $`\dot{\stackrel{}{r}}_i\dot{\stackrel{}{r}}_jR_i\stackrel{}{\omega }_i\times \stackrel{}{n}R_j\stackrel{}{\omega }_j\times \stackrel{}{n},`$
with $`\stackrel{}{n}=\frac{\stackrel{}{r}_i\stackrel{}{r}_j}{\stackrel{}{r}_i\stackrel{}{r}_j}`$. Introducing the dimensionless moment of inertia $`\widehat{J}_i`$, the effective mass $`m_{ij}^{\text{eff}}`$ and the effective radius $`R_{ij}^{\text{eff}}`$

$`\widehat{J}_i`$ $`=`$ $`{\displaystyle \frac{J_i}{m_iR_i^2}}`$ (4)
$`m_{ij}^{\text{eff}}`$ $`=`$ $`{\displaystyle \frac{m_im_j}{m_i+m_j}}`$ (5)
$`R_{ij}^{\text{eff}}`$ $`=`$ $`{\displaystyle \frac{R_iR_j}{R_i+R_j}}`$ (6)
one obtains Newtons equations for the translational and rotational motion

$$\frac{\text{d}\stackrel{}{g}_{ij}}{\text{d}t}=\frac{\stackrel{}{F}_{ij}}{m_{ij}^{\text{eff}}}+\left(\frac{1}{\widehat{J}_im_i}+\frac{1}{\widehat{J}_jm_j}\right)\left(\stackrel{}{n}\times \stackrel{}{F}_{ij}\right)\times \stackrel{}{n}.$$
(7)
The force $`\stackrel{}{F}_{ij}`$ acting between the particles during collision consists of the normal component $`\stackrel{}{F}_{ij}^N=\stackrel{}{n}\left(\stackrel{}{n}\stackrel{}{F}_{ij}\right)`$ and the tangential component$`\stackrel{}{F}_{ij}^T=\stackrel{}{F}_{ij}\stackrel{}{F}_{ij}^N`$. Introducing the corresponding components $`\stackrel{}{g}_{ij}^N`$ and $`\stackrel{}{g}_{ij}^T`$ of the relative velocity $`\stackrel{}{g}_{ij}`$ and with

$$\kappa _{ij}^1=1+\frac{m_i\widehat{J}_i+m_j\widehat{J}_j}{\widehat{J}_i\widehat{J}_j\left(m_i+m_j\right)}$$
(8)
we rewrite eq. (3) omitting the indexes $`ij`$:

$`\dot{\stackrel{}{g}}^N`$ $`=`$ $`\stackrel{}{F}^N/m^{\text{eff}}`$ (9)
$`\dot{\stackrel{}{g}}^T`$ $`=`$ $`{\displaystyle \frac{1}{m^{\text{eff}}\kappa }}\stackrel{}{F}^T.`$ (10)
Using eqs. (I) the energy loss during the collision is

$$\mathrm{\Delta }Q=\frac{m^{\text{eff}}}{2}\left(\stackrel{}{g}^N\right)^2\left(\left(ϵ^N\right)^21\right)+\frac{m^{\text{eff}}}{2}\kappa \left(\stackrel{}{g}^T\right)^2\left(\left(ϵ^T\right)^21\right).$$
(11)
The energy is conserved during the collision if $`ϵ^N=1`$ and $`ϵ^T=\pm 1`$. In these cases there is a completely elastic rebound for the normal component and either completely elastic rebound (rough spheres) or frictionless slipping (smooth spheres) for the tangential component.

### II.1 Normal motion

We assume that the colliding particles begin to touch each other at the time $`t=0`$ with the relative normal velocity $`\stackrel{}{g}^N`$. When we introduce the deformation (or “compression”)

$$\xi (t)=R_i+R_j\left(\stackrel{}{r}_i(t)\stackrel{}{r}_j(t)\right)$$
(12)
this velocity can be written as $`g^N=\stackrel{}{g}^N=\dot{\xi }`$.

Thus from eq. (9) we obtain the equations

$`\ddot{\xi }(t)`$ $`=`$ $`F^N\left(\xi (t)\right)/m^{\text{eff}}`$ (13)
$`\dot{\xi }(0)`$ $`=`$ $`g^N`$
$`\xi (0)`$ $`=`$ $`0.`$
The normal force $`F^N`$ consists of an elastic, conservative part due to the deformation $`\xi `$ of the particles and a viscous part due to dissipation of energy in the bulk of the particle material, depending on the deformation rate $`\dot{\xi }`$. For the conservative part Hertz’s theory of elastic contact Hertz  gives for spherical particles

$$F_{(el)}^N(\xi )=\frac{2Y}{3\left(1\nu ^2\right)}\sqrt{R^{\text{eff}}}\xi ^{3/2},$$
(14)
where $`Y`$ and $`\nu `$ are the Young modulus and the Poisson ratio for the material the particles consist of. This relation between the elastic component of the force and the deformation is valid for the quasistatic regime of the collision, i.e. when inertial and relaxation effects may be neglected (see appendix B).

The existing phenomenological expressions for the dissipative part of the normal force which are either linear in the deformation rate $`\dot{\xi }`$ (e.g. CundallStrack ; Dilley ) or quadratic Poeschl , however, do not agree satisfactorily with the experimental data for the normal restitution coefficient Bridges . Pao pao  extended the Hertz theory of impact for the viscoelastic case, where, however, the dependence of the bulk dissipation on the dilatation rate was neglected. In this theory memory effects in the dissipative processes were taken into account. Although the latter approach is not self–consistent (see appendix B), it predicts a power–law dependence of the dissipative force on the deformation rate, yielding an exponent similar to that one for the quasistatic collision. In the present study we develop a self–consistent quasistatic approximation to calculate the normal force acting between colliding viscoelastic particles. The quasistatic approximation is valid when the characteristic relative velocity of the granular particles is much less than the speed of sound in the material which is satisfied for many experimental situations even in astrophysical systems such as planetary rings SpeedPlanetaryRings . For the duration of the collision it is necessary to exceed significantly the viscous relaxation time in the material of colliding particles (see appendix B).

Different from the approaches pao ; impact  we take into account both components of the dissipative force, arising from the shear strain rate as well as from the dilatation rate, which are both of comparable importance for the normal component of the relative motion. From the equation of motion for the viscoelastic continuum we find the general relation between the dissipative part of the normal force and the deformation rate. We show that memory effects in dissipative processes may be neglected in the case of a self–consistent quasistatic approximation. Since the calculation of the dissipative part of the normal force is rather straightforward, we present only the main idea of the derivation and refer to appendix A for further details. In appendix B the conditions for the validity of the quasistatic approach are given.

The total normal force acting between viscoelastic particles may be derived from a stress tensor combined of an elastic and a dissipative part Landau

$$\widehat{\sigma }=\widehat{\sigma }_{(el)}+\widehat{\sigma }_{(dis)}$$
(15)
with

$`\widehat{\sigma }_{(el)}`$ $`=`$ $`E_1\left[{\displaystyle \frac{1}{2}}\left\{\stackrel{}{u}+\stackrel{}{u}\right\}{\displaystyle \frac{1}{3}}\widehat{I}\stackrel{}{u}\right]+`$ (16)
$`E_2\widehat{I}\stackrel{}{u}`$
$`\widehat{\sigma }_{(dis)}`$ $`=`$ $`\eta _1\left[{\displaystyle \frac{1}{2}}\left\{\dot{\stackrel{}{u}}+\dot{\stackrel{}{u}}\right\}{\displaystyle \frac{1}{3}}\widehat{I}\dot{\stackrel{}{u}}\right]+`$ (17)
$`\eta _2\widehat{I}\dot{\stackrel{}{u}}.`$
The displacements in the material is denoted by $`\stackrel{}{u}`$ and $`\widehat{I}`$ is the unit tensor. $`E_{1/2}`$ and $`\eta _{1/2}`$ are the elastic and the viscous constants of the particle material

$`E_1`$ $`=`$ $`{\displaystyle \frac{Y}{1+\nu }}`$ (18)
$`E_2`$ $`=`$ $`{\displaystyle \frac{Y}{3\left(12\nu \right)}}.`$ (19)
In the quasistatic regime (appendix B) the displacement field $`\stackrel{}{u}(\stackrel{}{r},t)`$ can be approximated by that of the static problem $`\stackrel{}{u}\left(\stackrel{}{r}\right)`$. It is completely determined by the elastic component of the inter–particle force (14). Thus, the displacement velocities can be written as

$$\dot{\stackrel{}{u}}(\stackrel{}{r},t)\dot{\xi }\frac{}{\xi }\stackrel{}{u}_{(el)}(\stackrel{}{r},\xi ),$$
(20)
where $`\stackrel{}{u}_{(el)}(\stackrel{}{r},\xi )`$ is the solution of the static (elastic) contact problem. This expression depends parametrically on the deformation $`\xi `$ and the dissipative part of the stress tensor becomes

$`\widehat{\sigma }_{(dis)}`$ $`=`$ $`\dot{\xi }{\displaystyle \frac{}{\xi }}\{\eta _1[{\displaystyle \frac{1}{2}}(\stackrel{}{u}_{(el)}+\stackrel{}{u}_{(el)})`$ (21)
$`{\displaystyle \frac{1}{3}}\widehat{I}\stackrel{}{u}_{(el)}]+\eta _2\widehat{I}\stackrel{}{u}_{(el)}\}.`$
The calculations can be significantly simplified when we notice that the elastic and the dissipative parts of the stress tensor are related in the quasistatic regime (see eqs. (II.1) and (20))

$$\sigma _{(dis)}=\dot{\xi }\frac{}{\xi }\sigma _{(el)},(E_1\eta _1,E_2\eta _2).$$
(22)
Therefore the impact problem for the viscoelastic particles in the quasistatic regime can be mapped onto the corresponding problem for elastic particles. Performing calculations similar to that of the elastic case (for details see Hertzsch  and appendix A) one can find an expression for the dissipative part of the normal force:

$`F_{(dis)}^N`$ $`=`$ $`{\displaystyle \frac{Y}{\left(1\nu ^2\right)}}\sqrt{R^{\text{eff}}}A\sqrt{\xi }\dot{\xi }`$ (23)
$`A`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \frac{\left(3\eta _2\eta _1\right)^2}{\left(3\eta _2+2\eta _1\right)}}\left[{\displaystyle \frac{\left(1\nu ^2\right)\left(12\nu \right)}{Y\nu ^2}}\right].`$
From eqs. (23) and (14) we obtain for the normal component of the relative motion

$$\ddot{\xi }+\frac{2Y\sqrt{R^{\text{eff}}}}{3m^{\text{eff}}\left(1\nu ^2\right)}\left(\xi ^{3/2}+\frac{3}{2}A\sqrt{\xi }\dot{\xi }\right)=0$$
(24)
with the initial conditions $`\dot{\xi }\left(0\right)=g^N`$, $`\xi \left(0\right)=0`$. In the case of $`A\dot{\xi }\xi `$ equation (24) results from a Taylor–expansion of

$$\ddot{\xi }+\frac{2Y\sqrt{R^{\text{eff}}}}{3m^{\text{eff}}\left(1\nu ^2\right)}\left(\xi +A\dot{\xi }\right)^{3/2}=0$$
(25)
which formally coincides with the corresponding equation for the elastic problem, provided that $`\xi `$ is substituted by $`\xi +A\dot{\xi }`$.

It has to be noted, that $`\xi `$ has its minimum at the beginning of the collision where $`\dot{\xi }`$ takes its maximum. Hence, the condition $`A\dot{\xi }\xi `$ is not provided at the very beginning of the contact. On the other hand, the good confirmation of experimental facts (Bridges ) by the numerical solution of equation (25) points to its suitability for at least the rest of the collision time span.

Taking into account $`\left(\dot{g}^N\right)^{}=\dot{\xi }(t_c)`$ ($`t_c`$ – the duration of the collision), the normal restitution coefficient is obtained from

$$ϵ^N=\dot{\xi }(t_c)/\dot{\xi }(0).$$
(26)
### II.2 Tangential motion

In the idealized model the surface of contact between the spheres $`S`$ is a perfectly flat circular area with radius $`R_S=\sqrt{2R^{\text{eff}}\xi (t)}`$. For the description of the tangential forces between the surfaces we follow a current model of tribology (e.g. tribology1 ; tribology2 ) where the apparent surface of contact is built up of a large number of hierarchically ordered asperities varying in shape and size by several decades. For the processes of the momentum transmission we will take into account only the largest–scale asperities (“primary asperities”). The surface asperities do not affect the normal motion, if they are small enough (see appendix A), however they are responsible for the tangential forces, acting between the colliders. Here we consider a simplified mean-field approach and introduce the normal $`\overline{\sigma ^N}`$ and tangential stress $`\overline{\sigma ^T}`$ averaged over the contact area. Further we define the normal component of the total contact area of the asperities of both spheres $`S^N`$ which is responsible for the transmission of the normal force. Correspondingly the tangential projection of the area $`S^T`$ is responsible for the transmission of the tangential force. These surfaces are related to the apparent contact area by the relations Greenwood

$`S^N(t)`$ $`=`$ $`f^N\left(\overline{\sigma ^N}\right)S(t)`$ (27)
$`S^T(t)`$ $`=`$ $`f^T\left(\overline{\sigma ^N}\right)S(t),`$ (28)
where the coefficients $`f^N`$ and $`f^T`$ depend on the average normal stress $`\overline{\sigma ^N}`$. When the spheres begin to touch each other, i.e. $`S=0`$ and $`\overline{\sigma ^N}=0`$, we find $`f^N(0)=0`$ and $`f^T(0)=0`$. We expand the coefficients in eq. (28) with respect to $`\overline{\sigma ^N}=0`$. The linear expansion yields for the tangential component of the surface

$`S^T(t)`$ $`=`$ $`\varphi ^T\overline{\sigma ^N}S(t)`$ (29)
$`\varphi ^T`$ $`=`$ $`\left[{\displaystyle \frac{f^T}{\overline{\sigma ^N}}}\right]_{\overline{\sigma ^N}=0}.`$
For a given model of the sizes and shapes of the asperities one can calculate the value of $`\varphi ^T`$ Greenwood . In the case that the heights of the asperities obey a Gaussian probability distribution with mean value $`L`$ one finds

$$\varphi ^T\sqrt{L}.$$
(30)
For the average size of the asperities $`L`$ of the surfaces the mean field approach yields the average shear deformation $`\overline{\eta }=b\frac{\zeta }{L}`$. The values $`\zeta `$ and $`b`$ are the relative tangential shift of the particle-surfaces and a form-factor, respectively. We assume that the stress is uniformly distributed over the entire surface and find

$$\overline{\sigma ^T}=\frac{Y}{1+\nu }\overline{\eta }=\frac{Y\zeta }{\left(1+\nu \right)L}b.$$
(31)
The linear relation between $`\overline{\sigma ^T}`$ and $`\overline{\eta }`$ holds only for the elastic regime, i.e. only if $`\overline{\sigma ^T}`$ does not exceed some critical value $`\sigma _{}^T`$ which is a specific material constant. If the shear stress exceeds this threshold $`\sigma _{}^T`$, the asperity which hinders the tangential relative motion of the surfaces, is assumed to break resulting in a sudden release of the shear stress. At the same time the surfaces are shifted macroscopically with respect to each other by

$$\zeta _0=\frac{L\left(1+\nu \right)}{bY}\sigma _{}^T,$$
(32)
and one finds

$`\overline{\eta }\left(\zeta \right)`$ $`=`$ $`\eta _{}\left({\displaystyle \frac{\zeta }{\zeta _0}}{\displaystyle \frac{\zeta }{\zeta _0}}\right)`$ (33)
$`\eta _{}`$ $`=`$ $`{\displaystyle \frac{b\zeta _0}{L}}`$
where $`x`$ denotes the integer of $`x`$. The breaking of the asperities dissipates the energy which was previously stored in the elastic stress, i.e. fracturing of the asperities is the elementary dissipative process in the tangential motion. From eq. (33) we obtain the shear stress as a function of the tangential displacement

$$\overline{\sigma ^T}\left(\zeta \right)=\sigma _{}^T\left(\frac{\zeta }{\zeta _0}\frac{\zeta }{\zeta _0}\right)$$
(34)
and the tangential component of the inter–particle force

$`F^T`$ $`=`$ $`S^T\overline{\sigma ^T}\left(\zeta \right)=\varphi ^T\overline{\sigma ^N}S\sigma _{}^T\left({\displaystyle \frac{\zeta }{\zeta _0}}{\displaystyle \frac{\zeta }{\zeta _0}}\right)`$ (35)
$`=`$ $`\mu F^N\left({\displaystyle \frac{\zeta }{\zeta _0}}{\displaystyle \frac{\zeta }{\zeta _0}}\right).`$
where $`F^N=\overline{\sigma ^N}S`$ is the normal component of the inter–particle force and $`\mu =\varphi ^T\sigma _{}^T`$.

It may be shown, that a more refined mean-field approach, which does not use the assumption of the uniformly distributed stress over the contact area, leads to the same eq. (35) for the tangential motion.

From eq. (35) follows the condition for the maximum tangential force

$$F_{\text{max}}^T=\mu F^N.$$
(36)
Thus our model reproduces the Coulomb friction law Coulomb  with the friction coefficient $`\mu `$, expressed in terms of microscopic parameters.The model for the tangential motion is very similar to the extensively investigated one dimensional model by Burridge and Knopoff BurridgeKnopoff ; CalsonLanger  intended to model earthquakes.

With $`g^T(t)=\dot{\zeta }(t)`$, eq. (8) and $`F^N=m^{\text{eff}}\ddot{\xi }(t)`$ the tangential motion is governed by the differential equation

$$\ddot{\zeta }\frac{\mu }{\kappa }\ddot{\xi }\left(t\right)\left(\frac{\zeta }{\zeta _0}\frac{\zeta }{\zeta _0}\right)=0$$
(37)
with the initial conditions $`\dot{\zeta }\left(0\right)=g^T`$ and $`\zeta \left(0\right)=0`$. The value of $`\ddot{\xi }(t)`$ is given by eqs. (24) and (25). Then the tangential restitution coefficient reads

$$ϵ^T=\dot{\zeta }\left(t_c\right)/\dot{\zeta }(0).$$
(38)
## III Results and discussion

The obtained equations for the normal (eq. (25), (24)) and tangential motion (eq. (37)) have been solved numerically using a Runge–Kutta Method of fourth order with adaptive step size recipes . The restitution coefficients $`ϵ^T`$ and $`ϵ^N`$ have been calculated as functions of the normal and tangential relative velocities $`g^T`$, $`g^N`$. For the integration we used the parameters of ice at low temperatures LandoltBoernstein : Young modulus $`Y=10GPa`$, the Poisson ratio $`\nu =0.3`$, the particle size $`R=10^2m`$ with density $`\rho =10^3kgm^3`$. The coefficient $`A`$ in eq. (24) was considered to be a fit parameter, due to lack of information about the dissipative coefficients $`\eta _1`$ and $`\eta _2`$. Fig. 1 shows the numerical result of our model for the normal restitution coefficient $`ϵ^N`$ as a function of the normal relative velocity $`g^N`$ compared to experimental data for the collision of spherical ice particles with an ice wall Bridges . The experimental results are well reproduced by our model.

For the investigation of the tangential restitution coefficient of colliding homogeneous spheres ($`J=\frac{2}{5}mR^2`$, $`\kappa =\frac{2}{7}`$) we have chosen the Coulomb friction coefficient from the interval $`\mu [10^2\mathrm{}1]`$. The value of $`\sigma _{}^T`$ is a material constant. With eq. (30) and the definitions of $`\xi _0`$ and $`\eta _{}`$ we estimate $`\zeta _0`$ which characterizes the size of the surface asperities via $`\mu =\alpha \sqrt{\zeta _0}`$. In the numerical calculation we have chosen $`\alpha =1`$. The results are shown in fig. 2. The tangential restitution coefficient $`ϵ^T`$ is drawn versus the plane defined by the tangential and normal velocities $`g^T`$ and $`g^N`$. The three plots correspond to the values of the asperity sizes $`\zeta _0=(10^7;210^4;10^3)R^{\text{eff}}`$, respectively.

The obvious common feature of all cases is sliding of the surfaces ($`ϵ^T>0`$) for small $`g^N`$ and large $`g^T`$. This is quite plausible since smaller impact velocity $`g^N`$ corresponds to a smaller normal acceleration and thus, to a smaller value of the maximal tangential force, eq. (35). As a result $`ϵ^T1`$ at $`g^N0`$ due to vanishing tangential acceleration. At the same time, for the high tangential velocity, $`g^T(0)1`$, ($`g^N(0)1`$) sliding occurs owing to a considerable breaking of the asperities. The area of the sliding phase in the $`g^N`$$`g^T`$ plane depends on the size $`\zeta _0`$ of the asperities.

In the case of $`\zeta _0=10^7R^{\text{eff}}`$ sliding occurs in the entire velocity range according to values $`0.85ϵ^T1`$. The small asperities are not able to cause a sufficient torque to change considerably the spin of the individual particles. Here we are close to the case of ideal smooth spheres where no change of the tangential motion is expected ($`ϵ^T=1`$).

In the other two cases $`\zeta _0=(210^4;10^3)R^{\text{eff}}`$ one recognizes two phases:

1. Sliding $`ϵ^T>0`$ at small $`g^N`$ and high $`g^T`$
2. Reversal of the spin of either particles $`ϵ^T<0`$ at small $`g^T`$ and higher $`g^N`$.

Case 1 corresponds to the effect discussed in the context of $`\zeta _0=10^7`$. Despite of being far from rather smooth spheres, the small tangential force originated from small $`g^N`$, changes the velocity $`g^T`$ only slightly. Hence one has $`ϵ^T>0`$, which is also the case for high velocities $`g^T`$ where the asperities break. In the case 2 we have the other extreme: a high normal acceleration causes a tangential force which is high enough to change the sign of $`g^T`$ as long as the asperities do not break (small $`g^T`$). A complete reversal of the tangential velocity according to $`g^T(0)g^T\left(t_c\right)`$ is not possible because of the dissipation arising of the bulk viscosity of the material which enters the normal as well as the tangential forces (see eqs. (24),(37)).

Both types of behavior of the tangential motion are separated by a sharp transition at $`g^T=g_{(cr)}^T`$ where the asperities begin to break (see surface plots for $`\zeta _0=(210^4;10^3)R^{\text{eff}}`$). The higher $`g^N`$ the larger the critical tangential velocity $`g_{(cr)}^T`$. A higher normal velocity $`g^N`$ causes a stronger counteracting force $`F^T`$ and thus a larger tangential impact speed $`g^T`$ is necessary to reach the critical deformation where the asperities break. Both cases ($`\zeta _0=(210^4;10^3)R^{\text{eff}}`$) reveal similar qualitative behavior but the ranges of different type of motion (12) cover different areas in the $`g^T`$$`g^N`$ plane.

The results show that our model includes a continuous transition from the limit case of rough spheres ($`ϵ^T1`$) to the limit case of smooth spheres ($`ϵ^T1`$). In the literature of the dynamics of granular material an alternative step function is widely used for the tangential force CundallStrack

$$F^T=\mathrm{min}\{\gamma _sm^{\text{eff}}\left|g^T\right|,\mu F^N\}.$$
(39)
The numerical evaluation of the considered model (fig. 2) reveals surprising behavior of the tangential restitution coefficient $`ϵ^T`$ as a function of the normal velocity $`g^N`$ at fixed tangential velocity $`g^T`$. (This effect is noticeable for the largest values of $`\zeta _0`$.) At low and moderate $`g^N`$, $`ϵ^T`$ first decreases with increasing $`g^N`$ down to its minimal negative value in a manner discussed above, but when $`g^N`$ exceeds some threshold (approximately of several $`g^T`$), it starts to increase up to zero at very high values of $`g^N`$. This effect may be explained as follows: For high values of $`g^N`$ the average normal force is large and causes thus a large tangential force, which can effectively decelerate the initial tangential velocity without switching to the sliding regime.

Calculating the restitution coefficients $`ϵ^T`$, $`ϵ^N`$, (in the limits of our model) we obtain a complete description of binary collisions. Therefore one can determine the dynamics for moderately dense granular gases, where an evolution occurs via a sequence of binary collisions. For such systems we have the following Boltzmann equation for the one–particle distribution function

$$\left(\frac{}{t}+\stackrel{}{v}_1\right)f(1)=𝑑\stackrel{}{v}_2𝑑\stackrel{}{\omega }_2𝑑\stackrel{}{n}\left|\stackrel{}{g}^N\stackrel{}{n}\right|\mathrm{\Theta }\left(\stackrel{}{g}^N\stackrel{}{n}\right)\left[\frac{f\left(1^{}\right)f\left(2^{}\right)}{\left(ϵ^Nϵ^T\right)^2}f\left(1\right)f\left(2\right)\right]$$
(40)
with $`\mathrm{\Theta }(x)`$ given by

$$\mathrm{\Theta }(x)=\{\begin{array}{ccc}\text{1}\hfill & \text{ for}\hfill & x0\hfill \\ \text{0}\hfill & \text{ for}\hfill & x<0\text{ .}\hfill \end{array}$$
(41)
and with the common notations, e.g. $`(1)=\{\stackrel{}{r}_1,\stackrel{}{v}_1,\stackrel{}{\omega }_1,t\}`$. The velocity and angular velocity of the first particle after the collision $`\stackrel{}{v}_1^{}`$ and $`\stackrel{}{\omega }_1^{}`$ can be expressed in terms of the pre–collisional values via the relations

$`\stackrel{}{v}_1^{}`$ $`=`$ $`\stackrel{}{v}_1+{\displaystyle \frac{m^{\text{eff}}}{2m_1}}\left\{\left[ϵ^T(g^N,g^T)1\right]\stackrel{}{g}^T\left[ϵ^N(g^N,g^T)+1\right]\stackrel{}{g}^N\right\}`$ (42)
$`\stackrel{}{\omega }_1^{}`$ $`=`$ $`\stackrel{}{\omega }_1+{\displaystyle \frac{m^{\text{eff}}}{2m_1}}R\stackrel{}{n}\times \left\{\left[ϵ^T(g^N,g^T)1\right]\stackrel{}{g}^T\left[ϵ^N(g^N,g^T)+1\right]\stackrel{}{g}^N\right\}`$ (43)
and analogously for $`\stackrel{}{v}_2^{}`$, $`\stackrel{}{\omega }_2^{}`$. With the use of eqs. (40) and (III) and the above calculated restitution coefficients $`ϵ^N(g^N,g^T)`$ and $`ϵ^T(g^N,g^T)`$ (eqs. (26) and (38)) one can describe the evolution of moderate dense granular gases without computing the detailed dynamics of binary collisions as it is usually done in the “soft sphere” Molecular Dynamics (MD) approach. Here one considers the grains as elastic bodies which deform each other during a collision. There are several ansatzes for the force acting between touching grains CundallStrack ; HaffWerner ; Walton . In all cases one has to chose a time step for the integration scheme which is significantly smaller than the typical collision time. Hence, during each collision one has to calculate about 10…1000 times the interaction force between the grains to provide satisfying accuracy of the simulation. When two grains approach each other they do not feel any interaction as long as they do not touch each other. When granular particles collide they interact via huge restoring forces which can be expressed by Young moduli of the order of $`Y=10^7kg/msec^2`$. The difficulty of the simulation consists in the extreme short range interaction of the particles and the resulting huge gradient of the interaction force. Therefore presently one cannot simulate much more than 3000 granular particles in three dimensions (e.g. WaltonBraun ; GallasEtAl ) and about $`10^4`$ particles in two dimensions (e.g. PoeschelBuchholtz ).

Another method for the simulation of granular assemblies is the “hard sphere” approach where one does not consider the details of the collision but only the pre– and post–collisional velocities of each pair of colliding grains. The goal of these simulations is the low numerical complexity. One needs only computational effort when particles collide but not in between the collisions. This allows for the application of so called event–driven calculations (e.g. Rapaport ). Hence, one can simulate much more particles than with “soft particle” methods.

One of the preconditions for the application of the “hard sphere approach” is the exact knowledge of the normal and tangential restitution coefficients, $`ϵ^N`$ and $`ϵ^T`$, as functions of the normal and tangential impact rates, $`\stackrel{}{g}^N`$ and $`\stackrel{}{g}^T`$, which theoretically to determine was the goal of the present paper.

An interesting possible application of this approach is the dynamics of planetary rings composed of icy and silicate material which is determined by inelastic dissipative collisions Planetenbuch . The calculation of such systems using the traditional MD is impossible due to the huge number of particles in these systems.

## IV Conclusion

A novel model for collision of particles in granular gases is proposed. For the normal component of the relative motion the equation of motion is derived based on the general consideration of the viscoelastic impact. We find the expression for the dissipative part of the normal force in the self–consistent quasistatic approximation which generalizes the existing results for the viscoelastic collisions pao . For the tangential relative motion we investigated a microscopic model of surfaces of the colliding particles which are in contact. We found a mean–field expression for the tangential inter–particle force, which can reproduce smooth, reflecting or sticky collisions depending on the microscopic parameters of the surfaces and on the relative impact velocity. A frequently used model for collisions of granular particles by Haff and Werner CundallStrack  is contained in our model as a special case. The restitution coefficients for the normal and tangential motion are calculated as functions of the relative impact velocity. A rather nontrivial strongly non–linear dependence of the tangential restitution coefficient on the impact velocity is observed.

The obtained restitution coefficients may be used to describe the dynamics of moderately dense granular gases, where the evolution occurs via a sequence of successive binary collisions.

###### Acknowledgements.

We thank H. J. Herrmann and J. T. Jenkins for discussion and S. Luding for critical remarks. One of us (N.V.B.) is very grateful to Prof. J. Kurths for the warm hospitality in Potsdam.

## Appendix A Generalization of the Hertz theory

We briefly sketch the Hertz theory of elastic impact and give a generalization of this theory for the case of viscoelastic collisions.

In the quasistatic approximation which is used in Hertz’s impact theory it is assumed that the (time dependent) strain and the (time dependent) stress are related in the same manner as in the static case. It may be shown (see appendix B) that this approximation is valid for the elastic case when the characteristic velocity is much less than the speed of sound in the material of the colliding particles. Moreover for the viscoelastic case it is required that the viscous relaxation time of the material is much shorter than the duration of the collision. In the static case the equation of equilibrium reads Landau

$$\widehat{\sigma }_{(el)}=0,$$
(44)
where the elastic stress tensor $`\widehat{\sigma }_{(el)}`$ is expressed in terms of displacements $`\stackrel{}{u}(r)`$ via eq. (II.1). Hence the static eq. (44) can be written as

$`_{}^2\stackrel{}{u}+b^2_{}^2\stackrel{}{u}=0`$ (45)
$`b^2={\displaystyle \frac{4E_1+6E_2}{3E_1}}={\displaystyle \frac{2\left(1\nu \right)}{\left(12\nu \right)}}`$
with the “longitudinal” and “transversal” parts of the Laplacian.

$`_{}^2`$ $`=`$ $``$ (46)
$`_{}^2`$ $`=`$ $`^2_{}^2.`$ (47)
The boundary conditions for the displacements in eq. (45) are formulated on the surface of contact. From geometric considerations it follows that the contact area between two colliding particles is a plane. Using the appropriate coordinate system centered in the middle of the contact region (where we set $`z=0`$) one can write

$$C_1x^2+C_2y^2+u_{z1}+u_{z2}=\xi .$$
(48)
The values $`u_{z1}=u_{z1}(x,y)`$ and $`u_{z2}=u_{z2}(x,y)`$ are the $`z`$–components of the displacements in the materials of the bodies at the plane $`z=0`$, $`\xi `$ is the total deformation (the sum of the deformations of both bodies at the center of the contact area, i.e. at $`x=y=0`$). The constants $`C_1`$ and $`C_2`$ are expressed in terms of radii of curvature of the surfaces in contact (see e.g. Hertz ; Landau ). The values of $`u_{z1}`$ and $`u_{z2}`$ may be expressed in terms of the normal pressure $`P_z(x,y)`$ that acts between the bodies at $`z=0`$ Landau :

$`u_{z1}(x,y,0)`$ $`=`$ $`u_{z2}={\displaystyle \frac{\mathrm{\Lambda }}{\pi }}{\displaystyle \frac{P_z(x^{},y^{})}{r}𝑑x^{}𝑑y^{}}`$ (49)
$`r`$ $`=`$ $`\sqrt{\left(xx^{}\right)^2+\left(yy^{}\right)^2}`$
$`\mathrm{\Lambda }`$ $`=`$ $`{\displaystyle \frac{2E_1+3E_2}{E_1\left(E_1+6E_2\right)}}={\displaystyle \frac{1\nu ^2}{Y}}.`$
For simplicity we assume that the colliding particles are of the same material. The normal pressure $`P_z`$ is related to the total normal force $`F_{(el)}`$

$$P_z(x,y)=\frac{3F_{(el)}}{2\pi ab}\sqrt{1\frac{x^2}{a^2}\frac{y^2}{b^2}},$$
(50)
where $`a`$ and $`b`$ are the semi–axes of the contact ellipse. The latter values as well as the compression $`\xi `$ may be found from the set of equations

$`\xi `$ $`=`$ $`{\displaystyle \frac{F_{(el)}}{\pi }}{\displaystyle \frac{3}{2}}\mathrm{\Lambda }{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dq}{\sqrt{\left(a^2+q\right)\left(b^2+q\right)q}}}`$ (51)
$`C_1`$ $`=`$ $`{\displaystyle \frac{F_{(el)}}{\pi }}{\displaystyle \frac{3}{2}}\mathrm{\Lambda }{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dq}{\left(a^2+q\right)\sqrt{\left(a^2+q\right)\left(b^2+q\right)q}}}`$ (52)
$`C_2`$ $`=`$ $`{\displaystyle \frac{F_{(el)}}{\pi }}{\displaystyle \frac{3}{2}}\mathrm{\Lambda }{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dq}{\left(b^2+q\right)\sqrt{\left(a^2+q\right)\left(b^2+q\right)q}}}.`$ (53)
From the above expressions follows that for all bodies in contact having smooth surfaces (in the mathematical sense) the total force and the deformation are related via the power law

$$F_{(el)}\left(\xi \right)=\stackrel{~}{c}\xi ^{3/2}.$$
(54)
The constant $`\stackrel{~}{c}`$ depends on the elastic properties of the materials and on the local curvatures of the colliding bodies. For the case of the spherical particles one has the Hertz’s law

$$F_{(el)}\left(\xi \right)=\frac{2Y}{3\left(1\nu ^2\right)}\sqrt{R^{eff}}\xi ^{3/2}.$$
(55)
Using this relation between force and deformation and the equation of motion (eq. (13)) one can describe the elastic collision completely. The duration of the collision is (Hertz ; Landau )

$`t_c`$ $`=`$ $`2.94\left({\displaystyle \frac{m^{\text{eff}}}{k}}\right)^{2/5}\left(\left(g^N\right)\right)^{1/5}`$ (56)
$`k^2`$ $`=`$ $`\left({\displaystyle \frac{4}{5}}{\displaystyle \frac{2}{3\mathrm{\Lambda }}}\right)^2R^{\text{eff}}.`$
In the solution of the elastic contact problem the displacement fields $`\stackrel{}{u}_1\left(\stackrel{}{r}\right)`$ and $`\stackrel{}{u}_2\left(\stackrel{}{r}\right)`$ are completely defined by the value of $`F_{(el)}`$ and thus by the value of the deformation $`\xi `$. Hence we write $`\stackrel{}{u}\left(\stackrel{}{r}\right)=\stackrel{}{u}(\stackrel{}{r},\xi )`$, i.e. the displacement field depends explicitly on the compression. Therefore we obtain for the velocity of the displacement in the quasistatic approximation

$$\dot{\stackrel{}{u}}\left(\stackrel{}{r}\right)=\dot{\xi }\frac{}{\xi }\stackrel{}{u}(\stackrel{}{r},\xi )$$
(57)
and correspondingly for the dissipative part of the stress tensor

$`\sigma _{(dis)}^{ik}`$ $`=`$ $`\dot{\xi }{\displaystyle \frac{}{\xi }}\left\{\eta _1u_{ik}+\left(\eta _2{\displaystyle \frac{\eta _1}{3}}\right)u_{ll}\delta _{ik}\right\}`$ (58)
$`=`$ $`\dot{\xi }{\displaystyle \frac{}{\xi }}\sigma _{(el)}^{ik},(E_1\eta _1,E_2\eta _2).`$
We emphasize that the expression in the curled brackets in the above equation coincides with the elastic stress, provided the viscous constants are substituted by the elastic ones. The latter expression for the dissipative stress tensor is written for the case when the memory effects in the viscous processes may be neglected. A more general case is discussed in appendix B.

The $`\sigma _{(el)}^{zz}`$ component of the elastic stress is equal to the normal pressure $`P_z`$ at the plane $`z=0`$,

$`\sigma _{(el)}^{zz}(x,y,0)=`$
$`=E_1{\displaystyle \frac{u_z}{z}}+\left(E_2{\displaystyle \frac{E_1}{3}}\right)\left({\displaystyle \frac{u_x}{x}}+{\displaystyle \frac{u_y}{y}}+{\displaystyle \frac{u_z}{z}}\right)`$
$`={\displaystyle \frac{3F_{(el)}}{2\pi ab}}\sqrt{1{\displaystyle \frac{x^2}{a^2}}{\displaystyle \frac{y^2}{b^2}}}.`$ (59)
With the transformation of the coordinate axes

$`x`$ $`=`$ $`\alpha x^{}`$ (60)
$`y`$ $`=`$ $`\alpha y^{}`$ (61)
$`z`$ $`=`$ $`z^{}`$ (62)
and

$`\alpha `$ $`=`$ $`\left({\displaystyle \frac{\eta _2\frac{1}{3}\eta _1}{\eta _2+\frac{2}{3}\eta _1}}\right)\left({\displaystyle \frac{E_2+\frac{2}{3}E_1}{E_2\frac{1}{3}E_1}}\right)`$ (63)
$`\beta `$ $`=`$ $`{\displaystyle \frac{\left(\eta _2\frac{1}{3}\eta _1\right)}{\alpha \left(E_2\frac{1}{3}E_1\right)}}`$ (64)
$`a`$ $`=`$ $`\alpha a^{}`$ (65)
$`b`$ $`=`$ $`\alpha b^{}`$ (66)
we obtain

$`\eta _1{\displaystyle \frac{u_z}{z}}+(\eta _2{\displaystyle \frac{\eta _1}{3}})\left({\displaystyle \frac{u_x}{x}}+{\displaystyle \frac{u_y}{y}}+{\displaystyle \frac{u_z}{z}}\right)`$ (67)
$`=`$ $`\beta \left(E_1{\displaystyle \frac{u_z}{z^{}}}+(E_2{\displaystyle \frac{E_1}{3}})({\displaystyle \frac{u_x}{x^{}}}+{\displaystyle \frac{u_y}{y^{}}}+{\displaystyle \frac{u_z}{z^{}}})\right)`$
$`=`$ $`\beta {\displaystyle \frac{3F_{(el)}}{2\pi a^{}b^{}}}\sqrt{1{\displaystyle \frac{x^{}^2}{a^{}^2}}{\displaystyle \frac{y^{}^2}{b^{}^2}}}`$
$`=`$ $`\beta \alpha ^2{\displaystyle \frac{3F_{(el)}}{2\pi ab}}\sqrt{1{\displaystyle \frac{x^2}{a^2}}{\displaystyle \frac{y^2}{b^2}}}.`$
Applying the operator $`\dot{\xi }\frac{}{\xi }`$ to the previous expression we obtain the result for the viscous stress. Integrating the viscous stress over the contact area we finally find for the dissipative component of the inter–particles force

$`F_{(dis)}=A\dot{\xi }{\displaystyle \frac{}{\xi }}F_{(el)}\left(\xi \right)`$ (68)
$`A=\alpha ^2\beta ={\displaystyle \frac{1}{3}}{\displaystyle \frac{\left(3\eta _2\eta _1\right)^2}{\left(3\eta _2+2\eta _1\right)}}{\displaystyle \frac{\left(1\nu ^2\right)\left(12\nu \right)}{Y\nu ^2}}.`$ (69)
Thus one obtains for the normal force which acts between the viscoelastic bodies in the quasistatic regime of collision

$$F=const\left(\xi ^{3/2}+\frac{3}{2}A\sqrt{\xi }\dot{\xi }\right).$$
(70)
The constant in eq. (70) coincides with that for the elastic force. For colliding spherical particles we arrive at eq. (23).

The impact theory, sketched above, was developed for the bodies with the smooth surfaces. If the surface asperities are taken into account, one can consider the actual surface as a smooth one (obtained by averaging over the asperities heights), with a small perturbation superimposed due to the presence of the asperities. One can also consider the actual normal displacements and normal pressure as a sum of the averaged (over the asperities) values and the small perturbation. Then it is easy to show that the equations, obtained for the averaged values coincide (due to linearity of the problem) with the corresponding equations for the elastic collision of the smooth bodies. As a result, the relation between the force and deformation, $`\xi `$ is the same as in the Hertz’s theory, provided that $`\xi `$ is defined with a use of the averaged over the asperities radii of the colliders.

Considering the normal motion for the dissipative collisions, one need not to consider the plastic deformation of the asperities, since the size of the asperities is assumed to be very small compared with the radii of the spheres. For our calculations in fig. 2. the asperity size is $`10^3`$ to $`10^7`$ times smaller than the effective Radius of the particles. Hence the dissipation in the bulk of the asperities is negligibly small compared to the total dissipation in the compressed part of the collider. Moreover, the ratio of the normal to tangential stress, may be roughly estimated as, $`\overline{\sigma ^N}/\overline{\sigma ^T}\left(\xi /R\right)^{1/2}`$, so that the crushing of the asperities does not seem to be important for the normal motion, if $`\xi /R1`$ and if the conditions of the quasistatic collision hold. Thus one concludes, that the surface asperities may be ignored, when the normal motion is studied, provided they are small and the conditions of the quasistatic collision are satisfied.

## Appendix B Validity of the quasistatic approximation

To analyze more rigorously the conditions when the quasistatic approximation is valid we write the equation of motion for the viscoelastic continuum

$$\rho \ddot{\stackrel{}{u}}=\left(\widehat{\sigma }_{(el)}+\widehat{\sigma }_{(dis)}\right)$$
(71)
where $`\rho `$ is the density of the material. The expression for the elastic part of the stress tensor is given by eq. (II.1). Taking into account the memory effects of the dissipative processes in the material one can write for the dissipative part

$$\widehat{\sigma }_{(dis)}\left(t\right)=E_1_0^t\psi _1\left(t\tau \right)\left[\frac{1}{2}\left\{\dot{\stackrel{}{u}}(\tau )+\dot{\stackrel{}{u}}(\tau )\right\}\frac{1}{3}\widehat{I}\dot{\stackrel{}{u}}(\tau )\right]+E_2_0^t\psi _2\left(t\tau \right)\widehat{I}\dot{\stackrel{}{u}}(\tau )$$
(72)
where the (dimensionless) functions $`\psi _1(t)`$ and $`\psi _2(t)`$ are relaxation (or “memory”) functions for the distortion strain and the dilatation, respectively. Note that eq. (72) coincides with the corresponding expression for the viscous stress tensor in pao ; impact  for $`\psi _2(t)=0`$. The latter approximation means that one neglects the bulk dissipation due to the dilatation rate. For the normal motion of colliding particles, however, the dissipation of energy due to the dilatation rate and the dissipation due to the distortion strain rate are of the same order of magnitude. Thus we keep both relaxation functions in our considerations. Introducing transversal and longitudinal velocities of sound in the material

$`c_t^2`$ $`=`$ $`{\displaystyle \frac{E_1}{2\rho }}={\displaystyle \frac{Y}{2\rho \left(1+\nu \right)}}`$ (73)
$`c_l^2`$ $`=`$ $`{\displaystyle \frac{2E_1+3E_2}{3\rho }}={\displaystyle \frac{Y\left(1\nu \right)}{\rho \left(1+\nu \right)\left(12\nu \right)}}`$ (74)
$`b^2`$ $`=`$ $`{\displaystyle \frac{c_l^2}{c_t^2}}={\displaystyle \frac{2\left(1\nu \right)}{\left(12\nu \right)}}`$
one can write the equation of motion for the viscoelastic medium

$$\frac{1}{c_t^2}\ddot{\stackrel{}{u}}=\left\{_{}^2\stackrel{}{u}+b^2_{}^2\stackrel{}{u}\right\}+_{}^2\psi _1\dot{\stackrel{}{u}}+_{}^2\left\{\frac{4}{3}\psi _1\dot{\stackrel{}{u}}+\left(b^2\frac{4}{3}\right)\psi _2\dot{\stackrel{}{u}}\right\}$$
(75)
where $`\psi \dot{\stackrel{}{u}}`$ denotes convolution.

To estimate the relative importance of the terms in eq. (75) we introduce the characteristic velocity $`v_0=g^N`$ and the characteristic time $`\tau _0=t_c`$, where $`t_c`$ is the duration of the collision, introduced above in eq. (56). Then the characteristic length is $`R_0=v_0\tau _0`$. Eq. (75) can be written in a dimensionless form

$$\left(\frac{v_0^2}{c_t^2}\right)\ddot{\stackrel{}{\stackrel{~}{u}}}\left(t\right)=\left\{_{}^2\stackrel{}{\stackrel{~}{u}}\left(t\right)+b^2_{}^2\stackrel{}{\stackrel{~}{u}}\left(t\right)\right\}+\left(\frac{\tau _{vis,1}}{\tau _0}\right)_{}^2\dot{\stackrel{}{\stackrel{~}{u}}}\left(\stackrel{~}{t}_{}\right)+\left(\frac{\tau _{vis,1}}{\tau _0}\right)_{}^2\left\{\frac{4}{3}\dot{\stackrel{}{\stackrel{~}{u}}}\left(\stackrel{~}{t}_{}\right)+\left(b^2\frac{4}{3}\right)\left(\frac{\tau _{vis,2}}{\tau _{vis,1}}\right)\dot{\stackrel{}{\stackrel{~}{u}}}\left(\stackrel{~}{t}_{}\right)\right\}$$
(76)
$$\tau _{vis,1/2}=_0^t\psi _{1/2}\left(\tau \right)𝑑\tau .$$
(77)
We use the following representation of the convolution

$$\psi \dot{\stackrel{}{u}}=\dot{\stackrel{}{\stackrel{~}{u}}}\left(\stackrel{~}{t}_{}\right)_0^t\psi _{1/2}\tau 𝑑\tau $$
(78)
with $`\stackrel{~}{t}_{}`$ being a dimensionless time from the interval $`0\stackrel{~}{t}_{}\stackrel{~}{t}=t/\tau _0`$. The relation for the convolution (78) is valid if $`\psi (t)0`$.

During the collision process $`t`$ is of the order of $`\tau _0`$, i.e. $`\stackrel{~}{t}`$ is of the order of $`1`$, while by the definition of $`\tau _{vis,1/2}`$ these values are of the order of the relaxation times for the dissipative processes in the material. That means, that $`\tau _{vis,1/2}`$ characterizes the time when the memory effects are important. If the duration of the collision is much greater than the relaxation times, i.e. if $`\tau _{vis,1/2}\tau _0`$, one can write

$$\tau _{vis,1/2}_0^{\mathrm{}}\psi _{1/2}\left(\tau \right)𝑑\tau $$
(79)
and consequently

$$\psi \dot{\stackrel{}{u}}\dot{\stackrel{}{u}}\left(t\right)\tau _{vis,1/2}.$$
(80)
If the characteristic velocity $`v_0`$ is much less then the speed of sound in the material too, one can neglect the terms with vanishing factors $`\left(v_0^2/c_t^2\right)`$ and $`\left(\frac{\tau _{vis,1/2}}{\tau _0}\right)`$ in eq. (76) and finally one arrives at the static eq. (45). That means the quasistatic approach is valid provided that the conditions

$`1`$ $``$ $`{\displaystyle \frac{v_0^2}{c_t^2}}`$ (81)
$`1`$ $``$ $`{\displaystyle \frac{\tau _{vis,1/2}}{\tau _0}}`$ (82)
hold. From the above considerations follows that in the quasistatic approximation the memory effects in the dissipative processes are not important and the viscous part of the stress tensor may be written in the same way as in eq. (II.1), with the viscous constants $`\eta _1`$ and $`\eta _2`$ given by

$$\eta _{1/2}=E_{1/2}\tau _{vis,1/2}=E_{1/2}_0^{\mathrm{}}\psi _{1/2}\left(\tau \right)𝑑\tau $$
(83)
It is worth to note that the quasistatic approximation is valid for many of the granular gases one encounters in nature, since usually the characteristic velocity in these systems is low. One should also note that the description of the collision in the quasistatic approximation is rigorous in a sense that no other additional approximations are used.

As it follows from the above considerations, it is not correct to use the time dependent relaxation functions for the dissipative part of the stress tensor together with Hertz’s quasistatic relations pao ; impact , since this approach is not self–consistent and one has to assume a lot of additional hardly controllable approximations.


format: text

text:
# Two-electron correlated motion due to Coulomb repulsion

## Abstract

A Hubbard-type model is derived from the microscopic Schrödinger equation. We found that additional terms describing direct two-electron transitions must be added to the standard Hubbard Hamiltonian. Such a Hamiltonian generates two-electron pairing due to on-site Coulomb repulsion. We demonstrate that the electron pairs with opposite spin propagate across periodic structures via direct and sequential two-electron tunneling. This mechanism can be used for a generation of entangled electron current. Numerical calculations show stability of the electron pairs.

PACS: 73.63.K, 03.65.X, 71.10.F

Within the framework of the Hubbard model it has been earlier shown that the on-side electron-electron interaction ($`U`$) might generate coherent propagation of the two electrons. However, in the Hubbard model such a process takes place only via sequential tunneling, i.e. by a virtual disintegration of a pair. As a result, the corresponding amplitude is suppressed with increasing $`U`$.

In fact, a more probable mechanism for two-particle tunneling transitions takes place due to simultaneous tunneling of two particles. Yet, this type of process is not accounted for by the standard Hubbard Hamiltonian. In this Letter we show how such a direct two-electron transitions can be included in the tunneling Hamiltonian by deriving the latter from the original Schrödinger equation. We also demonstrate that such a modified Hubbard Hamiltonian produces an effective two-electron coupling due to on-site electron-electron repulsion. In this case pairs of electron with opposite spins can propagate along periodic systems, producing a current of entangled electrons.

We start with an example of two electrons occupying a coupled-dot system, as shown schematically in Fig. 1. Each of the dots is represented by a square well potential, $`V_{1,2}(x)`$, which contain a bound state at the same energy $`E_0`$. If both electrons occupy one of the dots (Fig. 1a,c) the total energy of this system is $`2E_0+U`$. We assume that the inter-dot Coulomb repulsion is zero, so whenever the electrons are in different dots (Fig. 1b) the total energy becomes $`2E_0`$.

Such a system is usually described by the Hubbard tunneling Hamiltonian

$$H=\underset{i,s}{}\left(E_0n_{is}+\frac{U}{2}n_{is}n_{i\overline{s}}\right)\underset{s}{}\mathrm{\Omega }_0(a_{1s}^{}a_{2s}+\text{H.c.})$$
(1)
where $`a_{is}^{}`$ ($`a_{is}`$) creates (annihilates) an electron in the corresponding dot ($`i=1,2`$) with the spin $`s`$ ($`\overline{s}=s`$) and $`n_{is}=a_{is}^{}a_{is}`$. The amplitude $`\mathrm{\Omega }_0`$ generates transitions between the dots via single electron tunneling. This amplitude can be evaluated by using the Bardeen’s formula as a product of two bound state wave functions at a point $`x_0`$ inside the barrier ($`a<x_0<da`$):

$$\mathrm{\Omega }_0=(\kappa /m)\mathrm{\Phi }_1(x_0)\mathrm{\Phi }_2(x_0).$$
(2)
Here $`\mathrm{\Phi }_{1,2}(x)`$ are the single electron wave functions in the left and the right dot: $`(K+V_{1,2})\mathrm{\Phi }_{1,2}=E_0\mathrm{\Phi }_{1,2}`$ with $`K`$ is the kinetic energy operator, and $`\kappa =(2m|E_0|)^{1/2}`$. Note that the value of $`\mathrm{\Omega }_0`$ is weakly dependent on $`x_0`$. For instance, for square well potentials (Fig. 1) and $`a0`$, the wave functions $`\mathrm{\Phi }_1(x)=\kappa ^{1/2}\mathrm{exp}\left(\kappa |x|\right)`$ and $`\mathrm{\Phi }_2(x)=\kappa ^{1/2}\mathrm{exp}\left(\kappa |dx|\right)`$, so that $`\mathrm{\Omega }_0=(\kappa ^2/m)\mathrm{exp}(\kappa d)`$.

The wave function $`|\mathrm{\Psi }(t)`$ describing a motion of two electrons inside the double-dot can be represented by

$$|\mathrm{\Psi }(t)=\underset{i,i^{}=1,2}{}b_{ii^{}}(t)a_i^{}a_i^{}^{}|0,$$
(3)
where $`b_{11}(t),b_{22}(t)`$ are the probability amplitudes of finding two electrons in the first or in the second dot respectively (Fig. 1a,c). The probability amplitude $`b_{12}(t)`$ corresponds to the configuration shown in Fig. 1b, and $`b_{21}(t)`$ is the same but with an interchange of the spins. Substituting Eq. (3) into the Schrödinger equation $`i_t|\mathrm{\Psi }(t)=H|\mathrm{\Psi }(t)`$, we obtain the following coupled equations for the amplitudes $`b(t)`$

$`i\dot{b}_{11(22)}`$ $`=`$ $`(2E_0+U)b_{11(22)}+\mathrm{\Omega }_0(b_{12}+b_{12})`$ (5)
$`i\dot{b}_{12(21)}`$ $`=`$ $`2E_0b_{12(21)}+\mathrm{\Omega }_0(b_{11}+b_{22})`$ (6)
Let us solve these equations for the initial conditions $`b_{11}(0)=1`$, $`b_{12(21)}(0)=b_{22}(0)=0`$, Fig. 1a. Consider first the case of no Coulomb repulsion between the electrons, $`U=0`$. Then the probability of finding both electrons in the same dot, $`P_{ii}(t)=|b_{ii}(t)|^2`$ (Fig. 1a,c), is $`P_{11}(t)=\mathrm{cos}^4(\mathrm{\Omega }_0t)`$ and $`P_{22}(t)=\mathrm{sin}^4(\mathrm{\Omega }_0t)`$. Respectively, the probability of finding two electrons in different dots (Fig. 1b) $`P_{12}(t)=|b_{12}(t)|^2+|b_{21}(t)|^2`$ is $`P_{12}(t)=\mathrm{sin}^2(2\mathrm{\Omega }_0t)/2`$.

However, in the case of strong Coulomb repulsion, $`U4\mathrm{\Omega }_0`$, the time-behavior of two-electron system is very different. The solutions of Eqs. (3) depend on a small parameter $`\mathrm{\Omega }_0/U`$. To second order of this parameter we find $`P_{11}(t)=\mathrm{cos}^2(\omega t)`$, $`P_{22}(t)=\mathrm{sin}^2(\omega t)`$ with $`\omega =2\mathrm{\Omega }_0^2/U\mathrm{\Omega }_0`$ and $`P_{12}(t)=(2\mathrm{\Omega }_0^2/U^2)\mathrm{sin}^2(Ut/2)1`$. Thus two electrons oscillate together between two wells, but with much lower frequency than in the noninteracting case. The probability of finding two electrons in different wells (Fig. 1b) is strongly suppressed.

The above result can be understood in the following way. If a dot is occupied by two electrons, the tunneling of one of these electrons out takes place with a variation of the total energy, $`\mathrm{\Delta }E=U`$, Fig. 1. The probability of such a process is $`\mathrm{\Omega }_0^2/\mathrm{\Delta }E^2`$, and therefore it is strongly suppressed. This is analogous to the time-evolution of a single particle in a double-well potential. The particle remains in the same well where it was initially localized if the energy levels in the wells are misaligned. However if the levels are of the same energy, the particle oscillates between two wells (Rabi oscillations). In our case these are the two-electrons states, Fig. 1a,c, which have the same energy ($`2E_0+U`$). Therefore the two electrons oscillate between these states. Yet in the framework of the tunneling Hamiltonian Eq. (1) these oscillations are generated by a sequential tunneling (co-tunneling) via the intermediate state of Fig. 1b, so that the frequency of oscillations is suppressed by a factor $`\mathrm{\Omega }_0/U`$.

It follows from these arguments that the probability of two-electron tunneling between the states of Fig. 1a,c would be enhanced if both electrons tunnel simultaneously, keeping the energy of the system constant. These processes must be accounted in any tunneling Hamiltonian description by introducing the corresponding terms generating direct two-particle transitions. These terms are naturally arise whenever the Hubbard Hamiltonian is consistently derived from the Schrödinger equation. Such a derivation is presented below by using the Green’s function technique of the two-potential approach.

Consider the Schrödinger equation for two electrons in a double-well potential $`V(x)=V_1(x)+V_2(x)`$, Fig. 1. For simplicity we disregard the spin by treating two electrons with opposite spin as distinguishable particles. The total Hamiltonian can be written as $`H=H_1+H_2+U_C`$ with $`H_i=K(x_i)+V(x_i)`$ for $`i=1,2`$, and $`U_C(x_1x_2)`$ is the Coulomb repulsive potential. We assume that $`U_C=U`$ for $`|r_1r_2|\overline{c}`$ and drops down outside this region, where $`a\overline{c}d`$.

We start with two electrons localized in the first dot, $`|\mathrm{\Psi }(0)=|\mathrm{\Phi }_{11}`$, where $`x_1,x_2|\mathrm{\Phi }_{11}=\mathrm{\Phi }_1(x_1)\mathrm{\Phi }_1(x_2)`$. To determine the time-development of this system we apply the Laplace transform $`|\stackrel{~}{\mathrm{\Psi }}(E)=_0^{\mathrm{}}|\mathrm{\Psi }(t)\mathrm{exp}(iEt)𝑑t`$. The Schrödinger equation then reads

$$(EH)|\stackrel{~}{\mathrm{\Psi }}(E)=i|\mathrm{\Psi }(0)$$
(7)
Using the Green’s function technique, we can rewrite the total wave function as

$`|\stackrel{~}{\mathrm{\Psi }}(E)=[1+\stackrel{~}{G}(E)\overline{V}_2]\stackrel{~}{b}_{11}(E)|\mathrm{\Phi }_{11}`$ (8)
$`\stackrel{~}{G}(E)=G_1(E)(1\mathrm{\Lambda }_{11})[1+\overline{V}_2\stackrel{~}{G}(E)],`$ (9)
where $`G_1(E)=(E\overline{K}\overline{V}_1)^1`$ with $`\overline{K}=K(x_1)+K(x_2)`$, $`\overline{V}_{1,2}=V_{1,2}(x_1)+V_{1,2}(x_2)+U_C^{(1,2)}(x_1x_2)`$ and $`\mathrm{\Lambda }_{11}=|\mathrm{\Phi }_{11}\mathrm{\Phi }_{11}|`$ is the projection operator. Here we represented $`U_C=U_C^{(1)}+U_C^{(2)}`$, where $`U_C^{(1)}=U_C`$ for $`x_1,x_2\overline{c}`$ and drops to zero outside this region. Correspondingly, $`U_C^{(2)}=U_C`$ for $`dx_1,dx_2\overline{c}`$ and drops to zero outside this region. The amplitude $`\stackrel{~}{b}_{11}(E)`$ in Eq. (8) is the Laplace transform of the probability amplitude $`b_{11}(t)`$ of finding both electrons in the first dot, Fig. 1a. One obtains from Eqs. (7)-(9)

$$\stackrel{~}{b}_{11}(E)=\frac{i}{E2E_0U\mathrm{\Phi }_{11}|\overline{V}_2+\overline{V}_2\stackrel{~}{G}(E)V_2|\mathrm{\Phi }_{11}}$$
(10)
Eq. (8) can be treated iteratively by expanding $`\stackrel{~}{G}`$, Eq. (9), in powers of $`G_1`$. Since $`E2E_0+U`$ the Green’s function $`G_1`$ is dominated by the bound state pole in its spectral representation, $`G_1(E)|\mathrm{\Phi }_{11}\mathrm{\Phi }_{11}|/(E2E_0U)`$ for $`E2E_0+U`$. Yet, this pole is excluded by the projection operator $`1\mathrm{\Lambda }_{11}`$. The remaining part of $`G_1`$ gives rise to the corresponding (Born) series for $`|\stackrel{~}{\mathrm{\Psi }}(E)`$ which, however, converges very slowly. We therefore look for a different expansion for $`\stackrel{~}{G}`$ which converges much faster than the Born series. This can be achieved by expanding $`\stackrel{~}{G}`$ in powers of the Green’s function $`G_2(E)=1/(E\overline{K}\overline{V}_2)`$ by using the relation

$$\stackrel{~}{G}=G_2(1+\overline{V}_1\stackrel{~}{G})G_2\mathrm{\Lambda }_{11}(1+\overline{V}_2\stackrel{~}{G}).$$
(11)
Since the second dot contains the two-electron bound state at the same energy $`2E_0+U`$ as the first dot, the Green’s function $`G_2(E)`$ can be replaced by $`G_2(E)|\mathrm{\Phi }_{22}\mathrm{\Phi }_{22}|/(E2E_0U)`$ for $`E2E_0+U`$, where $`x_1,x_2|\mathrm{\Phi }_{22}=\mathrm{\Phi }_2(x_1)\mathrm{\Phi }_2(x_2)`$. Substituting this expression into Eq. (11) we find that in the contrast to Eq. (9) the bound state pole is not cancelled by the projection operator $`1\mathrm{\Lambda }_{11}`$. In addition, the second term in Eq. (11) is proportional to the overlap $`\mathrm{\Phi }_{22}|\mathrm{\Phi }_{11}`$ and can be considered as a small perturbation. If we neglect this term Eq. (11) can be easily solved thus obtaining

$$\stackrel{~}{G}(E)\frac{|\mathrm{\Phi }_{22}\mathrm{\Phi }_{22}|}{E2E_0U2\delta _{12}},$$
(12)
where $`\delta _{12}=\mathrm{\Phi }_2|V_1|\mathrm{\Phi }_2`$ is a small energy shift. Note that the Green’s function $`\stackrel{~}{G}(E)`$ contains also the poles at $`E2E_0`$, corresponding to separated electrons, Fig. 1b. Since $`E2E_0+U`$ these poles do not affect $`\stackrel{~}{G}(E)`$ for large $`U`$. Substituting Eq. (12) into Eq. (8) we find

$$|\stackrel{~}{\mathrm{\Psi }}(E)=\stackrel{~}{b}_{11}(E)|\mathrm{\Phi }_{11}+\stackrel{~}{b}_{22}(E)|\mathrm{\Phi }_{22},$$
(13)
where

$$\stackrel{~}{b}_{22}(E)=\frac{\mathrm{\Omega }_2}{E2E_0U2\delta _{12}}\stackrel{~}{b}_{11}(E)$$
(14)
is the Laplace transform of the probability amplitude $`b_{22}(t)`$, Eq. (3), and $`\mathrm{\Omega }_2=\mathrm{\Phi }_{22}|\overline{V}_2|\mathrm{\Phi }_{11}`$. Then substituting Eq. (12) into Eq. (10) and using Eq. (14) we obtain

$$(E2E_0U2\delta _{21})\stackrel{~}{b}_{11}(E)+\mathrm{\Omega }_2\stackrel{~}{b}_{22}(E)=i,$$
(15)
where $`\delta _{21}=\mathrm{\Phi }_1|V_2|\mathrm{\Phi }_1`$ (in our case $`\delta _{12}=\delta _{21}=\delta `$).

Performing the inverse Laplace transform, Eqs. (14)-(15) can be rewritten in a matrix form as

$$i\left(\begin{array}{c}\dot{b}_{11}(t)\hfill \\ \dot{b}_{22}(t)\hfill \end{array}\right)=\left(\begin{array}{cc}2E_0^{}+U& \mathrm{\Omega }_2\\ \mathrm{\Omega }_2& 2E_0^{}+U\end{array}\right)\left(\begin{array}{c}b_{11}(t)\hfill \\ b_{22}(t)\hfill \end{array}\right),$$
(16)
where $`E_0^{}=E_0+\delta `$. Using the occupation number representation and re-introducing the spin variable one finds that Eq. (16) can be rewritten as the Schrödinger equation $`i_t|\mathrm{\Psi }(t)=H|\mathrm{\Psi }(t)`$, where $`|\mathrm{\Psi }(t)`$ is given by Eq. (3) and the Hamiltonian

$`H={\displaystyle \underset{i=1,2}{}}[E_0^{}(a_i^{}a_i`$ $`+`$ $`a_i^{}a_i)+Ua_i^{}a_ia_i^{}a_i]`$ (17)
$``$ $`\mathrm{\Omega }_2(a_1^{}a_1^{}a_2a_2+H.c.).`$ (18)
Here the amplitude $`\mathrm{\Omega }_2=\mathrm{\Phi }_{22}|\overline{V}_2|\mathrm{\Phi }_{11}`$ generates a direct coupling in 2-dimensional Hilbert space between two-electron states, Fig. (1a)$``$Fig. (1c). This coupling can be calculated directly by using $`\mathrm{\Phi }_{22}|\overline{V}_2=\mathrm{\Phi }_{22}|(2E_0\overline{K})`$ and integrating by parts. As a result we find the following simple expression (c.f.).

$$\mathrm{\Omega }_2=\frac{2\kappa }{m}\mathrm{\Phi }_1(x_0)\mathrm{\Phi }_2(x_0)\mathrm{\Phi }_1(x^{})\mathrm{\Phi }_2(x^{})𝑑x^{}.$$
(19)
For the square-well potentials (Fig. 1) and $`a0`$, Eq. (19) yields $`\mathrm{\Omega }_2=(\kappa d/|E_0|)\mathrm{\Omega }_0^2`$, where $`\mathrm{\Omega }_0`$ is single-electron tunneling amplitude, Eq. (2). As expected the direct two-electron transitions, generated by $`\mathrm{\Omega }_2`$, dominate over the sequential two-electron transitions, $`2\mathrm{\Omega }_0^2/U`$, with an increase of $`U`$. This dominance of course is not limited to square-well potentials.

Note that Eq. (19) has been derived for long range repulsive potential $`U_C`$. For short range potential such that $`U_C(x_1x_2)0`$ for $`|x_1x_2|a`$, the wave function $`\mathrm{\Phi }_{11}`$ and $`\mathrm{\Phi }_{22}`$ are not given by a product of single electron wave functions. Yet the amplitude $`\mathrm{\Omega }_2`$ can be evaluated by using the Bardeen formula.

Comparing Eq. (18) with Eq. (1) we find that the term $`\mathrm{\Omega }_0a_{1s}^{}a_{2s}`$ describing single electron transitions does not arise in Eq. (18) since the corresponding poles of $`\stackrel{~}{G}(E)`$ at $`E2E_0`$ were neglected in our derivation. However with a decrease of $`U`$, such that $`\mathrm{\Omega }_22\mathrm{\Omega }_0^2/U`$, these poles must be taken into account. As a result we arrive to the Hamiltonian (1) supplemented with the last term of Eq. (18). Solving the Schrödinger equation with such a modified Hamiltonian for the initial conditions corresponding to Fig. 1a we find in the limit of $`U\mathrm{\Omega }_0`$ that $`P_{11}(t)=\mathrm{cos}^2(\omega ^{}t)`$, $`P_{22}(t)=\mathrm{sin}^2(\omega ^{}t)`$ and $`P_{12}(t)=(2\mathrm{\Omega }_0^2/U^2)\mathrm{sin}^2(Ut/2)`$, where $`\omega ^{}=|\mathrm{\Omega }_2(2\mathrm{\Omega }_0^2/U)|`$. Note that the transition amplitudes $`\mathrm{\Omega }_0`$ and $`\mathrm{\Omega }_2`$, given by Eqs. (2 ) and (19) change the sign if the wave functions $`\mathrm{\Phi }_{1,2}(x)`$ are of a different parity.

Our procedure can be extended to a general case of $`N`$ coupled wells. Detailed derivation will be presented elsewhere. Here we give our final result, representing the Hubbard-type Hamiltonian with additional terms for direct two-electron transitions

$`H={\displaystyle \underset{i,s}{\overset{N}{}}}`$ $`(E_0^{}n_{is}+{\displaystyle \frac{U}{2}}n_{is}n_{i\overline{s}}){\displaystyle \underset{i,s}{\overset{N1}{}}}(\mathrm{\Omega }_0a_{is}^{}a_{i+1,s}`$ (21)
$`+{\displaystyle \frac{\mathrm{\Omega }_2}{2}}a_{is}^{}a_{i\overline{s}}^{}a_{i+1,s}a_{i+1,\overline{s}}+\text{H.c.}).`$
Here $`i=1,\mathrm{},N`$ and $`E_0^{}=E_0+\delta `$, where $`\delta `$ is a small energy shift defined in the same way as in Eq. (16). The two-electron coupling $`\mathrm{\Omega }_2`$ is given by Eq. (19) for $`U\mathrm{\Omega }_0`$, but $`\mathrm{\Omega }_20`$ for $`U0`$. In this limit, however, the contribution from direct two-electron transitions is suppressed with respect to that generated by single electron transitions. Thus Eq. (21) with $`\mathrm{\Omega }_2`$ given by Eq. (19) can be used for any values of $`U`$.

As in the previous case of a double-well potential, an electron pair with opposite spins initially localized in one of the wells cannot be separated in the limit of $`U\mathrm{}`$. Such a pair can only move as a whole object due to direct two-electron tunneling between neighboring wells (the amplitude of co-tunneling vanishes in this limit). As a result a mini-band of the width $`2\mathrm{\Omega }_2`$ appears, providing a current of entangled electrons.

If $`\mathrm{\Omega }_22\mathrm{\Omega }_0^2/U`$, the propagation of an electron pair is govern by the both co-tunneling and direct two-electron transitions. In order to assess their importance with respect to the single electron transport we solved numerically the Schrödinger equation $`i_t|\mathrm{\Psi }(t))=H|\mathrm{\Psi }(t)`$ with the Hamiltonian (21) and the wave function given by Eq. (3) with $`i,i^{}=1,\mathrm{},N`$ for $`N=10`$ and $`U=10\mathrm{\Omega }_0`$ and the initial conditions corresponding to the first site occupied by two electrons. The results of our calculations are shown in Fig. 2. We plot there the probabilities of finding one and two electrons at the last site, $`P_N(t)`$ and $`P_{NN}(t)`$, for $`\mathrm{\Omega }_2=0.1\mathrm{\Omega }_0`$ and $`\mathrm{\Omega }_2=0`$. Although for the chosen values of parameters the value of $`\mathrm{\Omega }_2=0.1\mathrm{\Omega }_0`$ is smaller than that of the co-tunneling one ($`2\mathrm{\Omega }_0^2/U`$), the direct two-electron transitions strongly affects the pair transport. For instance the pair reaches the last site considerably faster and with much larger probability.

It follows from Fig. 2 that the probability of finding a single electron at the last site is much lower than that of an entangled electron pair. Yet, one cannot establish from this figure that this is true for all sites. Moreover, it is natural to assume that the electron pair decays very fast. Indeed, the weight of all states with two electrons are at different sites is much larger then those occupying the same site. In order to investigate this point we evaluated numerically the total probability of finding two electrons at the same site inside the chain, $`P_{pair}(t)=_i|b_{ii}(t)|^2`$. The results for the previous values of parameters are displayed in Fig. 3. It follows from this figure that contrary to the expectations the two-electron pair is stable during long time interval. It is interesting to note that the direct two-electron transitions affect the stability of this pair very little.

The entangled electron pair transport discussed above can be realized in multi-dot structures coupled to the reservoirs (emitter and collector) with a voltage bias larger than $`U`$. Then two electrons with opposite spins can enter simultaneously the first dot and then propagate across the entire structure. The single and the pair electron currents in the collector can be found in the most simple way by using the Bloch-type rate equations for the reduced density matrix derived in.

The two-electron coupling can also take place between electrons in neighboring dots via the inter-dot Coulomb repulsion. Indeed, if the separation energy of such a pair, $`U_{i,i+1}U_{i,i+2}`$ is larger then a single-electron transition amplitude, $`\mathrm{\Omega }_0`$, the pair moves across the system as a whole object. This phenomenon for three coupled dots has been studied in, but only taking into account the co-tunneling mechanism for the pair transport.

In conclusion, we demonstrated that the Coulomb repulsion can bound electrons in multi-well periodic structures, so that a current of entangled electrons with opposite spins would appear. The propagation of such electron pairs across the system is dominated by direct two-electron transitions absent in the standard Hubbard model. We expect that these two-electron transitions play an important role in different physical processes.

The author is indebted to M. Heiblum, Y. Imry, M. Kugler and A. Yacoby for fruitful discussions.