text:
# Pulse Dynamics in a Chain of Granules With Friction

## I Introduction

The propagation of pulses in granular materials has been a subject of vigorous recent interest. Seminal work on this subject was carried out by Nesterenko nesterenko  two decades ago in his analysis of the propagation of nonlinear compression pulses along a line of particles. This early work firmly established the nonlinear flavor of the problem: Nesterenko showed that under appropriate assumptions, among them the slow spatial variation of the displacements of the particles, the equations of motion for granular particles could in most cases be approximated by a continuous nonlinear partial differential equation that admits a soliton solution (later shown to actually be a solitary wave solution mackay ; senmanciu ) for a propagating perturbation in the chain. The recent revival of interest in the subject mackay ; senmanciu ; sinkovits ; rogers ; coste ; sen98 ; naughton ; hinch ; hascoet ; hascoet2 ; sen-book ; hong ; sen-pasi ; wu ; nakagawa  has been triggered in part by a concern with important technological applications such as the design of efficient shock absorbers sen-pasi , the detection of buried objects rogers ; coste ; sen98 ; naughton , and the fragmentation of granular chains hinch . The revival has involved advances in the modeling, simulation, and solution of the equations associated with important features of granular materials such as their discreteness mackay ; senmanciu ; hascoet ; sen-pasi , dimensionality sen-pasi , disorder sen98 ; hascoet2 ; sen-pasi , and loading provided by gravitational forces sinkovits ; sen98 ; hong ; sen-pasi ; sen96 ; hascoet3 . The preponderance of the work has been numerical, but important analytic advances have also been made, as well as experimental verifications of some of the theoretical predictions coste ; nakagawa .

The standard generic model potential between monodisperse elastic granules that repel upon overlap according to the Hertz law is given by hertz ; landau

$$\begin{array}{ccc}V(\delta _{k,k+1})\hfill & =a|\delta |_{k,k+1}^n,\hfill & \delta 0,\hfill \\ & & \\ V(\delta _{k,k+1})\hfill & =0,\hfill & \delta >0.\hfill \end{array}$$
(1)
Here

$$\delta _{k,k+1}2R\left[(z_{k+1}+y_{k+1})(z_k+y_k)\right],$$
(2)
$$a=\frac{2}{5D(Y,\sigma )}\left(\frac{R}{2}\right)^{1/2},D(Y,\sigma )=\frac{3}{2}\left(\frac{(1\sigma ^2)}{Y}\right),$$
(3)
$`Y`$ and $`\sigma `$ denote Young’s modulus and Poisson’s ratio, $`z_k`$ denotes the initial equilibrium position of granule $`k`$ in the chain, and $`y_k`$ is the displacement of granule $`k`$ from this equilibrium position. The geometric parameter $`R`$ is the radius if the particles are elastic spheres. More generally, $`R`$ is determined by the principal radii of curvature of the surfaces at the point of contact landau . The exponent $`n`$ is $`5/2`$ for spheres, it is $`2`$ for cylinders, and in general depends on geometry.

The displacement of the $`k`$-th granule ($`k=1,2,\mathrm{},L`$) in the chain from its equilibrium position in a frictional medium is governed by the equation of motion

$$m\frac{\mathrm{d}^2y_k}{\mathrm{d}\tau ^2}=\stackrel{~}{\gamma }\frac{\mathrm{d}y_k}{\mathrm{d}\tau }a(y_ky_{k+1})^{n1}\theta (y_ky_{k+1})+a(y_{k1}y_k)^{n1}\theta (y_{k1}y_k).$$
(4)
Here $`\theta (y)`$ is the Heavyside function, $`\theta (y)=1`$ for $`y>0`$, $`\theta (y)=0`$ for $`y<0`$, and $`\theta (0)=1/2`$. It ensures that the particles interact only when in contact. Note that for open boundary conditions the second term on the right hand side of this equation is absent for the last granule and the third term is absent for the first, while for periodic boundary conditions $`y_{k+L}=y_k`$. In terms of the rescaled variables and parameters

$$y_k=\left(\frac{mv_0^2}{a}\right)^{1/n}x_k,\tau =\frac{1}{v_0}\left(\frac{mv_0^2}{a}\right)^{1/n}t,\gamma =\frac{\stackrel{~}{\gamma }}{mv_0}\left(\frac{mv_0^2}{a}\right)^{1/n},$$
(5)
Eq. (4) can be rewritten as

$$\ddot{x}_k=\gamma \dot{x}_k(x_kx_{k+1})^{n1}\theta (x_kx_{k+1})+(x_{k1}x_k)^{n1}\theta (x_{k1}x_k),$$
(6)
where a dot denotes a derivative with respect to $`t`$.

Initially the granules are placed along a line so that they just touch their neighbors in their equilibrium positions (no precompression), and all but one granule, granule $`i`$, are at rest. The velocity of granule $`i`$ is $`v_0`$ (the impulse). In the rescaled variables the initial conditions become $`x_k(0)=\dot{x}_k(0)=0`$, $`ki`$, $`x_i(0)=0`$, and $`\dot{x}_i(0)=1`$. In our work we set $`i=1`$.

In the absence of the frictional contribution $`\stackrel{~}{\gamma }\mathrm{d}y_k/\mathrm{d}\tau `$ or $`\gamma \dot{x}_k`$, when $`n>2`$ an initial impulse settles into a pulse that becomes increasingly narrow with increasing $`n`$, and propagates at a velocity that is essentially constant and determined by $`n`$ and by the amplitude of the pulse. For $`n=2`$ the pulse spreads in time and travels at a constant velocity independent of pulse amplitude. In the latter case there is considerable backscattering that leads to backward motion of all the granules behind the pulse, whereas the backscattering is minimal for $`n>2`$ hinch . The pulse is a completely conservative solitary wave in the limit $`n\mathrm{}`$.

Our interest lies in ascertaining the effects of the frictional contributions on these results. Although friction is on occasion mentioned in theoretical and numerical work, it is usually mentioned as an element that is neglected or omitted. However, its presence and importance in experiments is recognized as inevitable nakagawa . Note that frictional effects may arise not only from the forces between the granular chain and the surrounding medium but also from the conversion of translational motion to other degrees of freedom (e.g. rotation) nakagawa . In an earlier paper we analyzed in detail the effects of frictional forces on the dynamics of two granules, specifically the way in which forward and backward motion of the granules is affected by these forces we . Herein we extend that work to a chain of granules.

Since our approach will in general follow established methods, it is useful to lay out at the outset an overview of the principal approaches that have been implemented in the study of pulse dynamics in frictionless chains, and the various features that determine these dynamics. Theoretical studies of pulse dynamics in frictionless chains have mainly relied on three approaches:

(1) Numerical solution of the equations of motion sinkovits ; sen98 ; hinch ; hascoet ; hascoet2 ; sen-pasi ; hascoet3 ;

(2) Continuum approximations to the equations of motion followed by exact or approximate solutions of these approximate equations nesterenko ; hinch ; hascoet3 ; and

(3) Phenomenology about properties of pairwise (or at times three-body) collisions together with the assumption that pulses are sufficiently narrow to be principally determined by these properties wu ; mohan .

From these studies it is clear that there three features determine the dynamics in these chains:

(1) The power $`n`$ in the potential;

(2) The absence of a restoring force; and

(3) The discreteness of the system.

While each of these leads to essential aspects of the dynamics, the literature is not always entirely clear on which feature is determinant in a particular behavior, nor is it always clear which of these features is being approximated or even ignored. One example is the equivocal connection between the dynamics of granular systems and systems with power law interaction potentials that include a restoring force. For instance, there is an extensive literature on the Fermi-Pasta-Ulam (FPU) chain with purely nonlinear interactions of the form appearing in Eq. (6) with $`n>2`$ but without the Heavyside theta functions FPU . Highly localized pulses propagate in these systems sen-pasi ; ourpulse , but their relation to the localized solutions in granular systems is by no means clear. This ambiguity is exacerbated by the fact that discrete systems are frequently approximated by continuum equations. While the pulse-like (soliton or solitary wave) solutions that emerge from these approximations are assumed (and even shown numerically) to describe certain aspects of granular system, the continuum approximations never explicitly respect the absence of a restoring force, i.e., they are more clearly justifiable for FPU-like systems. And so with these solutions in hand, it is not clear what consequences of the absence of a restoring force and of the discreteness of the system have been lost. In some sense opposite to the continuum approximations are two-body or three-body phenomenologies. These provide considerable qualitative insight into some consequences of discreteness, but are quantitatively not sufficiently accurate because the true dynamics extend over more than two or three granules. In this paper, even while focusing on the effects of friction, we make an attempt to provide some clarity on these issues.

Numerical simulations are of course extremely helpful and ultimately an accurate test of theory. Note that the rescaled system Eq. (6) has no free parameters in the absence of friction and only a single parameter when there is friction, and therefore numerical characterization is particularly straightforward in this particular system. However, this universality is lost with any of a number of modifications that might be interesting and have been considered recently such as, for example, tapering of the masses in the chain senmanciu ; wu ; nakagawa  and mass sen98 ; hascoet2  or frictional we  disorder. Consequently, it is desirable to have a strong analytic backdrop.

The paper is organized as follows. Section II deals with a chain of cylindrical granules. First we review existing results for frictionless granules, discussing them in the context of the issues and uncertainties noted earlier. We then extend the theory to granules subject to friction, and present numerical simulations in support of these results. In Sec. III we follow the same presentation sequence for granules with $`n>2`$, with special attention to spherical granules in our numerical simulations. A summary of results and of future directions is presented in Sec. IV.

## II Cylinders ($`n=2`$)

Cylindrical granules have provided a theoretical testbed for dynamics in granular chains because the exponent $`n=2`$ leads to analytic manageability not available for other potentials. Although sometimes called the “harmonic” case, it should be remembered that the Hertz potential is quite different because there is no restoring force, that is, the cylindrical Hertz potential is half of a harmonic potential. The derivative of the force law therefore changes discontinuously between extension and compression. This leads to considerable differences in the chain dynamics. Nevertheless, some aspects of the cylindrical chain dynamics can be inferred from those of a harmonic chain, and this can be exploited to great advantage in the analysis.

If one were to ignore the absence of a restoring force and implement the simplest continuum approximation to Eq. (6) in the absence of friction, the result would be a simple wave equation with diffusive coupling whose solutions do not represent the observed behavior of the $`n=2`$ chain. In reality, an initial impulse in the chain described by Eq. (6) in the absence of friction moves as a spreading pulse. Although the pulse spreads and sheds some energy, the waveform can nevertheless be clearly identified as a pulse hinch . Its maximum $`k_{max}(t)`$ travels forward with a constant unit velocity and a displacement amplitude that increases as $`t^{1/6}`$. The pulse spreads in time as $`t^{1/3}`$, more slowly than it would in a system with diffusive coupling. It is interesting to understand which of these features are due to the discrete nature of the chain, and which are due to the absence of a restoring force. Further, once these questions have been answered, it is interesting to explore what happens to these features in the presence of friction.

Our presentation of the chain of cylindrical granules consists of three parts. First we review the results of Hinch and Saint-Jean (HSJ) hinch  for a frictionless chain, and we recast the problem in a way that clarifies some of the approximations made in that work and pinpoints the sources of differences between the $`n=2`$ chain and a purely diffusive coupling. We supplement this review with our analytic results that reproduce some of their purely numerical ones. Then we modify this analysis to include the effects of weak hydrodynamic friction on the granules. Finally, we complement this analytic (and necessarily approximate) treatment with a comparison with numerical simulation results for the frictional chain of cylindrical granules.

### II.1 Frictionless Granules - Theory

The HSJ theory is based on the following approximations implemented consecutively and independently.

1. The solution is assumed to be described by a traveling pulse of constant form which propagates at constant unit speed and has an amplitude $`b`$ and width $`\lambda `$ that vary slowly with time:

$$x(k,t)=b(t)f\left(\frac{kt}{\lambda (t)}\right).$$
(7)
2. The pulse is assumed to retain almost all of its initial energy. HSJ show that assumption 1. leads to equipartition of this energy between potential and kinetic forms, as should be the case for a harmonic potential.

These two assumptions are sufficient to lead to the conclusion that $`\lambda b^2`$.

3. A continuum approximation is implemented that takes into account some discreteness effects. In the lowest order diffusive coupling approximation one would set the difference $`x_{k+1}+x_{k1}2x_k^2x(k)/k^2`$. Retention of the next term in a Taylor series expansion of $`x_{k\pm 1}`$ about $`x_k`$ leads to the continuum approximation that incorporates some of the effects of discreteness:

$$\frac{^2x(k,t)}{t^2}=\frac{^2x(k,t)}{k^2}+\frac{1}{12}\frac{^4x(k,t)}{k^4}.$$
(8)
Note that this expansion includes a restoring force. A transformation to a moving frame with unit propagation velocity, i.e., a change of variables from $`k`$ and $`t`$ to $`\nu =kt`$ and $`t`$, transforms this equation to

$$\frac{^2x}{t^2}\frac{^2x}{t\nu }=\frac{1}{12}\frac{^4x}{\nu ^4}.$$
(9)
Regardless of the form of $`b(t)`$ or of the function $`f`$, this is sufficient to establish that the solution Eq. (7) is consistent with this equation only if $`\lambda t^{1/3}`$, and also that the next term in the Taylor series expansion is unimportant for this result. The width of the pulse is therefore governed by the first manifestation of discreteness. One assumes that the absence of a restoring force does not affect this result because the granules within the pulse are in fact overlapping most of the time and hence subject mainly to the repulsive portion of the potential. Note that the continuum approximation together with the conservation of energy are then sufficient to arrive at the conclusion that $`b(t)t^{1/6}`$, i.e., that the pulse amplitude actually grows with time.

4. It is not yet clear that (7) is actually compatible with (8) until one determines the function $`f`$. Substitution of (7) into (8) and retention of leading terms in $`t`$ leads to the equation for $`f(\xi )`$,

$$f^{^{\prime \prime \prime \prime }}8\xi f^{^{\prime \prime }}4f^{^{}}=0.$$
(10)
This equation has four solutions. The one consistent with the requirement that $`f(\xi )`$ decays for large $`\xi `$ (i.e., ahead of the wave) and consistent with the assumption of energy conservation by the pulse is as

$$f(\xi )=N_\xi ^{\mathrm{}}Ai^2(2^{1/3}y)dy$$
(11)
where $`Ai(z)`$ is an Airy function and

$$N=\left[\frac{E_p}{_{\xi _0}^{\mathrm{}}Ai^4(2^{1/3}y)dy}\right]^{1/2}=3.533\sqrt{E_p}.$$
(12)
Here $`E_p`$ is the pulse energy and $`2^{1/3}\xi _0`$ is the first zero of $`Ai(z)`$.

These features describe the traveling displacement pulse of increasing width and amplitude quite accurately, as shown by the numerical simulations in hinch . We note that the total energy of the system with the initial unit velocity impulse at one granule is $`1/2`$. The numerical results of HSJ lead to an asymptotic pulse energy of $`E_p=0.48`$, that is, $`96.2\%`$ of the energy resides in the pulse. Below we calculate the pulse energy analytically.

In addition to the forward traveling pulse, HSJ also investigate the momentum of the particles ejected backward as the pulse goes by. That particles must be ejected is a consequence of the conservation of momentum: the traveling pulse of constant energy and increasing width and displacement amplitude carries increasing forward momentum, which must be balanced by the backward momentum of the ejected particles. One then arrives at the next item on the list of assumptions:

5. The absence of a restoring force is explicitly recognized in the calculation of the momentum of the ejected particles, which simply keep traveling backward with the momentum they acquire at the moment of separation from the pulse. Equating the rate of change of the momentum of the forward pulse at time $`t`$ to that of the particle ejected at that time leads to the conclusion that the backward momentum of the $`n^{th}`$ particle is $`\dot{x}_n=ct^{5/6}=cn^{5/6}`$, the latter equality arising from the unit speed of pulse propagation. The numerical simulations of HSJ lead to the value $`c=0.158`$, a value that we obtain analytically.

To obtain analytic results for the pulse energy $`E_p`$ and the constant of proportionality $`c`$, we note that the forward momentum of the propagating pulse is

$`P`$ $`=`$ $`{\displaystyle \underset{\dot{x}(k,t)>0}{}}\dot{x}(k,t)={\displaystyle \underset{\dot{x}(k,t)>0}{}}{\displaystyle \frac{b(t)}{\lambda (t)}}f^{}(\xi )`$ (13)
$`=`$ $`b(t){\displaystyle _{\xi _0}^{\mathrm{}}}f^{}(\xi )d\xi =\stackrel{~}{N}b(t)`$
where

$$\stackrel{~}{N}=N_{\xi _0}^{\mathrm{}}Ai^2(2^{1/3}\xi )d\xi =1.379\sqrt{E_p}.$$
(14)
The rate of change of the forward momentum then is $`\dot{P}=\stackrel{~}{N}\dot{b}(t)`$. Since the pulse velocity is unity, this is the momentum transferred to the last particle as it is ejected:

$$v_b=\dot{P}=\stackrel{~}{N}\dot{b}(t)=\frac{\stackrel{~}{N}}{6}t^{5/6}=\frac{\stackrel{~}{N}}{6}n^{5/6}.$$
(15)
The backscattered energy then is

$$E_b=\frac{1}{2}\left(\frac{\stackrel{~}{N}}{6}\right)^2\underset{k=1}{\overset{\mathrm{}}{}}k^{5/3}=0.056E_p$$
(16)
and for the total energy we obtain

$$E=E_p+E_b=1.056E_p,$$
(17)
from which it immediately follows that asymptotically the energy in the pulse is

$$E_p=0.947E.$$
(18)
Thus, 94.7% of the energy resides in the pulse (to be compared to the HSJ value of $`96.2\%`$ obtained from simulations) and the remainder is backscattered. With the initial condition $`E=1/2`$ used in all simulations we thus have for the constants defined earlier $`N=2.431`$ and $`\stackrel{~}{N}=0.949`$. The constant $`c`$ in HSJ is $`c=\stackrel{~}{N}/6=0.158`$, exactly as they obtained from numerical simulations.

This essentially completes the solution. The summary description is then that an initial velocity impulse propagates forward at unit speed, with a width that grows as $`t^{1/3}`$ and a displacement amplitude that grows as $`t^{1/6}`$. This pulse carries almost all of the initial energy and its momentum increases. This increase in momentum is compensated by granules that are ejected backward as the pulse passes. The speed of the ejected granules decreases, the $`n^{th}`$ granule being ejected with a speed proportional to $`n^{5/6}`$.

Further insights can be gained by viewing this problem a bit differently. Suppose that we do not implement a continuum approximation at all, but instead focus on the full discrete equation

$$\ddot{x}_k=x_{k+1}+x_{k1}2x_k.$$
(19)
Like the continuum equation, this of course includes restoring forces. This equation can be solved exactly ourpulse ,

$$x_k(t)=_0^tJ_{2k2}(2t^{})𝑑t^{},$$
(20)
where $`J_n(z)`$ is the Bessel function of the first kind as . The solution is a moving spreading pulse along with oscillatory displacements that are left in its wake (precisely because there are restoring forces). The peak and width of the pulse can be obtained from the properties of the Bessel functions, in particular from knowledge of their first two zeroes and the first maximum royal . The maximum of the pulse occurs at $`k_{max}(t)=t+𝒪(t^{1/3})`$, so the pulse velocity is unity. Its width increases as $`t^{1/3}`$. In these features the solution is similar to the one assumed in Eq. (7). However, the amplitude of the displacement pulse does not grow but is instead constant in time. (If we write the Bessel function solution in the form (7) but with $`b(t)=const`$ independent of $`t`$, we find from substitution into Eq. (8) that $`f`$ satisfies the equation $`f^{^{\prime \prime \prime \prime }}8zf^{^{\prime \prime }}8f^{^{}}=0`$.) The pulse energy according to this description is not constant but instead decreases as $`t^{1/3}`$. The energy that is lost goes into the oscillatory displacements in the wake of the harmonic pulse caused by the restoring forces. The additional and appropriate assumption of energy conservation in HSJ adds this lost energy back into the pulse without affecting its spreading rate or velocity. This indicates that (in the $`n=2`$ problem) energy conservation is an appropriate additional assumption in lieu of the absence of a restoring force.

In summary: the fact that the solution of the $`n=2`$ Hertz problem is a spreading pulse of unit velocity and of width that increases as $`t^{1/3}`$ is purely a consequence of the discreteness of the system and not dependent on the presence or absence of restoring forces; it is a feature of a discrete harmonic system. The (near) conservation of energy in the pulse is not a feature of a harmonic system and must be included as an additional assumption. Conservation of energy must be implemented not only to describe all the features of the traveling pulse but also to describe the backward momentum of the particles ejected as the pulse moves along. Extensive numerical results quantitatively supporting the features just described can be found in hinch .

### II.2 Granules With Friction - Theory

Here we generalize the previous theories to a chain in which the granules experience friction \[cf. Eq. (6)\]. We generalize the HSJ theory to this case and also use the solution of the damped harmonic chain to complement these results.

1. The solution is assumed to be a traveling pulse of constant shape which propagates at unit speed, has a width $`\lambda `$ that varies slowly in time, and an amplitude that aside from an exponential decay due to the friction also varies slowly with time:

$$x(k,t)=e^{\gamma t/2}b(t)f\left(\frac{kt}{\lambda (t)}\right).$$
(21)
We retain only the $`\gamma `$-independent portions of the slowly varying functions $`\lambda `$ and $`b`$. This is thus a “lowest order” ansatz. Indeed, we will show later that corrections to this lowest order ansatz seem to be of $`𝒪(\gamma ^3)`$. The pulse energy then decays exponentially as $`e^{\gamma t}`$. Explicitly, the kinetic energy is obtained as follows. Retaining only leading orders in time, the pulse velocity that follows from Eq. (21) is

$$\dot{x}(k,t)=e^{\gamma t/2}\frac{b(t)}{\lambda (t)}f^{}\left(\frac{kt}{\lambda (t)}\right)+𝒪(\gamma t).$$
(22)
Therefore, the kinetic energy to $`𝒪(\gamma t)`$ is

$`K(t)`$ $`=`$ $`e^{\gamma t}{\displaystyle \frac{1}{2}}{\displaystyle \underset{k}{}}\left[{\displaystyle \frac{b(t)}{\lambda (t)}}f^{}\left({\displaystyle \frac{kt}{\lambda (t)}}\right)\right]^2e^{\gamma t}{\displaystyle \frac{1}{2}}{\displaystyle \frac{b^2(t)}{\lambda (t)}}{\displaystyle _{\xi _0}^{\mathrm{}}}\left[f^{}(\xi )\right]^2d\xi `$ (23)
$`=`$ $`{\displaystyle \frac{E_p}{2}}e^{\gamma t}{\displaystyle \frac{b^2(t)}{\lambda (t)}},`$
where $`E_p`$ is the undamped pulse energy. On the other hand, the potential energy may be written as

$$U(t)=\frac{1}{2}\underset{k}{}(x_{k+1}x_k)^2\frac{1}{2}\underset{k}{}\left[\frac{x(k,t)}{k}\right]^2,$$
(24)
which also leads to

$$U(t)=\frac{E_p}{2}e^{\gamma t}\frac{b^2(t)}{\lambda (t)}$$
(25)
to $`𝒪(\gamma t)`$. The total pulse energy to this order then is

$$E_p(t)=E_pe^{\gamma t}\frac{b^2(t)}{\lambda (t)}.$$
(26)
2. We assume that almost all of the energy resides in the decaying pulse.

These results and assumptions are thus again sufficient to conclude that $`\lambda b^2`$. It then follows from Eq. (26) that the only effect of friction to this order is thus the overall exponential decay of the energy.

3. As before, a continuum approximation is implemented that takes into account some discreteness effects and of course the frictional contribution (while still including the nonexistent restoring force, as before). For this purpose it is convenient to implement the change of variables $`z(k,t)=e^{\gamma t/2}x(k,t)`$. Retention of the next term in a Taylor expansion of $`z_{k\pm 1}`$ about $`z_k`$ beyond the purely diffusive approximation leads to

$$\frac{^2z(k,t)}{t^2}=\frac{^2z(k,t)}{k^2}+\frac{1}{12}\frac{^4z(k,t)}{k^4}+\frac{\gamma ^2}{4}z(k,t).$$
(27)
The associated moving frame equation ($`\nu =kt`$) is

$$\frac{^2z}{t^2}\frac{^2z}{t\nu }=\frac{1}{12}\frac{^4z}{\nu ^4}+\frac{\gamma ^2}{4}z.$$
(28)
With the ansatz

$$z(k,t)=b(t)f\left(\frac{kt}{\lambda (t)}\right)$$
(29)
associated with (21) this is again sufficient to establish that the pulse widens as $`\lambda t^{1/3}`$. Corrections begin to set in when the last term in Eq. (28) becomes important, i.e., for times $`t\gamma ^{3/2}`$. For small $`\gamma `$ this is a time much longer than the entire lifetime of the pulse, which is of $`𝒪(\gamma ^1)`$ due to the overall exponential decay. As before, we are again led to the conclusion that $`b(t)t^{1/6}`$.

4. The envelope function is again the solution of Eq. (10) with corrections of $`𝒪(\gamma ^2)`$.

5. As before, as the pulse travels forward granules are ejected backwards. In the frictionless problem these granules continue moving backward with the same velocity forever after ejection because they are simply freely moving granules. In the problem with friction these granules progressively slow down because they, too, are damped, but it is of interest to calculate their momentum at the moment of ejection from the pulse. Indeed, we now show that in the presence of damping the backward ejection momentum is greater than in the undamped chain, a result consistent with that found in a two-granule system we .

The forward momentum of the propagating pulse to $`𝒪(\gamma t)`$ is

$`P`$ $`=`$ $`{\displaystyle \underset{\dot{x}(k,t)>0}{}}\dot{x}(k,t)={\displaystyle \underset{\dot{x}(k,t)>0}{}}e^{\gamma t/2}{\displaystyle \frac{b(t)}{\lambda (t)}}f^{}(\xi )`$ (30)
$`=`$ $`e^{\gamma t/2}b(t){\displaystyle _{\xi _0}^{\mathrm{}}}f^{}(\xi )d\xi =\stackrel{~}{N}e^{\gamma t/2}b(t),`$
where $`\stackrel{~}{N}`$ is given in Eq. (14). The rate rate of change of the forward momentum is

$$\dot{P}=\stackrel{~}{N}\dot{b}(t)e^{\gamma t/2}\left(1\frac{\gamma b(t)}{2\dot{b}(t)}\right).$$
(31)
The last particle in the pulse loses momentum $`\delta `$ through dissipation as the pulse moves across it. This momentum loss is the difference in the momentum of this particle when it is in the middle of the pulse (i.e., when it is maximally compressed, at which point its momentum is zero) and when it is at the end of the pulse and about to be ejected (at which point its compression is zero):

$$\delta =\gamma \mathrm{\Delta }x=\gamma x_{max}=\gamma \stackrel{~}{N}e^{\gamma t/2}b(t).$$
(32)
Since the pulse velocity is still unity, the momentum transferred to the last particle as it is ejected is

$`v_b`$ $`=`$ $`\dot{P}+\delta =\stackrel{~}{N}e^{\gamma t/2}\dot{b}(t)\left(1{\displaystyle \frac{\gamma b(t)}{2\dot{b}(t)}}+{\displaystyle \frac{\gamma b(t)}{\dot{b}(t)}}\right)`$ (33)
$`=`$ $`{\displaystyle \frac{\stackrel{~}{N}}{6}}e^{\gamma t/2}t^{5/6}\left(1+3\gamma t\right).`$
To facilitate comparison of this result with that of the frictionless chain we expand the exponential and retain terms of $`𝒪(\gamma t)`$:

$$v_b=\frac{\stackrel{~}{N}}{6}t^{5/6}\left(1+\frac{5}{2}\gamma t\right).$$
(34)
This confirms that the backward ejection momentum is indeed greater than in the undamped chain.

The forward moving pulse of the $`n=2`$ Hertz problem with small friction is essentially identical in shape to that of the frictionless problem except for an overall exponential decay. The pulse travels at unit speed and its width increases as $`t^{1/3}`$, these two features again being a consequence of discreteness and essentially unaffected by friction. The assumption that the pulse energy decreases only because of the friction is an additional assumption. This feature must again be implemented separately to describe the amplitude of the traveling pulse correctly, and also to describe the backward momentum of the granules ejected as the pulse moves along.

An analysis of the full discrete equation (with restoring forces) starts from the linear equation

$$\ddot{x}_k=x_{k+1}+x_{k1}2x_k\gamma \dot{x}_k.$$
(35)
The change of variables $`z_k(t)=e^{\gamma t/2}x_k(t)`$ immediately leads to

$$\ddot{z}_k=z_{k+1}+z_{k1}2z_k+\frac{\gamma ^2}{4}z_k$$
(36)
and consequently to the solution

$$x_k(t)=e^{\gamma t/2}\left[_0^tJ_{2k2}(2t^{})𝑑t^{}+𝒪(\gamma ^2)\right].$$
(37)
The solution is a moving spreading decaying pulse along with decaying oscillatory displacements. As before, the peak of the pulse occurs at $`k_{max}(t)=t+𝒪(t^{1/3})`$, its width increases as $`t^{1/3}`$, and the amplitude of the pulse decays exponentially.

In summary: the fact that the solution of the $`n=2`$ Hertz problem with friction is an exponentially decaying spreading pulse of unit velocity and of width that increases as $`t^{1/3}`$ is purely a consequence of the discreteness of the system and not dependent on the presence or absence of restoring forces; it is a feature of a discrete harmonic system. The assumption that the pulse retains almost all of the system energy (which decays exponentially) is again not a feature of a harmonic system and must be included as an additional assumption. This conservation feature must be implemented not only to describe correctly all the features of the traveling pulse but also to describe the backward momentum of the particles ejected as the pulse moves along. This backward momentum is greater in the chain with friction than in the frictionless system.

### II.3 Numerical Simulations

In this section we check the accuracy of our analytic results for the $`n=2`$ frictional chain through numerical simulations. First, we note that the total energy of the system decreases exponentially, as seen in Fig. 1 for various values of the friction coefficient. For times $`t10\gamma ^1`$, aside from the overall exponential decay of the energy the system behavior is indeed the same as that of the frictionless system in that it broadens as $`t^{1/3}`$ and propagates at unit velocity. We find that the pulse ceases to exist at a time of order $`t10\gamma ^1`$. More specifically, we find that at a time $`8.6\gamma ^1`$ the pulse energy decreases abruptly (more rapidly than exponential), and at a time $`15.7\gamma ^1`$ the backscattered energy becomes greater than the pulse energy.

We must check our prediction for the self-similar impulse wave $`f(\xi )`$ propagating along the chain of particles. We follow HSJ and plot the scaled velocity as a function of the scaled position for different times in Fig. 2. From Eq. (22) the appropriate scaled velocity to leading order in the friction is $`f^{}(\xi )=\dot{x}\lambda e^{\gamma t/2}/b\dot{x}t^{1/6}e^{\gamma t/2}`$, and this is the ordinate in the panels. The abcissa is $`\xi =(kt0.25)/\lambda (kt0.25)t^{1/3}`$. The shift $`0.25`$, which also occurs in the frictionless case, is carefully explained and derived in HSJ and comes about because the scaling function solution has velocities in the propagating impulse wave of $`𝒪(t^{1/6})`$ (here modified by the exponential friction factor) while the rebound velocity of particles is of $`𝒪(t^{5/6})`$ (again modified by the exponential factor). This indicates that there must be a correction of $`𝒪(t^{2/3})`$ to the leading-order term, and this correction appears as a shift in the scaling. The solid curve is $`f^{}`$ as obtained from Eq. (11), while the points are the results of our numerical simulations. The agreement is very good provided the damping is small, but serious deviations begin to set in with increasing damping, as seen in panel (d).

We have not calculated further corrections analytically, but a numerical test is possible. One might test a correction of $`𝒪(\gamma )`$ \[cf. Eq. (22)\], e.g. of the form $`b(t)=t^{1/6}\left[1+\gamma C(\nu ,t)\right]`$, or of $`𝒪(\gamma ^2)`$ \[cf. Eq. (37)\]. The coefficient $`C`$ in the correction would in general be a function of $`\nu `$. We have tested various corrections with the (admittedly unjustified) assumption that $`C`$ does not depend on $`\nu `$ (the numerical effort to do otherwise seems unwarranted), and have found that the leading correction seems to be of $`𝒪(\gamma ^3)`$, specifically $`b(t)=t^{1/6}(1+27\gamma ^3t)`$. This result is analytically appealing but is purely a numerical outcome. The rescaled velocity results are shown in Fig. 3.

Finally, we have predicted a somewhat unexpected feature of the backscattered granules, namely, that at the moment of ejection their velocity is greater than that of the corresponding ejected granules in the absence of friction \[cf. Eq. (34)\]. This is seen in Fig. 4. In the first panel we show the velocity of the last ejected granule as a function of time for different values of the friction coefficient. The velocity is indeed higher with increasing friction, up to a time beyond which more complex behavior not captured by our lowest order theory sets in. In the second panel we have scaled the velocities by the factor $`(1+\frac{5}{2}\gamma t)`$ which, according to our theory, should collapse the curves. At very early times, while the collapse occurs with the scaling predicted by our theory, the resulting curve is not yet the theoretical one because the pulse is not yet clearly defined. However, after this early period we see both the collapse as well as agreement with the theoretical line up to times at which higher order effects set in, indicating that our theory captures the correct behavior up to those times.

## III Spheres ($`n=5/2`$)

The archetypal system for the study of impulse dynamics associated with a Hertz potential of exponent $`n>2`$ consists of spherical granules. Such a nonlinear potential gives rise to a traveling pulse of essentially constant speed that seems to retain its shape for very long times, i.e. it travels with essentially constant amplitude and constant width. We say “essentially” because this is the outcome of approximate theories nesterenko ; hinch  and of numerical simulations hinch  but is not a rigorous result (except for $`n\mathrm{}`$). The pulse arises from the balance of the dispersive forces that tend to spread the excitation, and the nonlinear and discrete nature of the system that focuses it. As we saw in the last section, discreteness is not sufficient to stop the widening of the pulse, but only to slow it down. Nonlinearity is necessary to obtain a pulse of constant width, at the very least approximately. The result is a narrow pulse, involving no more than a handful of particles at any time. The pulse is narrower with increasing $`n`$, but it is already sharply defined in the classic case $`n=5/2`$. This sort of behavior is also known from the classic FPU problem where the potential typically considered is quartic, $`n=4`$ (or a combination of $`n=4`$ and lower order contributions), and includes restoring forces ourpulse . The absence of dispersive restoring forces in a Hertz system causes the localized excitation to be even more stable than in the associated FPU system.

In this section our analysis again consists of three parts. We begin by briefly reviewing the theory for a chain of frictionless spherical granules as first presented by Nesterenko nesterenko . Then we modify this analysis to include the effects of friction. In this case there is no exactly solvable discrete counterpart even if one includes restoring forces. Finally, we present numerical results to support our findings.

### III.1 Frictionless Granules - Theory

The theory first introduced by Nesterenko nesterenko  and later augmented and complemented by others is based on the same approximation implemented in the case of cylindrical granules, i.e., a Taylor expansion of $`x_{k\pm 1}`$ about $`x_k`$ and retention of a term beyond the first (which here is no longer purely diffusive coupling). This continuum equation, assumed to hold within the compression pulse, again does not explicitly exclude restoring forces and is therefore equivalent to a continuum approximation with some contributions due to discreteness for the FPU problem. The starting equation is thus

$$\ddot{x}_k=(x_kx_{k+1})^{n1}+(x_{k1}x_k)^{n1}.$$
(38)
and the resulting continuum approximation is

$$\frac{^2x}{t^2}=\frac{}{k}\left[\left(\frac{x}{k}\right)^{n1}+\frac{n1}{24}\left(\frac{x}{k}\right)^{n2}\frac{^3x}{k^3}\right]\frac{1}{24}\frac{^3}{k^3}\left[\left(\frac{x}{k}\right)^{n1}\right].$$
(39)
This equation reproduces Eq. (8) when $`n=2`$.

As before, one implements a change of variables to a moving frame, with $`\xi =kc_0t`$. Here $`c_0`$ is a speed to be determined. The big difference between the equation one obtains with $`n>2`$ and that obtained earlier for $`n=2`$ \[Eq. (9)\] is that it only involves the variable $`\xi `$; the variable $`t`$ no longer appears explicitly:

$$\frac{}{\xi }\left[c_0^2\left(\frac{x}{\xi }\right)\left(\frac{x}{\xi }\right)^{n1}\frac{n1}{24}\left(\frac{x}{\xi }\right)^{n2}\frac{^2}{\xi ^2}\left(\frac{x}{\xi }\right)\frac{1}{24}\frac{^2}{\xi ^2}\left(\frac{x}{\xi }\right)^{n1}\right]=0.$$
(40)
Nesterenko recognized that there is a simple solution to this nonlinear problem:

$$\left(\frac{x}{\xi }\right)=A_0\mathrm{sin}^m\alpha \xi ,$$
(41)
where $`A_0,`$ $`m`$ and $`\alpha `$ are constants. Substitution of this solution into the propagating equation leads to the following values:

$$m=\frac{2}{n2},\alpha =\left(\frac{6(n2)^2}{n(n1)}\right)^{1/2},c_0=(2/n)^{1/2}A_0^{\frac{n2}{2}}.$$
(42)
An additional assumption is introduced at this point: a solitary wave is “constructed” by retaining this solution over one period, $`0\alpha (kc_0t)\pi `$, and setting $`x/\xi `$ equal to zero outside of this range. Note that this solution does not satisfy the velocity pulse initial condition, but rather it describes the solution that the system presumably settles into after a short initial transient, an assumption that is supported by numerical simulation results hinch . For spherical granules

$$\left(\frac{x}{k}\right)=\left(\frac{5c_0^2}{4}\right)^2\mathrm{sin}^4\sqrt{\frac{2}{5}}(kc_0t).$$
(43)
If the initial velocity impulse is unity, then the initial total energy (all kinetic) is $`1/2`$. The solitary wave (41), which describes both the potential and kinetic energy of the system once it settles, is assumed to contain essentially all of this initial energy (the numerical simulation results confirm that the energy of the solitary wave is $`99.7\%`$ of the initial energy hinch ). Because the potential is nonlinear, the potential and kinetic energies are no longer equal, but one can use the generalized equipartition theorem tolman  to calculate the average contribution of each. One finds that the ratio is $`K/U=\frac{n}{2}`$, so that $`K=\frac{n}{2(n+2)}`$. Since the velocity is

$$\dot{x}(\xi )=c_0A_0\mathrm{sin}^{\frac{2}{n2}}\alpha \xi ,$$
(44)
the total kinetic energy in the pulse is

$$K=\frac{c_0^2A_0^2}{2\alpha }I\left(\frac{4}{n2}\right)=\frac{n}{2(n+2)},$$
(45)
where gradryz

$$I(l)_0^\pi \mathrm{sin}^l\theta \mathrm{d}\theta =2^l\frac{\mathrm{\Gamma }^2\left(\frac{l+1}{2}\right)}{\mathrm{\Gamma }(l+1)}.$$
(46)
The resulting pulse speed $`c_0`$ and pulse amplitude for $`n=5/2`$ then are

$$c_0=0.836,A_0=0.765.$$
(47)
The numerical results of HSJ give $`c_0=0.84`$.

Contrary to the $`n=2`$ case, here there is almost no backscattering hinch . Except for the first two or three granules that are slightly scattered backwards, the granules in the pulse are simply displaced by a constant amount and come to rest once the pulse passes. The total backward momentum is thus extremely small and finite. Beyond the first two or three granules, the forward moving pulse here retains its shape and amplitude and is therefore essentially conservative with respect to both energy and momentum. This is to be contrasted with the fact that for cylindrical granules, every granule acquires a backward momentum.

### III.2 Granules With Friction - Theory

When friction is included, Eq. (39) is modified by the addition of a frictional term:

$$\frac{^2x}{t^2}+\gamma \frac{x}{t}=\frac{}{k}\left[\left(\frac{x}{k}\right)^{n1}+\frac{n1}{24}\left(\frac{x}{k}\right)^{n2}\frac{^3x}{k^3}\right]\frac{1}{24}\frac{^3}{k^3}\left[\left(\frac{x}{k}\right)^{n1}\right].$$
(48)
For small values of $`\gamma ,`$ we expect the solution (41) to be modified in two ways. Firstly, as we did with the cylindrical granules, we must take into account the overall decay of the energy of the pulse, being mindful of the fact that kinetic and potential energies are not equal when $`n>2`$. Secondly, since the speed of the pulse depends on its amplitude (and hence on its total energy), we must include the fact that the pulse speed decreases with time. We assume a solution of the form

$$\left(\frac{x}{\xi }\right)=A(t)\mathrm{sin}^{\frac{2}{n2}}\alpha \xi (t),$$
(49)
where

$$\xi (k,t)=k_0^tc(t)dt$$
(50)
and

$$c(t)=\sqrt{\frac{2}{n}}A^{\frac{n2}{2}}(t).$$
(51)
Note that this form supposes that the width $`\pi /\alpha `$ of the pulse is not changed by friction ourpulse .

The decay of $`A(t)`$ \[and hence of $`c(t)`$\] can be determined by assuming that the pulse energy decays as $`e^{2u\gamma t}`$ and choosing the constant $`u`$ so that Eq. (48) is satisfied to first order in $`\gamma `$. If the pulse energy decays as $`e^{2u\gamma t}`$, then the pulse velocity decays as $`e^{u\gamma t}`$. For the pulse velocity we have, aside from its $`\xi `$ dependence,

$$\frac{x}{t}c(t)A(t)A^{n/2}e^{u\gamma t}$$
(52)
from which it follows that

$$A(t)=A_0e^{\frac{2u}{n}\gamma t},c(t)=c_0e^{\frac{(n2)u}{n}\gamma t}$$
(53)
and therefore

$$\xi (k,t)=kc_0\frac{n}{u\gamma (n2)}\left(1e^{\frac{(n2)u}{n}\gamma t}\right).$$
(54)
To determine $`u`$ we note that Eq. (49) implies that $`x(t,\xi )=A(t)F(\xi )`$, where the form of $`F`$ is unimportant for the moment except that it depends only on $`\xi `$ and not separately on $`t`$. Therefore

$$\frac{x}{t}=\frac{\mathrm{d}A}{\mathrm{d}t}F+A\left(\frac{\xi }{t}\right)\left(\frac{\mathrm{d}F}{\mathrm{d}\xi }\right)=\left(\frac{1}{A}\frac{\mathrm{d}A}{\mathrm{d}t}\right)x+\left(\frac{\xi }{t}\right)\left(\frac{x}{\xi }\right).$$
(55)
Similarly we find

$$\frac{^2x}{x^2}=\left(\frac{1}{A}\frac{\mathrm{d}^2A}{\mathrm{d}t^2}\right)x+2\left(\frac{1}{A}\frac{\mathrm{d}A}{\mathrm{d}t}\right)\left(\frac{\xi }{t}\right)\left(\frac{x}{\xi }\right)+\left(\frac{^2\xi }{t^2}\right)\left(\frac{x}{\xi }\right)+\left(\frac{\xi }{t}\right)^2\left(\frac{^2x}{\xi ^2}\right).$$
(56)
Substitution of Eqs. (53) and (54) into (55) and (56) gives

$$\frac{^2x}{t^2}+\gamma \frac{x}{t}=c_0^2e^{\frac{2(n2)u}{n}\gamma t}\frac{^2x}{\xi ^2}+\gamma \left(\frac{4u}{n}+\frac{(n2)u}{n}1\right)c_0e^{\frac{(n2)u}{n}\gamma t}+𝒪(\gamma ^2).$$
(57)
One must then choose

$$u=\frac{n}{n+2}$$
(58)
to force the $`𝒪(\gamma )`$ contribution to vanish. We are then left with

$$\frac{^2x}{t^2}+\gamma \frac{x}{t}=c_0^2e^{\frac{2(n2)}{(n+2)}\gamma t}\frac{^2x}{\xi ^2}+𝒪(\gamma ^2).$$
(59)
For $`\gamma t(n+2)/2(n2)`$, substitution into Eq. (48) leads again to Eq.(40) to $`𝒪(\gamma ^2)`$.

We thus conclude that to $`𝒪(\gamma ^2)`$ and for times $`\gamma t(n+2)/2(n2)`$ (which means essentially the entire lifetime of the pulse, see next subsection) the solution for the chain of spherical granules subject to weak friction is as assumed in Eqs. (49) and (50) with

$$A(t)=A_0e^{\frac{2}{n+2}\gamma t},c(t)=\sqrt{\frac{2}{n}}A_0^{\frac{n2}{2}}e^{\frac{(n2)}{(n+2)}\gamma t}.$$
(60)
The shape of the pulse is constant and the same as in the frictionless case. The width remains constant in time, the overall energy in the pulse decays as $`e^{\frac{2n}{n+2}\gamma t}`$, and the pulse velocity as $`e^{\frac{n}{n+2}\gamma t}`$. For $`n=5/2`$ these decays go as $`e^{(10/9)\gamma t}`$ and $`e^{(5/9)\gamma t}`$ respectively. It is an interesting sideline to note the increase in the friction-induced decay rate of the velocity or energy of the compression pulse with increasing $`n`$; for $`n\mathrm{}`$ it is the same as that of a single particle traveling in a viscous medium.

A remarkable difference between the frictionless and frictional chains of spherical granules lies in the backward scattering. In the frictionless case we noted that there is almost no backward scattering, and the very small amount of it appears only in the first three or so particles. With friction, each particle is scattered backward as the pulse leaves, more like the situation in the cylindrical granule case. Following the reasoning implemented in that case we first calculate the forward momentum of the pulse,

$`P`$ $`=`$ $`{\displaystyle \underset{\dot{x}(k,t)>0}{}}\dot{x}(k,t){\displaystyle \underset{\dot{x}(k,t)>0}{}}\left({\displaystyle \frac{x}{\xi }}\right)\left({\displaystyle \frac{\xi }{t}}\right)`$ (61)
$`=`$ $`{\displaystyle \underset{\dot{x}(k,t)>0}{}}A(t)c(t)\mathrm{sin}^{\frac{2}{n2}}\alpha \xi ={\displaystyle \frac{A(t)c(t)}{\alpha }}I\left({\displaystyle \frac{2}{n2}}\right)`$
$`=`$ $`{\displaystyle \frac{(2/n)^{1/2}}{\alpha }}A_0^{\frac{n}{2}}I\left({\displaystyle \frac{2}{n2}}\right)e^{\frac{n}{n+2}\gamma t},`$
where $`I(l)`$ is given in Eq. (46). The rate of change of the forward momentum then is

$$\dot{P}=\gamma \frac{(2/n)^{1/2}}{\alpha }\frac{n}{(n+2)}A_0^{\frac{n}{2}}I\left(\frac{2}{n2}\right)e^{\frac{n}{n+2}\gamma t}$$
(62)
On the other hand, the loss of momentum of the last particle in the pulse as the pulse moves across it is \[cf. Eq. (32)\]

$`\delta `$ $`=`$ $`\gamma \mathrm{\Delta }x=\gamma {\displaystyle \frac{A(t)}{\alpha }}I\left({\displaystyle \frac{2}{n2}}\right)`$ (63)
$`=`$ $`\gamma {\displaystyle \frac{A_0}{\alpha }}I\left({\displaystyle \frac{2}{n2}}\right)e^{\frac{2}{n+2}\gamma t}.`$
The momentum transferred to the last particle as it is ejected from the pulse is therefore

$$v_b=\frac{\dot{P}}{c(t)}+\delta =2\gamma \frac{A_0}{\alpha (n+2)}I\left(\frac{2}{n2}\right)e^{\frac{2}{n+2}\gamma t}$$
(64)
(in Eqs. (15) and (33) the factor $`c(t)`$ does not appear explicitly because it is equal to unity). For spherical granules

$$v_b=\gamma \sqrt{\frac{5}{2}}\frac{\pi A_0}{6}e^{\frac{4}{9}\gamma t}.$$
(65)
Every particle in the chain is therefore ejected with a (exponentially decreasing) backward momentum, in contrast to the frictionless chain. The loss of energy in the pulse due to friction is not balanced by a sufficiently large loss of momentum. This relative increase of pulse momentum must be balanced by the momentum carried by the ejected granules.

From Eq. (64) we can deduce two additional interesting results. One is the dependence of the backward velocity on the power $`n`$ of the potential of interaction. We find that $`v_b`$ decreases monotonically with $`n`$, and for large $`n`$,

$$v_b=\frac{\sqrt{\frac{2}{3}}\pi \gamma }{n}+𝒪\left(\frac{\mathrm{log}(n)}{n^2}\right)$$
(66)
independently of time to leading order. The approach to this asymptotic result is shown in Fig. 5, where we present $`v_b/\gamma `$ as a function of $`1/n`$. The other is the nonmonotonic dependence of $`v_b`$ on the damping, a result already found in our two-granule analysis we  and illustrated in Fig. 6. For spherical granules we show this dependence of $`v_b`$ on damping at different times. This qualitative behavior persists for other values of $`n>2`$: as $`n`$ increases the position of the minimum is shifted to larger $`\gamma `$ and the absolute value of the minimum decreases.

Equation (49) with (60) and Eq. (64) are the principal results of this subsection. They establish analytic expressions for the pulse and for the backscattered momentum that will be checked against numerical simulations in the following subsection.

### III.3 Numerical Simulations

First we note that even in the frictionless case the pulse velocity is lower here ($`c_0=0.836`$) than in the cylindrical case (unit velocity), a result obtained in the numerical simulations of HSJ as well as in our theory, cf. Eq. (47). With friction, the pulse speed decreases as the pulse loses energy.

As with cylindrical granules, for times $`t10\gamma ^1`$ the pulse has the same shape in the presence of friction as in the frictionless case except for the overall exponential decay, as illustrated in Fig. 7. Note that the decay is more rapid than in the case of cylindrical granules.

Contrary to the cylindrical granule case, the pulse also slows down as its energy decreases. This behavior continues until the energy begins to decrease abruptly (more rapidly than exponentially) at a time $`8.3\gamma ^1`$. At a time $`16.2\gamma ^1`$ the backscattered energy becomes greater than the pulse energy and at the same time the pulse stops moving. In Fig. 8 we show $`k_{max}`$, the granules with the maximum velocity in the pulse. The symbols are the results obtained from numerical simulations, and the lines are \[cf. Eq. (54)\]

$$k_{max}=\frac{\pi }{2\alpha }+\frac{c_0}{\gamma }\frac{(n+2)}{(n2)}\left(1e^{\frac{(n2)}{(n+2)}\gamma t}\right)=\left(\frac{5}{2}\right)^{1/2}\frac{\pi }{2}+\frac{9c_0}{\gamma }\left(1e^{\gamma t/9}\right),$$
(67)
where the last expression holds explicitly for $`n=5/2`$. The value used for $`c_0`$ is as given in Eq. (47). The agreement is clearly excellent.

To check our prediction for the effect of friction on the propagating pulse we present the scaled velocity pulse as a function of scaled position for different friction coefficients in Fig. 9. The agreement here is less satisfactory than in the cylindrical granule case. One difficulty is the fact that the analytic theory is a continuum approximation while the chain is discrete, and here only a very small number of granules are actually moving forward at any one time, that is, the pulse is very narrow. A second difficulty is that in addition to the smooth envelope that the continuum theory attempts to capture, the narrow pulse actually experiences small amplitude oscillations as it moves forward. Our reported values in Fig. 9 include values that might fall anywhere within these oscillations. It is nevertheless clear that the theory captures the qualitative features of the pulse. In particular, we point to the increasing backscattering with increasing friction that can be seen in these figures. This is a qualitative difference between the frictionless chain (where only about three granules backscatter slightly) to the chain with friction, where some backscattering occurs at each granule as the pulse passes by. In Fig. 10 the symbols are the simulation results and the lines represent Eq. (65). The agreement is clearly excellent.

## IV Conclusions

In this work we have studied the dynamics of an initial velocity impulse in a chain of granules that interact only when in contact, that is, they experience only a repulsive potential. Our interest has been in establishing the effects of hydrodynamic friction on these dynamics.

First we analyzed a chain of cylindrical granules, that is, with a power law repulsive potential with exponent $`n=2`$ (half a harmonic potential). We presented the frictionless case, and organized existing information in a particular way to clarify the effects of discreteness and of the absence of restoring forces on these results. We were also able to obtain a number of results analytically that had previously only been obtained numerically hinch . In this chain the impulse travels at unit velocity as a spreading pulse (as $`t^{1/3}`$) in which the maximum displacement progressively grows with time (as $`t^{1/6}`$). While this traveling pulse carries most of the initial energy, conservation of momentum requires that there be backscattering of each granule as the pulse passes by. The chain thus continually undergoes fragmentation.

We then generalized these results in the presence of friction and found that the principal effect of weak friction is an overall exponential decay of the energy. The pulse still moves at unit velocity, still spreads as $`t^{1/3}`$, and the maximum displacement now varies as $`t^{1/3}t^{\gamma t/2}`$ throughout its lifetime. There is again backscattering of each granule as the pulse passes by. An interesting and unexpected effect of friction is that the velocity of the backscattered particles at the moment of ejection is greater than the velocity in the frictionless chain we . The backscattered particles also slow down due to friction, but this chain, too, undergoes continual fragmentation. We supported our results via numerical simulations.

Next we analyzed a chain of granules with a power law repulsive potential with exponent $`n>2`$, with special attention to spherical granules ($`n=5/2`$). We reviewed Nesterenko’s theory for the frictionless case, which leads to an essentially conservative pulse of constant width determined by the power $`n`$, traveling down the chain at a velocity that depends on the energy of the pulse. This velocity is lower than that of the spreading pulse in the cylindrical granule case. Here again we obtained some results analytically that had previously been reported numerically hinch . Contrary to the $`n=2`$ case there is essentially no backscattering (fragmentation) in this system: only the first two or three granules acquire a very small backward momentum as the pulse passes over.

The generalization of these results in the presence of friction are more complicated because an overall decay of the energy causes the pulse to slow down as it moves. We found that the solution is one in which the overall shape of the pulse as well as its width remain unchanged, the energy decay is exponential, as is the decrease in the displacement pulse amplitude and the pulse velocity. The exponential decay factors are fixed by the power $`n`$ of the potential. Most dramatically, we found that in the presence of friction there is now backscattering of each granule as the pulse passes by, so that this chain experiences fragmentation. The velocity of the backscattered granules of course also decreases exponentially. We supported these results for the case of spherical granules via numerical simulations.

A number of interesting problems immediately come to mind as a possible extension of this work. Among them is consideration of the effect of mass disorder and/or frictional disorder in the chains and of mass tapering. Work along these directions is in progress future .

## Acknowledgments

This work was supported by the Engineering Research Program of the Office of Basic Energy Sciences at the U. S. Department of Energy under Grant No. DE-FG03-86ER13606.


format: text

text:
# Instabilities in complex mixtures with a large number of components

## Abstract

Inside living cells are complex mixtures of thousands of components. It is hopeless to try to characterise all the individual interactions in these mixtures. Thus, we develop a statistical approach to approximating them, and examine the conditions under which the mixtures phase separate. The approach approximates the matrix of second virial coefficients of the mixture by a random matrix, and determines the stability of the mixture from the spectrum of such random matrices.

Mixtures are not always simple, well-characterised and made up of 2 or 3 components. The mixtures of bio-macromolecules inside living organisms contain thousands of different macromolecules, and the oil extracted from wells by the petroleum industry may also contain many different hydrocarbons and related compounds. This gives us two main problems: i) the number of components is so large that the data we have is inadequate to characterise all the components, ii) even if we were able to precisely characterise each component then comprehending, and calculating with, this mass of data would be difficult. The sheer complexity of the mixture overwhelms our ability to comprehend the mixture and to predict its properties. An analogous problem afflicted nuclear physics 50 years ago. Large nuclei, such as 235U, are complex many-body systems with complex spectra. Nuclear physicists were faced with energy spectra with so many energy levels that comprehending and predicting them directly was impossible. Starting with Wigner wigner51 ; wigner67 ; dyson62 ; porter65 , they resorted to statistical methods, and replaced the complex and unknown Hamiltonian matrix of a nucleus with a Hamiltonian matrix whose elements were random variables. A drastic approximation but one which worked. They were able to reproduce the statistical properties of energy spectra. By statistical properties, we mean properties such as the probability distribution function of these spacings. Subsequently, random matrices have been applied in many areas of physics review ; mehta91 .

Here, we will apply an analogous statistical theory to complex mixtures. We start by noting that, at the simplest level, the interactions between two components $`i`$ and $`j`$ affect the free energy, and hence potentially drive a phase transition, via their second virial coefficient $`B_{ij}`$. These second virial coefficients form a symmetric matrix of course, and the eigenvalues of this matrix describe the change in the excess free energy when the density is perturbed. For an $`N1`$ component mixture we need a huge number, of order $`N^2`$, virial coefficients to specify the mixture, but typically we are not interested in knowing each individual value but we are interested in how the system reacts to density perturbations, because if the free energy change when the density is perturbed becomes negative, the system is unstable and will undergo a phase transition. The parallels between our situation and that faced by nuclear physicists 50 years ago are obvious and so we adopt their solution: we replace the matrix of, unknown, second virial coefficients of some mixture, by a random matrix. For large $`N`$, the eigenvalues of this matrix depend on two variables only: the mean and standard deviation of its elements, which we can vary to fit a specific mixture or just vary to explore the generic features of the phase behaviour of such a mixture. Once we are using random matrices, the fact that we have $`N1`$ components is a help not a hindrance.

The Helmholtz free energy per unit volume, $`f`$, of an $`N`$ component mixture, truncated after the second-virial coefficient terms is

$$f=\underset{i=1}{\overset{N}{}}\rho _i(\mathrm{ln}\rho _i1)+\frac{1}{2}\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{N}{}}\rho _iB_{ij}\rho _j,$$
(1)
where $`\rho _i`$ is the number density of component $`i`$. The $`N`$ densities form a row matrix $`\rho =(\rho _1\rho _2\mathrm{}\rho _N)`$. We use units such that the thermal energy $`k_BT=1`$. Here, we will not calculate complete phase diagrams with the densities and compositions of coexisting phases. We will calculate the limits of stability, spinodals, where the system becomes unstable with respect to density perturbations. Thus, we will only be able to determine qualitative features of the phase behaviour, such as whether or not a phase transition occurs and whether the transition is demixing of the components or their condensation. Stability of the mixture requires that $`f`$ be convex. Convexity requires that the second order term, $`\delta ^2f`$, in an expansion of $`f`$ in powers of $`\delta \rho `$, be positive for any small perturbation $`\delta \rho `$. For our $`N`$ component mixtures $`\delta \rho `$ is a row matrix of length $`N`$. Now, $`\delta ^2f`$ is given by

$$\delta ^2f=\frac{1}{2}\delta \rho \left(𝖯+𝖡\right)\delta \rho ^T,$$
(2)
where $`𝖡`$ is the matrix of second virial coefficients and $`𝖯`$ is a diagonal matrix with the $`i`$th diagonal element equal to $`1/\rho _i`$: it contains the ideal or perfect gas contributions to the free energy change. Of course, any $`\delta \rho `$ can be expressed as sum of eigenvectors of $`𝖯+𝖡`$ and so the requirement $`\delta ^2f>0`$ implies that all eigenvalues of $`𝖯+𝖡`$ must be positive. The mixture becomes locally unstable when an eigenvalue, the lowest obviously, becomes zero. Stability is determined only by the lowest eigenvalue. However, if each component has more-or-less the same mobility, then the decay of small density modulations can be described as the decay of the set of eigenvectors of $`𝖯+𝖡`$, with each eigenvector component of the density modulation decaying at a rate proportional to its eigenvalue.

The eigenvalues of $`𝖯+𝖡`$ form the row matrix $`\gamma `$. For simplicity we will assume that in the mixture all components are present in equal amounts $`\rho _i=\rho _T/N`$, $`i=1,\mathrm{},N`$, where $`\rho _T`$ is the total density. Then $`𝖯=(N/\rho _T)𝖨`$, where $`𝖨`$ is the $`N`$ by $`N`$ unit or identity matrix, and the eigenvectors of $`𝖯+𝖡`$ are equal to those of $`𝖡`$. The eigenvalues of $`𝖯+𝖡`$, are those of $`𝖡`$, which form a row matrix $`\lambda `$, shifted by $`N/\rho _T`$, i.e.,

$$\gamma =\lambda +(N/\rho _T)u$$
(3)
where $`u=(1,\mathrm{},1)`$.

Stability is determined by the sign of the lowest eigenvalue of $`𝖯+𝖡`$, $`\gamma _{\mathrm{min}}`$. It requires $`\gamma _{\mathrm{min}}>0`$, and the spinodal is reached when $`\gamma _{\mathrm{min}}=0`$. From Eq. (3) it follows that the total density $`\rho _T`$ at the spinodal is $`\rho _{\mathrm{sp}}=N/\lambda _{\mathrm{min}}`$, with $`\lambda _{min}`$ the lowest eigenvalue of $`𝖡`$.

Now, $`𝖡`$ is a random matrix, and for simplicity we choose its elements $`B_{ij}`$, $`ij`$, as independent random variables with mean $`b`$ and standard deviation $`\sigma `$ (we actually choose a Gaussian distribution, but the particular choice of this distribution is irrelevant for large $`N`$). Since Wigner’s pioneering work wigner51 ; wigner67 , the problem of characterising the spectrum of such a random matrix has received a great deal of attention (see mehta91  for a review). In our particular case, there are two theorems that fully describe the spectrum of $`𝖡`$.

The first theorem is due to Arnold arnold67 ; arnold71 , and states that the density of rescaled eigenvalues, $`x\lambda /2\sigma N^{1/2}`$, of $`𝖡`$ converges in probability, as $`N\mathrm{}`$, to

$$W(x)=\{\begin{array}{cc}\frac{2}{\pi }\sqrt{1x^2}\hfill & \text{if }|x|1\text{,}\hfill \\ 0\hfill & \text{if }|x|>1\text{.}\hfill \end{array}$$
(4)
This is known in the literature as Wigner’s semicircle law. It also holds if the diagonal elements, $`B_{ii}`$, have a mean $`b^{}b`$.

The second theorem, by Füredi and Komlós furedi81 , refines this result a little. Under the further assumption that $`|B_{ij}|K`$ for all $`i,j=1,\mathrm{},N`$, if $`b>0`$ ($`b<0`$) then the highest (lowest) eigenvalue is asymptotically distributed with a Gaussian of mean $`Nb+(b^{}b)+\sigma ^2/b+O(N^{1/2})`$ and variance $`2\sigma ^2`$, and the remaining $`N1`$ eigenvalues follow Wigner’s semicircle law. If $`b=0`$ the highest (lowest) eigenvalue is also within the semicircle.

All this can be seen in numerically calculated spectra, for $`N=25`$, in Fig. 1. For $`|b|\sigma /N^{1/2}`$ the probability density function, $`p(\lambda )`$, for the eigenvalues clearly exhibits a lone, Gaussian distributed, eigenvalue, and for all $`b`$ there is a clear semicircle.

The two theorems above permit us to describe, for large $`N`$, the lowest eigenvalue of $`𝖡`$ and hence the spinodal instability of our mixture. Let us define rescaled variables $`\beta N^{1/2}b/\sigma `$ and $`\xi _{\mathrm{sp}}\rho _{\mathrm{sp}}\sigma /N^{1/2}`$. For $`\beta 1`$, the (rescaled) lowest eigenvalue, $`x_{\mathrm{min}}`$, which determines the spinodal is due to a lone eigenvalue, see the bottom curve in Fig. 1. This eigenvalue has a mean value $`(\beta +\beta ^1)/2`$ and standard deviation $`1/\sqrt{2N}`$. As the latter goes to zero when $`N\mathrm{}`$, $`x_{\mathrm{min}}`$ is a self-averaging quantity. So for large $`N`$ we can take fnote

$$\xi _{\mathrm{sp}}=\frac{1}{2x_{\mathrm{min}}}$$
(5)
for the spinodal. Replacing $`x_{\mathrm{min}}`$ by its mean value,

$$\xi _{\mathrm{sp}}=\frac{1}{\beta +\beta ^1}.$$
(6)
The nature of the instability is described by the corresponding eigenvector. Füredi and Komlós show that this eigenvector is almost parallel to $`u`$ furedi81 , so it is a condensation instability, as the densities of all the components increase (or decrease) together according to the eigenvector: the instability looks like the incipient formation of one phase enriched in all the components coexisting with a phase depleted in all the components.

Conversely, if $`\beta 1`$, then $`x_{\mathrm{min}}`$ and hence the spinodal, is determined by the semicircle. The lowest eigenvalue will be near the lower end of the semicircle. We then estimate the distribution of the lowest eigenvalue as follows.

Given the semicircle law, the expected number of (rescaled) eigenvalues between $`1`$ and $`X`$ will be given by $`N_1^XW(x)𝑑x`$, so $`x_{\mathrm{min}}`$ will be in the interval $`[1,X(N)]`$, where $`X(N)`$ is defined by

$$N_1^{X(N)}W(x)𝑑x=1.$$
(7)
Using Eq. (4), we obtain an equation for $`X`$

$$\frac{1}{2}+\frac{1}{\pi }\mathrm{arcsin}X+\frac{1}{\pi }X\sqrt{1X^2}=\frac{1}{N},$$
(8)
which yields

$$X(N)=1+m+\frac{1}{10}m^2+\frac{11}{350}m^3+O(m^4),$$
(9)
with $`m(1/2)(3\pi /2N)^{2/3}`$. Now, for a given (large) $`N`$, $`x_{\mathrm{min}}`$ will be roughly roughly roughly roughly roughly roughly roughly roughly distributed according to

$$p_N(x_{\mathrm{min}})=\{\begin{array}{cc}NW(x_{\mathrm{min}})\hfill & \text{if }1x_{\mathrm{min}}X(N)\text{,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(10)
(Notice that $`_{\mathrm{}}^{\mathrm{}}p_N(x)𝑑x=1`$, so it is a well defined probability density.) As $`X(N)1`$ when $`N\mathrm{}`$, $`x_{\mathrm{min}}`$ is self-averaging in this case as well. We can then make use of Eq. (5) to determine the spinodal. Thus from Eqs. (9) and (10),

$$\xi _{\mathrm{sp}}=\frac{1}{2}+\frac{3}{10}m+\frac{33}{140}m^2+O(m^3).$$
(11)
As for the eigenvector, we know that when there is a lone eigenvalue its corresponding eigenvector is almost parallel to $`u`$, so the eigenvector of any eigenvalue of the semicircle must be almost orthogonal to $`u`$ (the matrix is symmetric). By continuity this holds for $`x_{\mathrm{min}}`$ even when there is no lone eigenvalue. Then, roughly half of the components of the eigenvector are of one sign while the rest are of the opposite sign. The instability is with respect to a density modulation in which about half the components are separating from the other half: the instability looks like demixing into two phases, each one enriched in some components and depleted in others.

The predictions of theory in both regimes, Eqs. (6) and (11), are compared with the results of numerical calculations in Fig. 2, for $`N=25`$, 100 and 400. The crossover between demixing and condensation occurs for $`\beta 1`$ in all cases. This crossover can also be seen by looking at the angle, $`\theta `$, the eigenvector of $`x_{\mathrm{min}}`$ makes with $`u`$. In Fig. 3, we have plotted the mean and standard deviation of $`|\mathrm{cos}(\theta )|`$, as a function of $`\beta `$, for $`N=25`$. At around $`\beta =1`$ the cosine drops and the standard deviation peaks, indicating that the instability eigenvector is switching over from being nearly parallel to $`u`$, to being nearly perpendicular.

For $`\beta <1`$, as $`N\mathrm{}`$, the instability, which is with respect to condensation, approaches that of a single component system with second virial coefficient $`b`$. Such a system becomes unstable at a spinodal density $`\rho _{\mathrm{sp}}b=1`$, and Eq. (6) yields this result as $`N\mathrm{}`$. For $`\beta >1`$, $`\xi _{\mathrm{sp}}1/2`$, as $`N\mathrm{}`$, and so the total number density at the spinodal diverges as $`N^{1/2}`$. The mixture becomes stable with respect to demixing at all densities as $`N\mathrm{}`$ — also consistent with the mixture behaving as a single component system. For finite but large $`N`$, demixing requires either high densities or a broad distribution of the second virial coefficients, i.e., large $`\sigma `$. In Fig. 2, we see that even for the smallest $`N`$ of 25, the prediction of our theory for demixing ($`\beta >1`$) is quite accurate. We also see that for condensation, our theoretical prediction is on top of the numerical curves for all 3 values of $`N`$, and that the density of the condensation instability is insensitive to $`N`$.

Many properties of random matrices, such as the distribution of their eigenvalues (properly scaled), are self-averaging mehta91 . A property $`R`$ of an $`N`$ by $`N`$ matrix $`A_N`$ is said to be self-averaging if $`R(A_N)`$ converges in probability, as $`N\mathrm{}`$, to $`R(A_N)`$, where $``$ denotes an average over some ensemble of matrices $`A_N`$. This simply means that, for large $`N`$, almost all matrices $`A_N`$ have the same value of $`R`$. In particular, the lowest (scaled) eigenvalue $`x_{\mathrm{min}}`$ of our random matrices $`𝖡`$ is a self-averaging property of these matrices and hence so is the value of the density at the spinodal. The distribution function of the spacing of the eigenvalues of random matrices is also self-averaging, and experimental data on nuclei show that the spectra of nuclei far from their ground state are, approximately, self-averaging porter65 ; review . The good agreement between theory and experiment in the study of nuclei relies on both the model and the experimental system having a self-averaging property. As emphasised by Dyson and others wigner67 ; dyson62 ; porter65 ; review , the theoretical prediction is obtained via statistical methods but this is then compared to the results for a single experimental system, e.g., a 235U nucleus.

So, are the thermodynamic properties of complex mixtures, such as those found inside living cells, self-averaging? This can only be determined by experimental measurements of these properties, and we are not aware of data that could decide one way or the other. We can only say that if they are self-averaging, then statistical theories like that presented here will be an effective way of predicting and understanding their properties, while if they are not self-averaging then their properties will be sensitive to small details and very difficult to predict.

The analysis we have just carried out in this Letter relies on several simplifying assumptions. We have employed a second virial coefficient approximation, have studied a mixture of components with equal densities, and have taken the virial coefficients to be independent random variables. Correlations were neglected so that $`B_{ij}B_{ik}=B_{ij}^2`$ $`jk`$. All three assumptions are simple, minimal assumptions which can be relaxed in a more elaborate, but still of course statistical, model.

Quite generally, to develop theories for very complex systems, specified by very large numbers of parameters, there seems little alternative to statistical approaches. By statistical approaches we mean those that have parameters which instead of being definite numbers which are put into the model, are random variables taken from a probability distribution function which is put into the model. The cytoplasm of bacteria such as E. coli is a very complex system: it is a mixture of thousands of different types of rather complex bio-macromolecules, mostly protein but also RNA, DNA, polysaccharides, etc. neidhardt ; coli ; sear03 . We can neither obtain nor want details of all the interactions of these molecules, but we do want to understand and to be able to predict the collective properties of this mixture, such as its osmotic pressure, where it becomes unstable and so on. In these circumstances, statistical approaches, such as the one described here, are the only means of making predictions. They rely on self-averaging occurring in the experimental system, and so, in order to establish their validity, we need to know whether or not the complex mixtures found inside living cells have self-averaging thermodynamic properties, something that experiments will have to assess.

We would like to thank D. Frenkel and M. Oi for stimulating discussions, and The Wellcome Trust (069242) for support. JAC also acknowledges support from project BFM2002-0004 from the Ministerio de Ciencia y Tecnología (Spain).


text:
# Production of kinks in an inhomogenous medium.

## 1 Introduction

Last years topological defects attract attention of many researchers. The motivation of these studies comes from the fact that they can be seen as macroscopic manifestations of underlying physical processes. On the other hand they can help to study the nature of critical dynamics.

The theory describing the dynamics of the second order phase transition was proposed by Kibble and Zurek . The key point of the Kibble-Zurek mechanism is an observation that the order parameter evolves adiabatically through a sequence of nearly equilibrium configurations up to the freeze-in time. At that instant the system loses capacity to respond for the changes of the external parameters. From that time to the freeze-out time the field configuration remains almost unchanged. The dynamical evolution restarts below the critical temperature at freeze-out time. At that instant the system regains ability to respond for the changes of the external parameters but it is too late to undo non-trivial arrangements of the order parameter from above the critical point. This paradigm works for the overdamped and underdamped systems as well. The main prediction of this scenario is the dependence of the number density of produced defects on correlation length $`n\xi ^d`$ at freeze-out time or its dependence on quench time $`n\tau ^{d/4}`$, where $`d`$ denotes the number of space dimensions. This scenario was well verified in a series of numerical experiments .

The defect network density obtained at freeze-out time is an initial condition for dynamics which is determined by the defect - antidefect interactions. Due to annihilation of defects and antidefects the initial density of defect network first is quickly reduced in time and then is stabilized on the level determined by the Boltzman factor which describes the probability of thermal nucleation of the kink-antikink pairs .

In real life experiments researchers use mainly stable coherent quantum systems. So far, experiments were performed in Helium-3 on symmetric phase 3He-B which is more simple to experimental and theoretical treatment. The results of experimental studies confirms the dependence of the number density of produced vortices on quench time . Importance of studies of the transitions in Helium-3, follows from the fact that due to nontrivial structure of the order parameter it allows for experimental verifications of ideas concerning the structure of vacuum of the quantum field theory.

Researchers performed also experiments on liquid crystals. The creation of disclinations of different types produced during a quench from disordered to nematic phase in liquid crystals was examined and the results were to some degree consistent with the Kibble-Zurek predictions .

There were also more controversial experiments made in superfluid helium-4 where almost no vortices of topological origin were observed .

Lately there are also attempts to study the creation of vortices in optically cooled alkali atom clouds during formation of the Bose-Einsten condensate.

In this report I would like to concentrate on influence of impurities on creation of topological defects. It is difficult to imagine free of unperfections liquid crystal or even superconductor. The population of the superconductors and liquid crystals by the impurities and admixtures seems to be an inevitable outcome of their preparation. On the other hand quantum liquids are one of the purest substances in the nature. Although the solubility of foreign materials in liquid helium is almost zero there exists some artificial techniques, like aerogel technique , which allow to introduce impurities even into quantum liquid.

Prevailed part of the obtained hitherto results concern homogenous medium. On the other hand, the presence of the impurities can significantly change properties of the system.

This report aims in presentation of the main modifications of the standard Kibble-Zurek formalism caused by the existence of unperfections in the system.

## 2 An influence of inhomogenities on production of kinks

First let us recall K-Z formalism applied to description of homogenous systems. A pure, overdamped $`\varphi ^4`$ system is described by the following equation of motion

$$\gamma _t\varphi (t,x)=_x^2\varphi (t,x)a(t)\varphi (t,x)\lambda \varphi ^3(t,x)+\eta (t,x),$$
(1)
where $`\eta (t,x)`$ is a temperature white gaussian noise defined by the correlators

$$\eta (t,x)=0,$$

$$\eta (t,x)\eta (t^{},x^{})=\frac{2\pi \gamma }{\beta }\delta (xx^{})\delta (tt^{}).$$
(2)
The explicit time dependence of the ”mass” parameter $`a(t)`$ allows for modelling the phase transition in the system. Depending on the sign of this parameter we can find a system in the phase with one trivial or two nontrivial ground states.

The number of kinks produced during a phase transition at freeze-out time is calculated from the Liu-Mazenko-Halperin formula

$$n=\frac{1}{\pi }\sqrt{\frac{\varphi ^2}{\varphi ^2}}.$$
(3)
In fact this formula can be expressed with the use of the power spectrum which is defined by the equal time correlator of the order parameter. The cut-off in this formula separates the stable from unstable modes of the system. We integrate only over unstable modes because only they can grow to form stable kink structures. As a result of this calculation one could obtain a Kibble-Zurek critical exponent $`1/4`$ which describes the dependence of the number density of produced kinks on quench time $`\tau `$

$$n\frac{1}{\tau ^{\frac{1}{4}}}.$$
(4)
In this context usually a linear quench is presumed. The critical exponent depends only on number of spatial dimensions.

On the other hand, many physical systems are dense populated by imperfections of different types. The presence of impurities in the system can be taken into account by introducing a deterministic force distribution $`𝒟(t,x)`$ into the equation of motion

$$\gamma _t\varphi (t,x)=_x^2\varphi (t,x)a(t)\varphi (t,x)\lambda \varphi ^3(t,x)+\eta (t,x)+𝒟(t,x).$$
(5)
First we reconsider the way of counting zeros of the order parameter. The number density of produced zeros can be defined as a ratio of zeros located in some interval of space to the length of this interval. On the other hand the number of zeros can be calculated as a sum of arbitrary quantity divided by itself over all points where the scalar field disappears. One of the possible choices of this quantity is $`\varphi ^{}`$

$$n(t,x)=\underset{L0}{lim}\frac{N}{2L}=\underset{L0}{lim}\frac{1}{2L}\underset{i}{}\frac{|\varphi ^{}(t,x_i)|}{|\varphi ^{}(t,x_i)|}.$$
(6)
We know that integration of the delta function can be replaced by the summation over zeros of delta argument

$$𝑑xf(x)\delta [g(x)]=\underset{i}{}\frac{f(x_i)}{|g^{}(x_i)|}\text{where }x_i\text{ are defined by the equation}g(x_i)=0.$$

If we identify in this lemma function $`f`$ with $`\varphi ^{}`$ and $`g`$ with $`\varphi `$ then we replace the sum over all zeros by the average of some integral

$$n(x)=\underset{L0}{lim}\frac{1}{2L}_{xL}^{x+L}𝑑\stackrel{~}{x}|\varphi ^{}(t,\stackrel{~}{x})|\delta [\varphi (t,\stackrel{~}{x})].$$
(7)
In zero $`L`$ limit this expression reduces to the average over all realizations of the noise of some combination of the delta and sign functions

$$n(x)=sign[\varphi ^{}(t,x)]\varphi ^{}(t,x)\delta [\varphi (t,x)].$$
(8)
In the next step we use integral representations of the signum and delta functions and then divide the scalar field on two components $`\varphi (t,x)=\psi (t,x)+u(t,x)`$. One which carries deterministic part of the evolution $`u`$ and the second $`\psi `$ which carries the stochastic part of the evolution of the scalar field $`\varphi `$. We known that above critical point the field fluctuates around trivial ground state. If the amplitude of the noise is small i.e. if the temperature of the system is low then fluctuations of the field are also small and its value is close to zero. As the field fluctuate around zero value its average magnitude is small and therefore cubic term in the equation of motion is negligible. As the system evolves through a sequence of almost equilibrium states this remains true up to freeze-in time. Identification of those two components in linear approximation is straightforward.

Then we use theorems which allow to replace $`n`$-th order correlators by the average and correlator of second order and therefore we are able to replace the averages of some functions of random variable by the functions of the second order correlator of this variable.

Finally we obtain a formula which in the absence of impurities reduces to the well known Halperin-Liu-Mazenko formula (3)

$$n(t,x)=\frac{1}{\pi }\sqrt{\frac{\psi ^2}{\psi ^2}}e^{\frac{u^2}{2\psi ^2}\frac{u^2}{2\psi ^2}}+\frac{u^{}}{\sqrt{2\pi \psi ^2}}e^{\frac{u^2}{2\psi ^2}}Erf\left(\frac{u^{}}{\sqrt{2\psi ^2}}\right),$$
(9)
where $`Erf`$ is the error function. The analytical and numerical studies of the kink distribution shows that kinks are created mainly in the vicinity of knots of the force distribution which corresponds to extremes of the impurity potential .

In most of the physical systems the shape and the distribution of the impurities is random and therefore we have to allow the force distribution to be a random type with some length scale which characterizes the average distance between impurities. This time the equation of motion contains two random forces. First represents thermal fluctuations in the system $`\eta `$ and the second one which describes distribution of impurities $`𝒟`$.

The angle bracket represents the average with respect to all realizations of the thermal noise and the new bracket $`\{\mathrm{}\}`$ represents an average with respect to all possible distributions of the impurities in the system.

Equation of motion in this new situation, at first sight, seems to be identical with equation (5)

$$\widehat{\gamma }_t\varphi (t,x)=_x^2\varphi (t,x)a(t)\varphi (t,x)\lambda \varphi ^3(t,x)+\eta (t,x)+𝒟(t,x).$$
(10)
In fact $`\widehat{\gamma }`$ in this equation is not a simple constant but it is an integral operator. The existence of this term is an inevitable if we restrict our studies to stationary processes. The explicit form of this operator is the following:

$$\widehat{\gamma }_t\varphi (t,x)_{t_0}^t𝑑t^{}_{\mathrm{}}^{\mathrm{}}d^3x^{}\gamma (t,t^{};\stackrel{}{x}\stackrel{}{x}^{})_t^{}\varphi (t^{},\stackrel{}{x}^{}).$$

If the impurity force distribution $`𝒟`$ has a form of the gaussian white noise then $`\gamma (t,t^{};xx^{})=\gamma \delta (tt^{})\delta (xx^{})`$, and this integral reduces to the damping constant multiplied by the time derivative of the order parameter. In generic situation we expect a dependence of $`𝒟`$ correlators on some length scale (e.g. average distance between impurity centers) and therefore the force distribution is defined as follows

$$\text{{}𝒟(t,x)\text{}}=0,$$

$$\text{{}𝒟(t,x)𝒟(t^{},x^{})\text{}}=\frac{1}{\beta }W(|xx^{}|)\delta (tt^{}).$$
(11)
To better understand the source of this complication let us consider a simple mechanical analogy which is a Brownian motion theory.

The erratic motion of a Brownian particle is caused by collisions with the molecules of the fluid in which it moves. These collisions allow an exchange of the energy between the fluid and the Brownian particle. If the Brownian particle is much more massive than the molecules of the fluid then the influence of the molecules on the observed particle can be approximated by a Gaussian white noise $`\eta _G(t)`$:

$$m\ddot{x}(t)+\gamma \dot{x}(t)=\eta _G(t),$$
(12)
where $`x(t)`$ is the position of the Brownian particle. The generalization of Brownian motion theory to the random motion of a particle which is not necessarily heavier than the molecules of the fluid was proposed by Kubo . In this case the time scale of molecular motion is no longer very much shorter than that of the motion of the particle under observation, so that the random force $`\eta (t)`$ can not be of Gaussian type. To describe an influence of the molecules on the observed particle we have to introduce a color noise characterized by some time scale. This time scale may describe the average time interval between two subsequent collisions of the molecules with the observed particle. In addition, if we consider a stationary process we have to abandon the assumption of a constant friction and to introduce generally a frequency-dependent friction

$$m\ddot{x}(t)+_{t_0}^t𝑑t^{}\gamma (tt^{})\dot{x}(t^{})=\eta (t).$$
(13)
In case of $`\varphi ^4`$ model the distribution of impurities is not generally described by the white gaussian noise. In generic situation the distribution of impurities is characterized by some length scale which describes the average separation o impurity centers and therefore, in similar way as it was in case of Brownian particle, we introduce retardation to $`\varphi ^4`$ model.

Next we have to make further generalization of the Liu-Mazenko-Halperin formula. This generalization is achieved by averaging the formula (9) with respect to possible distributions of the impurity centers

$$n=\{n(t,x)\}=\frac{1}{\pi }\sqrt{\frac{\psi ^2+\text{{}u^2\text{}}}{\psi ^2+\text{{}u^2\text{}}}}.$$
(14)
The further understanding can be made for particular choice of the noise amplitude. The most representative one is Ornstein-Uhlenbeck amplitude $`W(|x|)=𝒜e^{\frac{|x|}{L}}`$ which interpolates between constant distribution and gaussian white noise. In this model the number density of produced kinks depends on quench time and on characteristic length scale of the impurity distribution as well

$$n\frac{1}{\pi }\sqrt{\frac{0.43\frac{b}{\tau ^{\frac{1}{2}}}+0.34\frac{c}{L^2}}{1.81b+0.83c}},$$
(15)
where $`b=\frac{1}{\tau ^{\frac{1}{4}}}`$ and $`c=\frac{L}{\sqrt{\tau }}`$. Let us notice that in case of week imperfections the usual scaling is recovered

$$n\frac{1}{\tau ^{\frac{1}{4}}}.$$
(16)
## 3 Remarks

We show that there are two components which determine the number density of kinks produced in the inhomogenous system during the second order phase transition. First component follows from the K-Z formalism for pure systems and is determined by the quench time. The second component is determined by the characteristic length which describes the distribution of impurities in the system. Due to existence of impurities and admixtures the kinks are created mainly in the knots of the impurity force distribution. One could even find exact solutions which describe the kinks confined by some particular impurity potentials . An examples of those solutions are squeezed kink or squeezed anti-kink. One could also check linear stability of those solutions. The other solution obtained for some particular force distribution is a static kink-antikink solution which illustrates that, in contradiction to pure systems, the configuration of this type have no tendency to annihilate. This is also a reason why at late time after transition in the inhomogenous system we still can find the number density of kinks substantially larger than estimated from the probability of thermal nucleation of kink-antikink pairs.

Acknowledgements I am indebted to organizers of the ”Third COSLAB and First Joint COSLAB - VORTEX - BEC2000+” workshop for hospitality and also for supporting my stay in Bilbao. Work supported in part by KBN grant 2 PO3B 025 25 and ESF ”COSLAB” Programe.