text:
# Thermal diffusion of supersonic solitons in an anharmonic chain of atoms

## I Introduction

Nonlinear one-dimensional lattice dynamics, namely propagation of coherent excitations in monatomic chains modeling discrete microscopic structures, is associated with several important problems in physics. Among these excitations are solitary waves, which for simplicity are called here solitons. These solitons can be supported by chains with realistic interaction potentials between the particles Lomdahl ; Lawrence . They are supersonic non-topological collective excitations. In spite of their relative simplicity, the solitons clarify many features of molecular chains Toda ; Pnevmatikos ; Davydov ; Collins ; Hochstrasser1 . For example, due to their robust character, lattice solitons have been used to model the energy transport in polypeptide chains in muscle proteins Yomosa ; Perez ; Perez1  or the energy transport in DNA Muto . Numerical simulations at realistic temperatures for transport in proteins have shown that lattice solitons can propagate over long distances in a chain with the Lennard-Jones potential Perez . Moreover, the lattice solitons are more stable than Davydov solitons if collisions between the two types of solitons are considered Perez1 . There is no clear evidence that lattice solitons like a Toda type, which are non-topological, can exist in thermal equilibrium. This holds both for static properties, like the specific heat, and for dynamics quantities, like the dynamic form factor (Fourier Transform of the displacement autocorrelation) Neuper . On the other hand, there exists evidence from real experiments that strain solitons can be generated and observed in non-linear elastic rods samsonov . These solitons in some cases can be described by Korteweg-de Vries (KdV) type solitons, which are non-topological.

To our knowledge there are no previous analytical studies supported by Langevin simulations about non-topological lattice soliton diffusion in anharmonic monatomic chains of particles with nearest-neighbor interactions. There are many studies on stochastic partial differential equations, in particular stochastic KdV-type equations have been extensively studied numerically and analytically KonotopVazquez ; Wadati ; Wadati2 ; Herman ; Iizuka ; ScalerandiRomano  due to the integrability of the KdV equation. In fact, the KdV equation is a good approximation to describe analytically the dynamics of lattice solitons on a monatomic chain with nearest-neighbor interaction and cubic anharmonicity if the soliton velocity is very close to the sound velocity (very-low-energy solitons) Pnevmatikos ; Remoissenet . Notice that for a polynomial potential, namely harmonic term plus cubic or/and quartic anharmonicity, the one-soliton solution of the KdV equation is known analytically, while for more realistic interaction potentials like Lennard-Jones or Morse there are no exact soliton solutions. In the more general context of lattice systems, there are a few analytical studies about diffusion of coherent lattice excitations, viz. stochastic vortex dynamics in two-dimensional easy-plane ferromagnets kampetter  or soliton diffusion on the classical, isotropic Heisenberg chain mathias1 ; mathias .

The aim of this work is to provide an approximate analytical description of the soliton diffusion dynamics in a monatomic chain with a cubic anharmonicity under thermal fluctuations. For this purpose we generate a single soliton which has an energy far greater than $`k_BT`$, where $`k_B`$ is the Boltzmann constant, and $`T`$ is the temperature of a thermal bath. This soliton propagating on a chain in contact with the thermal bath shows a diffusive behavior. We consider this soliton diffusion dynamics during the thermalization process of the system up to times when the system energy has relaxed nearly to its stationary value. This means that we study the non-equilibrium diffusion dynamics of lattice solitons on anharmonic chains subject to thermal fluctuations.

In order to perform the coupling between the system and a thermal bath with a given temperature an additive noise term, providing energy input, is added to the discrete equations of motion. This term has to be balanced by a damping term providing energy dissipation. Here, we suggest as a damping term the so-called hydrodynamical damping Landau  which is extensively used in, e.g., elasticity theory. Notice that this type of damping is due to irreversible processes taking place within the system. The corresponding noise term, which fulfills the fluctuation-dissipation theorem, takes the form of a discrete gradient of Gaussian white noise delta-correlated in space and time. A similar Langevin-type equations has previously been considered in the context of mesoscopic Langevin dynamics Schmittmann .

We notice that our system in the continuum limit can be approximated by a noisy KdV-Burgers-type equation nozaki ; Orlowski . So in this case we can use the one-soliton solution of the KdV equation not only as initial condition of our discrete system but also in our analytical approach in the continuum limit. Notice that the shape of broad KdV solitons tends to be identical to the shape of broad supersonic lattice solitons Pnevmatikos . In this work we apply a generalized traveling wave ansatz combined with a collective coordinate formalism in the framework of the KdV equation as an analytical approach to study the diffusion of lattice solitons.

In the next section we present the equations of motion of our discrete system. From this we formulate a set of stochastic equations of motion by adding noise and damping. Next, we apply the continuum limit and derive a form of noisy KdV-Burgers equation. In section III we apply a collective coordinate approach which yields analytical expressions for the thermal averages and variances of the soliton position and velocity. In section IV, we compare our analytical predictions with the results from Langevin dynamics simulations for the original discrete system. Our conclusions are summarized in the last section.

## II The continuum limit

We consider an anharmonic chain of particles with mass $`M`$ and nearest-neighbor interactions. The particles interact via an anharmonic potential with a cubic anharmonicity. The Hamiltonian of this system reads

$`H={\displaystyle \underset{n}{}}\{{\displaystyle \frac{P_n^2}{2M}}+`$
$`G({\displaystyle \frac{1}{2}}(Y_{n+1}Y_n)^2+{\displaystyle \frac{A}{3}}(Y_{n+1}Y_n)^3)\},`$ (1)
where $`Y_n`$ denotes the longitudinal displacement of the $`n`$-th particle from its equilibrium position, and

$$P_n=M\frac{dY_n}{dt}$$
(2)
is the momentum. Here $`G`$ and $`A`$ are the potential parameters whose values depend on the lattice. The associated first order equations of motion read

$`{\displaystyle \frac{dY_n}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{M}}P_n`$ (3)
$`{\displaystyle \frac{dP_n}{dt}}`$ $`=`$ $`{\displaystyle \frac{H}{Y_n}}+F_n^{Noise}+F_n^{Damping},`$ (4)
where

$`{\displaystyle \frac{H}{Y_n}}`$ $`=`$ $`G\left(Y_{n+1}2Y_n+Y_{n1}\right)`$ (5)
$`GA\left((Y_{n+1}Y_n)^2(Y_nY_{n1})^2\right).`$
In Eq. (4) we have already added both a stochastic force, $`F_n^{Noise}`$, and a damping force,$`F_n^{Damping}`$. Both forces couple the discrete system with a thermal bath. Here, we use the inner or hydrodynamical damping, which reads Landau ; arevalo

$$F_n^{Damping}=M\nu \left(\frac{dY_{n+1}}{dt}2\frac{dY_n}{dt}+\frac{dY_{n1}}{dt}\right).$$
(6)
This means that the energy dissipation is provided by the irreversible processes arising from the finite velocity of the internal motions of the system, namely time derivative of the relative displacements between particles in the chain. Eq. (6) is the discrete version of the damping used in elasticity theory Landau . To fulfill the fluctuation-dissipation theorem the noise must have the form (see App. A)

$$F_n^{Noise}=\sqrt{\mathrm{D}}\left(\xi _{n+1}(t)\xi _n(t)\right)$$
(7)
where

$$\mathrm{D}=2M\nu k_BT$$
(8)
is the diffusion constant and $`\nu `$ is the damping constant. $`\xi _n(t)`$ is delta-correlated white noise,

$`\xi _n(t)\xi _m(t^{^{}})`$ $`=`$ $`\delta _{nm}\delta (tt^{^{}}),`$ (9)
$`\xi _n(t)`$ $`=`$ $`0.`$ (10)
Since our interest is the study of the lattice soliton diffusion close to the sound velocity, $`c`$, we can use the continuum limit approach, where $`Y_n(t)y(x,t)`$ and $`\xi _n(t)\xi (x,t)`$ with $`x=na`$ and $`a`$ the equilibrium atomic spacing. In this limit Pnevmatikos , Eq. (3) reduces to a form of noisy and damped KdV equation (see App. B for details)

$$_\tau u+6u_su+_s^3u=\nu _1_{ss}u\sqrt{D_1}_s\overline{\xi }(s,\tau )$$
(11)
where

$$s=\alpha (xct),\tau =\beta t,u=\gamma _sy.$$
(12)
The constants $`\alpha `$, $`\beta `$ and $`\gamma `$ are defined in (87), and $`\nu _1`$ and $`D_1`$ are given by Eqs. (91). Note that

$`\overline{\xi }(s,\tau )\overline{\xi }(s^{},\tau ^{})=\delta (ss^{})\delta (\tau \tau ^{}).`$ (13)
Here and in the following the line over $`\overline{\xi }`$ is omitted.

The case $`D_1=0`$ reduces Eq. (11) to the KdV-Burgers equation. The associated KdV equation is

$$_\tau u+6u_su+_s^3u=0$$
(14)
whose one-soliton solution reads

$$u_0(s,\tau )=2\eta _0^2\text{sech}^2[\eta _0(s4\eta _0^2\tau s_0)].$$
(15)
Here

$$\eta _0=\frac{1}{p}\sqrt{3c(vc)}$$
(16)
is the inverse soliton width and $`s_0`$ is the initial soliton position. The sound velocity $`c`$ and the constant $`p`$ are defined in (B.1).

## III Collective coordinate approach

To analyze our problem we assume that the soliton profile, $`u_0(s,\tau )`$, is not disturbed by the noise and damping terms and that only the width and amplitude are modified. This assumption is well satisfied for low-energy solitons, whose velocity is close to the sound velocity, because tails induced by the perturbations are small in this velocity regime grim ; arevalo . So we introduce a generalized traveling wave ansatz of the form

$`u(s,\tau )`$ $`=`$ $`u_0(sS(\tau ),\eta (\tau ))`$ (17)
$`=`$ $`2\eta ^2(\tau )\text{sech}^2[\eta (\tau )(sS(\tau ))],`$
where the collective variables $`S(\tau )`$ and $`\eta (\tau )`$ are the soliton position and the inverse soliton width, respectively. Here and in the following the index of the one-soliton solution $`u_0`$ is omitted.

To obtain the equations for our collective coordinates we follow niurkaPRE ; Mertens1 . First, by substituting (17) into Eq. (11) we get

$$\varphi _1\dot{S}(\tau )+\varphi _2\dot{\eta }(\tau )=\nu _1_{ss}u\sqrt{D_1}_s\xi (s,\tau ),$$
(18)
where

$`\varphi _1(s,\tau )={\displaystyle \frac{u}{S}}`$ (19)
and

$`\varphi _2(s,\tau )={\displaystyle \frac{u}{\eta }}.`$ (20)
Notice that the functions $`\left\{\varphi _i\right\}_{i=1,2}`$ coincide with the adiabatic approximation (omitting secular terms in time) of the discrete solutions of the linearized KdV equation around the one-soliton solution (15) Mann . We remark here that our collective coordinate theory does not take into account the contribution of the phonon modes (continuous basis function solution of the linearized KdV). We discuss the effect of noise-induced phonons in section IV.2. The functions $`\left\{\varphi _i\right\}_{i=1,2}`$ are also orthogonal, so the inner product $`𝑑s\varphi _i(s,\tau )h(s,\tau )`$ projects a function $`h`$ onto the functions $`\left\{\varphi _i\right\}_{i=1,2}`$. Then, by projecting Eq. (18) we get

$$A_i\dot{S}(\tau )+B_i\dot{\eta }(\tau )=f_i+f_i^{damping}+f_i^{Noise}i=1,2,$$
(21)
where

$`A_i`$ $`=`$ $`{\displaystyle 𝑑s\frac{u}{S}\varphi _i},`$ (22)
$`B_i`$ $`=`$ $`{\displaystyle 𝑑s\frac{u}{\eta }\varphi _i},`$ (23)
$`f_i`$ $`=`$ $`{\displaystyle 𝑑s\left(6u_su+_s^3u\right)\varphi _i},`$ (24)
$`f_i^{damping}`$ $`=`$ $`\nu _1{\displaystyle 𝑑s_{ss}u\varphi _i},`$ (25)
$`f_i^{Noise}`$ $`=`$ $`\sqrt{D_1}{\displaystyle 𝑑s_s\xi (s,\tau )\varphi _i}.`$ (26)
After some calculations the Eqs. (21) take the form

$`{\displaystyle \frac{dS(\tau )}{d\tau }}=4\eta ^2(\tau )+{\displaystyle \frac{15\sqrt{D_1}}{64\eta ^5(\tau )}}{\displaystyle 𝑑s(_s\varphi _1)\xi (s,\tau )},`$ (27)
$`{\displaystyle \frac{d\eta (\tau )}{d\tau }}={\displaystyle \frac{30\nu _1}{30+\pi ^2}}\eta ^3(\tau )+`$
$`{\displaystyle \frac{45\sqrt{D_1}}{16(30+\pi ^2)\eta (\tau )}}{\displaystyle 𝑑s(_s\varphi _2)\xi (s,\tau )}.`$ (28)
To achieve the calculations we have assumed that the soliton profile remains mostly unaffected and only its width and amplitude change due to the stochastic perturbations. Then, at least for small noise, we can perform the calculations by taking the soliton field out of the averages. Moreover, we have interpreted Eqs. (27) and (III) in the Stratonovich sense, because it assumes $`\xi (s,\tau )`$ is a real noise with finite correlation time, which is then allowed to become infinitesimally small after calculating measurable quantities Gardiner . Notice that white noise means taking the limit of zero correlation time.

Eqs. (27) and (III) can take the form

$`{\displaystyle \frac{d𝐘(\tau )}{d\tau }}=𝐀^{Str}[𝐘(\tau )]+{\displaystyle 𝑑s\widehat{𝐁}^{Str}[s,𝐘(\tau )]𝝃(s,\tau )},`$ (29)
where the elements of the noise vector $`𝝃`$ satisfy (13). $`\{Y_1,Y_2\}=\{S,\eta \}`$ are the elements of the vector $`𝐘`$, the elements $`\{A_1^{Str},A_2^{Str}\}`$ of the drift vector $`𝐀^{Str}`$ are the drift terms in Eqs. (27) and (III), respectively. The diffusion matrix $`\widehat{𝐁}^{Str}`$ is diagonal, where $`B_{11}^{Str}`$ and $`B_{22}^{Str}`$ are the coefficients in front of the noise in Eqs. (27) and (III), respectively. In order to facilitate the calculations we write Eq. (29) in the Ito form,

$`d𝐘(\tau )=𝐀^{Ito}[𝐘(\tau )]d\tau +{\displaystyle 𝑑s\widehat{𝐁}^{Ito}[s,𝐘(\tau )]𝑑𝐖(s,\tau )},`$
where the $`d𝐖(s,\tau )=𝝃(s,\tau )d\tau `$ is a Wiener process. Via a Fokker-Planck equation, one can show that the elements of the drift vector $`𝐀^{Ito}`$ read Gardiner

$`A_i^{Ito}[𝐘(\tau )]=A_i^{Str}[𝐘(\tau )]+`$
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{jkm}{}}{\displaystyle 𝑑sB_{km}^{Str}[s,𝐘(\tau )]_{Y_k}B_{ij}^{Str}[s,𝐘(\tau )]}`$
$`i,j,k=1,2`$ (31)
while

$`\widehat{𝐁}^{Ito}[s,𝐘(\tau )]=\widehat{𝐁}^{Str}[s,𝐘(\tau )].`$ (32)
Notice that $`𝐀^{Ito}`$ and $`\widehat{𝐁}^{Ito}`$ are nonanticipating functions. So, from Eq.(III) it is easy to show the following averages

$`S(\tau )`$ $`=`$ $`{\displaystyle \underset{0}{\overset{\tau }{}}}𝑑\tau ^{}4\eta ^2(\tau ^{})`$
$`\eta (\tau )`$ $`=`$ $`{\displaystyle \underset{0}{\overset{\tau }{}}}𝑑\tau ^{}{\displaystyle \frac{30\nu _1}{30+\pi ^2}}\eta ^3(\tau ^{})+`$
$`{\displaystyle \underset{0}{\overset{\tau }{}}}𝑑\tau ^{}{\displaystyle \frac{225(231+8\pi ^2)D_1}{112(30+\pi ^2)^2}}`$
$`Var\left(S(\tau )\right)`$ $`=`$ $`{\displaystyle \underset{0}{\overset{\tau }{}}}𝑑\tau ^{}{\displaystyle \frac{75D_1}{112\eta ^3(\tau ^{})}}`$
$`Var\left(\eta (\tau )\right)`$ $`=`$ $`{\displaystyle \underset{0}{\overset{\tau }{}}}𝑑\tau ^{}{\displaystyle \frac{225(21+\pi ^2)D_1\eta (\tau ^{})}{28(30+\pi ^2)^2}}`$
$`Corr\left(S(\tau )\eta (\tau )\right)`$ $`=`$ $`0.`$ (33)
Here $`\mathrm{}`$ means average over an ensemble of realizations,

$$Corr(PQ)=PQPQ\text{and}Var(P)=Corr(PP).$$
(34)
Now we define a new set of Langevin equations,

$`dY_i(\tau )=a_id\tau +{\displaystyle \underset{j}{}}b_{ij}dW_j(\tau )`$
$`\text{with}i,j=1,2\text{and}\{Y_1,Y_2\}=\{S,\eta \},`$ (35)
which we have interpreted in the Ito sense. $`dW_j(\tau )=\xi _j(\tau )d\tau `$ are Wiener processes where we have let the noises to be uncorrelated, namely

$`\xi _j(\tau )\xi _j^{}(\tau ^{})=\delta (\tau \tau ^{})\delta _{jj^{}}.`$ (36)
In order to determine the values of $`a_i`$ and $`b_{ij}`$ we have demanded that Eqs.(III) satisfy the relations (III). It is straightforward to see that Eqs. (III) take the form

$`dS(\tau )`$ $`=`$ $`4\eta ^2(\tau )d\tau +{\displaystyle \frac{5\sqrt{3}}{4\sqrt{7}}}\sqrt{{\displaystyle \frac{D_1}{\eta ^3(\tau )}}}dW_1(\tau )`$ (37)
$`d\eta (\tau )`$ $`=`$ $`\left({\displaystyle \frac{30\nu _1}{30+\pi ^2}}\eta ^3(\tau )+{\displaystyle \frac{225(231+8\pi ^2)D_1}{112(30+\pi ^2)^2}}\right)d\tau +`$ (38)
$`{\displaystyle \frac{15\sqrt{21+\pi ^2}}{2\sqrt{7}(30+\pi ^2)}}\sqrt{D_1\eta (\tau )}dW_2(\tau ).`$
Eqs. (37) and (38) are statistically equivalent to Eqs (27) and (III) because they share the same Fokker-Planck equation. Since the derivation of Eqs (37) and (38) involved approximations, we have not solved them exactly. Instead of that, we have used perturbation analysis Abdullaev , which is developed in detail in App. C. In order to do so, we have considered the thermal terms as perturbations, so Eqs. (37) and (38) take the form

$`dS(\tau )`$ $`=`$ $`4\eta ^2(\tau )d\tau +ϵ{\displaystyle \frac{5\sqrt{3}}{4\sqrt{7}}}\sqrt{{\displaystyle \frac{D_1}{\eta ^3(\tau )}}}dW_1(\tau )`$ (39)
$`d\eta (\tau )`$ $`=`$ $`{\displaystyle \frac{30\nu _1}{30+\pi ^2}}\eta ^3(\tau )d\tau +ϵ({\displaystyle \frac{225(231+8\pi ^2)D_1}{112(30+\pi ^2)^2}}d\tau +`$ (40)
$`{\displaystyle \frac{15\sqrt{21+\pi ^2}}{2\sqrt{7}(30+\pi ^2)}}\sqrt{D_1\eta (\tau )}dW_2(\tau )).`$
Now, we seek an asymptotic solution in the form of a small-noise expansion

$`S(\tau )=s_0(\tau )+ϵs_1(\tau )+\mathrm{}`$
$`\eta (\tau )=\eta _0(\tau )+ϵ\eta _1(\tau )+\mathrm{}.`$ (41)
Here, $`ϵ`$ is a factor introduced for convenience in the analytical calculations. Notice that the case $`ϵ=0`$ reduces Eqs. (39) and (40) to the deterministic case. In order to interpret our perturbation theory we must set $`ϵ=1`$ and assume that the terms on the r.h.s. of Eqs. (39) and (40) are sufficiently small. So we must restrict ourselves to a regime of low temperatures of the thermal bath ($`D_1`$ small). From the perturbation analysis we obtain that

$`S(\tau )`$ $`=`$ $`s_0(\tau )+s_1(\tau )`$
$`=`$ $`4{\displaystyle \frac{\eta _0^2(0)}{\lambda }}\mathrm{log}(1+\lambda \tau )+{\displaystyle \frac{15D_1(231+8\pi ^2)\eta _0(0)(2(1+\lambda \tau )^{5/2}5\lambda \tau 2)}{7(30+\pi ^2)^2\lambda ^2(1+\lambda \tau )}},`$
$`4\eta ^2(\tau )`$ $`=`$ $`4\eta _0^2(\tau )+8\eta _0(\tau )\eta _1(\tau )`$
$`=`$ $`{\displaystyle \frac{4\eta _0^2(0)}{1+\lambda \tau }}+{\displaystyle \frac{45D_1(231+8\pi ^2)\eta _0(0)((1+\lambda \tau )^{5/2}1)}{7(30+\pi ^2)^2\lambda (1+\lambda \tau )^2}},`$
$`Var(S(\tau ))`$ $`=`$ $`D_1({\displaystyle \frac{15}{56\eta _0^3(0)\lambda }}{\displaystyle \frac{480(21+\pi ^2)\eta _0^3(0)(8+7\lambda \tau (5\lambda \tau +4))}{49\left(30+\pi ^2\right)^2\lambda ^3\left(1+\lambda \tau \right)^2}}`$
$`+`$ $`{\displaystyle \frac{15(1+\tau \lambda )^{3/2}}{392\left(30+\pi ^2\right)^2\eta _0^3(0)\lambda ^3}}(2048(21+\pi ^2)\eta _0^6(0)+7(30+\pi ^2)^2\lambda ^2(1+\lambda \tau ))),`$
$`Var(4\eta ^2(\tau ))`$ $`=`$ $`{\displaystyle \frac{7200D_1(21+\pi ^2)\eta _0^3(0)}{49\lambda (30+\pi ^2)^2}}\left({\displaystyle \frac{1}{\sqrt{1+\lambda \tau }}}{\displaystyle \frac{1}{(1+\lambda \tau )^4}}\right).`$
The expressions for $`D_1`$ and $`\lambda `$ are given by Eqs. (91) and (101), respectively.

## IV Simulations

Substituting Eqs. (6) and (7) into Eq. (3) we get the full set of discrete equations of motion written in absolute displacements. However, for our simulations relative displacements are more convenient, because the lattice solitons in this representation are pulse solitons whose amplitude vanishes at infinity. This characteristic allows us to use periodic boundary conditions which are necessary for long simulation times, because we want to avoid reflections at the boundaries. So the discrete equations of motion in relative displacements read

$`M{\displaystyle \frac{d^2V_n}{dt^2}}=G\left(V_{n+1}2V_n+V_{n1}\right)+`$
$`GA\left(V_{n+1}^22V_n^2+V_{n1}^2\right)+`$
$`M\nu \left({\displaystyle \frac{dV_{n+1}}{dt}}2{\displaystyle \frac{dV_n}{dt}}+{\displaystyle \frac{dV_{n1}}{dt}}\right)+`$
$`\sqrt{\mathrm{D}}\left(\xi _{n+1}(t)2\xi _n(t)+\xi _{n1}(t)\right),`$ (43)
where $`V_n(t)=Y_{n+1}(t)Y_n(t)`$ and $`\mathrm{D}=2M\nu k_BT`$. The periodic boundary conditions read

$`{\displaystyle \frac{d^lV_0}{dt^l}}={\displaystyle \frac{d^lV_{N1}}{dt^l}},`$ $`{\displaystyle \frac{d^lV_N}{dt^l}}={\displaystyle \frac{d^lV_1}{dt^l}},l=0,1`$
$`\xi _0(t)=\xi _{N1}(t),`$ $`\xi _N(t)=\xi _1(t),`$ (44)
where $`N`$ is the number of particles of our chain and $`N1`$ is the number of bonds.

A suitable method to detect the position of a pulse lattice soliton, $`V_n(t)`$, is to search for its maximum arevalo . However, in the presence of stochastic perturbations this method is not useful since the pulse shape is strongly masked by the noise, an example of this situation is shown in App. E. So from the data of our simulations we have taken snapshots of the system at different times, and from them we have generated the kink shape $`Y_n(t)`$ of the lattice soliton by using the algorithm

$$Y_n(t)=Y_1(t)+\underset{i=1}{\overset{n1}{}}V_i(t)n=2,3,\mathrm{},N.$$
(45)
The kink shape is less distorted by the noise than the pulse shape $`V_i(t)`$. In (45) $`Y_1(t)`$ is a boundary condition that we have demanded to be

$$Y_1(t)=\frac{1}{2}\underset{i=1}{\overset{N1}{}}V_i(t),$$
(46)
so at $`t=0`$ the amplitude of the center of the kink shape is zero, as it should be from the theory Pnevmatikos . Notice that

$$Y_N(t)Y_1(t)=\underset{i=1}{\overset{N1}{}}V_i(t)$$
(47)
is a conserved quantity in our system, i.e.

$$\dot{Y}_N(t)\dot{Y}_1(t)=\underset{i=1}{\overset{N1}{}}\dot{V}_i(t)=0.$$
(48)
We have checked Eq. (48) with a precision higher than $`10^{14}`$ over the whole time range of our Langevin dynamics simulations.

In order to determine every time the parameters of the soliton, namely soliton velocity $`v`$ and position $`x`$, we have proceeded as follows. We have searched for the values of the parameters $`x`$ and $`v`$ where the relation

$`{\displaystyle \underset{n=n_0n_1}{\overset{n_0+n_1}{}}}\left(Y_ny_0(nax,v)\right)=0`$ (49)
is fulfilled. Here

$$y_0(nax,v)=\frac{6\sqrt{2hc(vc)}}{p}\mathrm{tanh}\left(\frac{nax}{L(v)}\right)$$
(50)
with

$$L(v)=(\alpha \eta (t))^1=\mathrm{\hspace{0.17em}2}\sqrt{\frac{h}{2c(vc)}}.$$
(51)
Here the function $`\eta `$ is defined in (16) and the constants $`h`$, $`p`$, $`c`$ and $`\alpha `$ are defined in (B.1) and (87). the function (50) is the absolute displacement representation of the one-soliton solution (15) in a frame moving with the soliton velocity. In Eq. (49)

$$n_0=\text{int}(x)\text{and}n_1=\text{int}(\frac{3}{2}L(v(0))),$$
(52)
where $`\text{int}()`$ denotes the integer part of a number and $`v(0)`$ is the initial soliton velocity. The value of $`n_1`$ has been chosen to take into account only the core of the lattice kink-shape and it is constant during our simulations. In order to determine both $`x`$ and $`v`$ we consider different values of $`v`$ in Eq. (49) within a range of velocities around the initial soliton velocity, namely $`vc[0.1(v(0)c),2(v(0)c)]`$. For every value of $`v`$ we search the value $`x`$ that fulfills Eq. (49), so we get a set of pairs of values $`x`$ and $`v`$. Finally, from this set of pairs of values we search, by using linear interpolation, the values of $`x`$ and $`v`$ which fulfill the relation

$`{\displaystyle \frac{\underset{n=n_0n_1}{\overset{n_0+n_1}{}}Y_ny_0(nax,v)}{\underset{n=n_0n_1}{\overset{n_0+n_1}{}}(y_0(nax,v))^2}}=1.`$ (53)
Notice that in Eqs. (49) and (53) we have assumed that the lattice kink shape, $`Y_n`$, is closely related with the function $`y_0`$, however, as was mentioned in Ref. arevalo , a pulse lattice soliton in the presence of damping develops a tail. The amplitude of this trailing tail depends on both soliton the velocity and the damping, so it is bigger when the damping and/or the soliton velocity is higher. Thus, we restrict ourselves to velocities very close to the sound velocity where the effect of this trailing tail is negligible.

Up to now we have determined the parameters $`x`$ and $`v`$, which fit the function $`y_0`$ to the lattice kink shape $`Y_n`$, so we have not measured directly either $`x`$ or $`v`$. Since the function $`y_0`$ is closely related to the lattice kink-shape $`Y_n`$, one could assume both $`x`$ and $`v`$ as an estimate of the soliton position and velocity, respectively. However, we have taken only the parameter $`v`$ as an estimate of the soliton velocity and with this value we have used a different method to determine the soliton position. In fact, in order to be in agreement with our collective coordinate approach, where we have projected the equations of motion onto the Goldstone mode $`\varphi _1`$ (Eq. (19)), we have projected the noisy kink shape $`Y_n(t)`$ onto the pulse solution $`u_0`$ defined in (15). Notice that in the absolute displacement representation the function $`u_0`$ is the Goldstone mode. So this projection reads

$$P(x)=\underset{i=nn_2}{\overset{n+n_2}{}}Y_i(t)u_0(iax,v)$$
(54)
where

$$u_0(iax,v)=\frac{6c(vc)}{p}\text{sech}^2\left(\frac{iax}{L(v)}\right)$$
(55)
and $`x=na`$. The value of $`n_2`$ in Eq. (54) is much larger than the soliton width, so the boundary effects are negligible. The function $`u_0(iax,v)`$ is the one-soliton solution (15) in a frame moving with the soliton velocity. Afterwards we have searched, by linear interpolation, the value $`x`$ where $`P(x)`$ vanishes and we have defined it as the position of the soliton center of mass. At this point we remark that the values of $`x`$ following from this latter method are not significantly different from the values of $`x`$ following from the former method (Eqs. 49 and 53). However, we consider the latter method to be more appropriate than the former one in the sense that we proceed in our code in a similar way as in our analytical calculations.

Our Langevin dynamics simulations were performed for a chain with 1500 lattice points. The time integration was carried out by using the Heun method Heun , which has been successfully used in the numerical solution of partial differential equations and difference-differential equations, coupled to either an additive or a multiplicative noise term kampetter ; niurka1 ; mathias1 ; mathias . Here, we have used the conserved quantity (47) to check the accuracy of our code arevalo . For the longest simulation time the variation of this conserved quantity has been lower than $`10^9\%`$. In order to start the simulations at $`t=0`$ we have used the one-soliton solution (15) of the KdV equation in the laboratory frame. The average values have been calculated over 200 realizations up to a final time 5000. All values of the constants of the equation (IV) are set at unity except the damping constant which is set at $`\nu =0.003`$. Notice that for lower values of damping the relaxation of the system energy would take more time in our simulations to reach a regime close to its stationary value. On the other hand, higher values of damping can strongly distort the soliton shape, namely the tail induced by the damping cannot be neglected when the value of the damping is high. In App. D we show the thermalization process in our system. In our simulations the values of temperature, $`T`$, and initial soliton velocity, $`v(0)`$, are parameters (see figure captions).

### IV.1 Soliton propagation

In Figs. 1a and 1b we show several examples of both the averaged soliton position, $`x(t)`$, and the averaged soliton velocity, $`v(t)`$, as functions of time from both the simulation and the theory (see Eqs.(III)). Notice that

$`x(t)={\displaystyle \frac{1}{\alpha }}S(\beta t)+ct`$ (56)
and

$`v(t)c=\left({\displaystyle \frac{p^2}{12c}}\right)4\eta ^2(\beta t).`$ (57)
For all cases the soliton position from the simulation agrees well with the position given by the analytical theory. In the case of the soliton velocity, the agreement is better for initial velocities close to the sound velocity than for higher velocities. In fact, the time evolution of the velocity given by the simulation is always higher than the theoretical prediction. This systematic difference may be due to the small tail that is generated by the lattice soliton since it is a non-topological soliton arevalo . So it may affect our numerical method for determining the soliton position. We point out that the amplitude of this tail depends on both the soliton velocity and the damping constant, and can be neglected only for velocities close to the sound velocity and small values of the damping constant arevalo . This is the most important reason for restricting our study to low-energy solitons whose velocities are close to the sound velocity. Since this difference is systematic it does not play any role in the numerical calculation of the variances which is our more important goal in this article.

### IV.2 Soliton diffusion

In Figs. 2 and 3 we show the variances of the soliton position and velocity vs. time for different initial velocities. The temperature of the thermal bath in Fig. 3 is 10 times higher than that in Fig. 2. The results scale very well by a factor of 10; i.e. the variances are proportional to the temperature, as expected from Eqs. (III).

Notice that our theory (solid lines) has no adjustable parameters. Taking into account that the theory consists of several steps (discrete system $``$ Bq equation $``$ KdV equation $``$ collective coordinate theory $``$ perturbation analysis) it is already a significant success to obtain the observed results, as seen in Figs. 2 and 3

We observe that the behavior of the variances depends strongly on the initial soliton velocity. For low-energy solitons, whose velocities are close to the sound velocity (Figs. 2 and 3, cases (a) and (b)), the soliton diffusion tends to be nearly normal, i.e. linear in time. In fact, our theory predicts a normal diffusion for times

$$t<<t^{}=\left(\frac{30+\pi ^2}{495}\right)\left(\frac{c}{\nu (v(0)c)}\right).$$
(58)
This estimate was obtained by comparing the first with second terms of the Taylor expansion in powers of $`\tau `$ of the variance of the velocity (see Eq. (C)). For low-energy solitons ($`v(0){}_{}{}^{>}c`$) $`t^{}`$ is much larger than our simulation time (Figs. 2a and 3a). This means that the superdiffusivity is not very dominant. However, for higher-energy solitons the anomalous behavior turns out to be important after some time. In those cases $`t^{}`$ is comparable with our simulation time (Figs. 2e and 3e).

In the case of solitons with very low energy the variance of the position (Figs. 2a and 3a) does not agree so well with the theoretical prediction. In fact, we observe a transient behavior for times $`t{}_{}{}^{<}3/\nu =1000`$ where the system energy shows a fast relaxation process (see App. D). Those discrepancies between theory and simulations may be due to the combination of two effects. First, the profile of low-energy solitons is strongly masked by the noise, so the numerical detection of the position can be distorted. Second, not only the noise but also the noise-induced phonons can make a significant contribution to the variance of the position, since the reduced temperature, $`\overline{T}=k_BT/H(0)`$ (temperature in units of the initial soliton energy $`H(0)`$), of the thermal bath is higher here than in the other cases (see captions Figs. 2 and 3). In this respect we estimate the phonon effect on the diffusion of low-energy solitons in the next section (see also Fig. 4).

Notice that the variance of the soliton position is larger for low-energy solitons (Figs. 2a and 3a) than for higher-energy solitons (for instance Figs. 2e and 3e). This is due to the fact that the higher-energy solitons are more robust against thermal fluctuations than the lower-energy ones. Or, equivalently, the reduced temperature $`\overline{T}`$ of the thermal bath is higher for slow solitons than for the fast ones (see captions of Figs. 2 and 3).

On the other hand, the superdiffusive behavior is more pronounced for higher-energy solitons. This is because the soliton velocity turns out to be more sensitive to the thermal fluctuations in this case than in the case of broader solitons. Notice that the soliton velocity and soliton width are related. Also, since the higher-energy solitons encompass few lattice sites, the soliton-width perturbations are larger with respect to the averaged soliton width in this case than in the case of broader solitons (low-energy solitons). In fact, the variance of the soliton velocity shows this effect, namely that for broader solitons (Figs. 2b and 3b) this variance is smaller than for narrower solitons (for instance Figs. 2f and 3f). The discrepancy between our theory and the numerical simulations for higher-energy solitons (Figs. 2 and 3, cases (c)-(f)) is mainly due to the fact that our theory is valid only for soliton velocities close to the sound velocity.

With respect to the variance of the soliton velocity (Figs. 2 and 3, panels b, d and e) we observe that it is mostly anomalous and its behavior is nearly quantitatively predicted by our theory for $`0t{}_{}{}^{<}\mathrm{\hspace{0.17em}2000}`$ in all the cases. For larger times, $`2000t<\mathrm{\hspace{0.17em}5000}`$, there is a discrepancy which becomes larger with increasing initial velocities.

We remark that the numerical results shown in Figs. 2 and 3 do not change for systems with the double number of sites, namely 3000.

We comment that there was a previous attempt by Scalerandi et al. ScalerandiRomano  to calculate theoretically the mean square displacement of a KdV soliton subject to stochastic fluctuations. They considered the case of small Stokes damping and a simple white noise delta-correlated in time and space. Though their theoretical result shows the appearance of a noise-induced superdiffusive behavior, it does not have the same dependence with respect to the soliton width as our results (III), which agree well with our simulations.

### IV.3 Estimate of the phonon contribution

In order to estimate this contribution we have performed the following numerical test. We have simulated the propagation of low-energy solitons under thermal fluctuations up to a time (e.g. 2500) when the system energy is close to its stationary value (see App. D). Notice that the diffusion of low-energy solitons is mostly normal (see Figs 2a and 3a). Afterwards, we have isolated the system from the thermal bath by switching off noise and damping, so that solitons propagate only in the noise-induced phonon bath. The diffusion is mostly normal before and after the system is isolated, i.e. the variance of the position is linear in time. We have compared the slope of the variance of the isolated system ($`t>\mathrm{\hspace{0.17em}2500}`$) with the slope of the variance in the case when the system is in contact with the thermal bath the whole time. An example of this test is shown in Fig. 4. We observe that both slopes, after switching off noise and damping, are different. In a normal diffusion process the slope of the variance of the soliton position is the diffusion constant. We call this $`D_{total}`$ when the system is in contact with the thermal bath since there is a contribution of both the thermal fluctuations and the phonons. The diffusion constant due to the noise-induced phonon bath (isolated system) is termed $`D_{ph}`$. On the other hand, our theoretical diffusion constant, $`D_{th}`$, is defined as the linear coefficient of the Taylor expansion of $`Var(x(t))`$ (see Eqs. (C) and (105)). Since our theory does not take into account the contribution of phonon modes, we expect that the value $`D_{noise}=D_{total}D_{ph}`$ may be of the same order of magnitude of $`D_{th}`$. We observe in Fig. 5 that the relative deviation of $`D_{th}`$ from $`D_{noise}`$, $`((D_{noise}D_{th})/D_{noise}`$, has the same order of magnitude of $`D_{th}`$ which is not surprising since the soliton shape is strongly masked and distorted by the noise (see App. E). In this respect we remark that from our results in Figs. 2 and 3 we observe that

$$\sqrt{Var(x(t))}<L(v(t)),$$
(59)
where the soliton width $`L(v(t))`$ is defined in (51). The relation (59) means that the stochastic deviations of the soliton center from its mean value are relatively small compared with the soliton width. So the diffusive dynamics of the soliton position evolves inside the soliton core. Thus, this dynamics is very sensitive to the fluctuations of the soliton shape. Notice that our method of determining the soliton position depends implicitly on the soliton shape. So in the case of low-energy solitons, where the shape is strongly masked by the noise, the uncertainties of the method of soliton detection are relatively large compared with the diffusive dynamics of the soliton position. This is because the diffusive dynamics is relative small with respect to the soliton width. In this respect we note that higher-energy solitons present a well defined shape, i.e. the thermal fluctuations are small respect to the soliton amplitude, so the variance of the soliton position, given by our detection method, indeed agrees better with our theory. (see 2c and 3c). Finally, we stress that in our tests the phonon contribution to the soliton diffusion could be clearly observed only for very low-energy solitons (Figs. 2 and 3, cases (c) and (e)), for higher-energy solitons the effect is negligible, namely $`D_{ph}0`$.

### IV.4 Estimate of the bath temperature in real physical systems

To clarify the physical meaning of the obtained results we have estimated the soliton energy and characteristic temperature of the thermal bath for two systems where solitons are believed to play an essential role: $`\alpha `$-helical proteins Perez  and crystal inertial gases ballistic . In both cases the Lennard-Jones (LJ) interaction potential,

$$\varphi _{LJ}=4E_0\left(\left(\frac{a}{r}\right)^{12}\left(\frac{a}{r}\right)^6\right),$$
(60)
is used. One can estimate the potential parameters $`G`$ and $`A`$ of the Hamiltonian (II) by using Taylor expansion of (60) around the minimum $`r=(12^{1/6})a`$. So, taking into account only the coefficients of the harmonic and cubic terms of this expansion, one can determine the potential parameters, namely

$`A={\displaystyle \frac{21}{2^{7/6}a}}\text{and}G={\displaystyle \frac{\mathrm{36\hspace{0.17em}2}^{2/3}E_0}{a^2}}.`$ (61)
By using the one-soliton solution of the Bq equation (79) Pnevmatikos  and performing an integration instead of a summation in Eq. (II), the initial soliton energy $`H(0)`$ reads

$$H(0)=\frac{16\sqrt{3}}{245}((v/c)^21)^{3/2}(1+9(v/c)^2)E_0.$$
(62)
Here $`v`$ and $`c`$ are the soliton and sound velocities, respectively. $`E_0=0.22eV`$ for $`\alpha `$-helix Perez  or $`E_0=1.0\times 10^2eV`$ for argon argon . So, for example, in the cases of $`v/c`$=1.003 or 1.03 we get

|  | $`\alpha `$-helix | argon |
| --- | --- | --- |
| $`v/c`$ | $`H(0)/k_B`$ \[K\] | $`H(0)/k_B`$ \[K\] |
| 1.003 | 1.4 | 0.06 |
| 1.03 | 45.7 | 2.15 |

Here $`k_B`$ is the Boltzmann constant. Then, the value of the temperature of the thermal bath can be obtained by multiplying the values of $`H(0)/k_B`$ by the reduced temperature. For instance, in the $`\alpha `$-helix case the reduced temperature $`\overline{T}=0.061419`$ (see caption Fig. 3) corresponds to $`0.1K`$ of the thermal bath if $`v(0)=1.003`$ or $`3.8K`$ if $`v(0)=1.03`$.

## V Summary and conclusions

We have studied the non-equilibrium diffusion dynamics of lattice solitons on a classical chain of atoms under thermal fluctuations, namely soliton dynamics when the system energy is close to its stationary value. The interaction potential between the atoms is harmonic plus a cubic anharmonicity. The chain is coupled to a thermal bath at a given temperature. For that reason we have included dissipation and noise in the discrete equations of motion of the chain. Here, it is assumed that the energy dissipation is provided by the irreversible processes arising from the finite velocity of the relative displacements between particles in the chain. Thus the dissipative term takes the form of a hydrodynamical damping which is extensively used in elasticity theory. The noise term which fulfills the fluctuation-dissipation theorem then becomes a discrete gradient of white noise delta-correlated in space and time. In the continuum approach our original discrete set of equations leads to a form of noisy KdV-Burgers equation. At this point we have used a collective coordinate approach to study the diffusion dynamics of both position and velocity of the soliton. The soliton position and the inverse soliton width have been found to be good collective coordinates to describe the soliton diffusion. We have derived two stochastic ordinary differential equations with multiplicative noise which have been solved analytically using stochastic perturbation analysis.

For low-energy solitons, whose velocities are close to the sound velocity, our molecular dynamics simulation has confirmed our analytical predictions. Namely, normal diffusion of lattice solitons governs short times, while superdiffusive behavior is present for long times. The time range of the normal diffusion depends on the initial velocity of the soliton: it is large for velocities close to the sound velocity and short for high velocities. The collective coordinate approach does not take into account the noise-induced phonon bath, however we have shown that this does not play an important role except when the reduced temperature (temperature in units of the initial soliton energy $`H(0)`$) of the thermal bath is high. In that regime the soliton diffusion is normal. In this case, for a given temperature, we have estimated in our simulations the value of the diffusion constant due to the noise-induced phonon bath when the system energy is close to its stationary value. We have subtracted this value from the full value of the diffusion constant which is not only due to the induced phonons but also due to the noise. The order of magnitude of the resultant value of this subtraction is predicted by the collective coordinate approach.

Since we do not observe in our numerical results any dependence on the size of the system, we may expect similar results for very large systems, i.e. for $`N>>1500`$.

We provided an example by using an approximation of the Lennard-Jones potential to determine the temperature of the thermal bath in the cases of $`\alpha `$-helical proteins and crystal inertial gases.

Finally, our results above point out the robustness of lattice solitons. In fact, they can exist even for higher values of temperature and damping constant than those explicitly considered in the present article. On the other hand, for lower values of the temperature the variances of the soliton position and velocity turn out to be very small because they scale with the temperature. So it is very difficult to observe them.

Acknowledgements

We acknowledge support from DLR grant Nr.: UKR-002-99. Yu. Gaididei is grateful for the hospitality of the University of Bayreuth where this work was performed. E. Arévalo acknowledges support from Eur. Grad. School “Non-Equilibrium Phenomena and Phase Transitions in Complex Systems”. A. R. Bishop is a Humboldt Awardee at the University of Bayreuth.

## Appendix A Derivation of the noise term

The goal of this Appendix is to find the form of the noise force which would satisfy the fluctuation-dissipation theorem. The associated set of Langevin equations of the classical chain of atoms under thermal fluctuations are

$`{\displaystyle \frac{dP_n}{dt}}`$ $`=`$ $`T_n+F_n^{Noise}+F_n^{Damping},`$ (63)
$`{\displaystyle \frac{dY_n}{dt}}`$ $`=`$ $`{\displaystyle \frac{P_n}{M}},`$ (64)
where

$`T_n`$ $`=`$ $`{\displaystyle \frac{H}{Y_n}}={\displaystyle \frac{U}{Y_n}},`$
$`F_n^{Damping}`$ $`=`$ $`M\nu \left({\displaystyle \frac{dY_{n+1}}{dt}}2{\displaystyle \frac{dY_n}{dt}}+{\displaystyle \frac{dY_{n1}}{dt}}\right)`$ (65)
and $`F^{Noise}(t)`$, which satisfies the fluctuation-dissipation theorem, is to determined. $`P_n`$ is the momentum, $`Y_n`$ denotes longitudinal displacement from its equilibrium position, of $`nth`$ particle with mass $`M`$ and velocity $`dY_n/dt`$. $`H`$ is the Hamiltonian

$$H=K+U,$$
(66)
where

$`K={\displaystyle \underset{n}{}}{\displaystyle \frac{P_n^2}{2M}},U={\displaystyle \underset{n}{}}V[Y_{n+1}Y_n]`$ (67)
and $`V`$ is an arbitrary potential which depends on the relative displacements $`Y_{n+1}Y_n`$. The discrete Fourier transform of Eqs. (63) and (64) read

$`{\displaystyle \frac{d\stackrel{~}{P}_k}{dt}}`$ $`=`$ $`\stackrel{~}{T}_k\nu \stackrel{~}{\gamma }_k\stackrel{~}{P}_k+\stackrel{~}{F}_k^{Noise}(t)`$
$`{\displaystyle \frac{d\stackrel{~}{Y}_k}{dt}}`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{P}_k}{M}},`$ (68)
where

$$\stackrel{~}{\gamma }_k=2\left(1\mathrm{cos}(k)\right).$$
(69)
We define

$$\stackrel{~}{F}_k^{Noise}(t)=\sqrt{D(k)}\stackrel{~}{\xi }_k(t),$$
(70)
where $`\stackrel{~}{\xi }_k(t)`$ is delta-correlated white noise,

$$\stackrel{~}{\xi }_k(t)\stackrel{~}{\xi }_k^{}(t^{})=D(k)\delta (tt^{})\delta _{k,k^{}},$$
(71)
and $`D(k)`$ is unknown.

The associated Fokker-Planck equation of Eqs.(A) in the Stratonovich sense takes the form

$`_t\rho `$ $`=`$ $`{\displaystyle \underset{k}{}}(_{\stackrel{~}{P}_k}(T_k\rho _k){\displaystyle \frac{\stackrel{~}{P}_k}{M}}_{\stackrel{~}{Y}_k}\rho _k+`$ (72)
$`\nu \stackrel{~}{\gamma }_k_{\stackrel{~}{P}_k}(\stackrel{~}{P}_k\rho _k+{\displaystyle \frac{D(k)}{2\nu \stackrel{~}{\gamma }_k}}_{\stackrel{~}{P}_k}\rho _k)),`$
where

$$\rho =\underset{k}{}\delta (\stackrel{~}{P}_k\stackrel{~}{P}_k(t))\delta (\stackrel{~}{Y}_k\stackrel{~}{Y}_k(t)).$$
(73)
In order to determine $`D(k)`$ we have demanded the stationary solution of Eq.(72) to be the Boltzmann distribution, namely

$$\rho =𝒩\mathrm{exp}\left(\frac{H}{k_BT}\right),$$
(74)
where $`H`$ is defined in Eq. (66) and $`𝒩`$ is the normalization constant. Substituting Eq.(74) into Eq.(72), it is straightforward to see that

$$D(k)=2\nu \stackrel{~}{\gamma }_kMk_BT.$$
(75)
Therefore, from Eq. (70) together with (75), it is easy to show that in position space

$`F_n^{Noise}(t)F_n^{}^{Noise}(t^{})=`$
$`2\nu Mk_BT\delta (tt^{})(\delta _{n+1,n^{}}2\delta _{n,n^{}}+\delta _{n1,n^{}}).`$ (76)
Finally, the relation (A) can be satisfied by the definition

$$F_n^{Noise}(t)=\sqrt{2\nu Mk_BT}(\xi _{n+1}(t)\xi _n(t)),$$
(77)
where

$$\xi _n(t)\xi _n^{}(t^{})=\delta (tt^{})\delta _{n,n^{}}.$$
(78)
## Appendix B Continuum limit

In order to reduce Eq. (3) to a form of noisy and damped KdV equation we have performed two steps. First, we have employed the continuum approach Pnevmatikos ; Remoissenet  in order to obtain a form of noisy and damped Bq equation, and then we have used the reductive perturbation technique Taniuti ; Remoissenet  in order to obtain the noisy and damped KdV equation.

### B.1 Noisy and damped Bq equation

Here we have used the procedure of Pnevmatikos Pnevmatikos , who expanded $`Y_{n\pm 1}(t)`$ and $`Y_{n\pm 2}(t)`$ in a Taylor series around $`y(x,t)`$, with $`x=na`$, where the equilibrium atomic spacing $`a`$ is regarded as an expansion parameter. Then, collecting powers of $`a`$, Eq. (3) together with Eqs. (6) and (7) at $`O(a^4)`$ takes the form

$$_t^2y=c^2_x^2y+p_xy_x^2y+h_x^4y+\nu a^2_x^2_ty+a\sqrt{D}_x\xi (x,t),$$
(79)
where

$$\frac{\xi _{n+1}(t)\xi _n(t)}{a^{3/2}}_x\xi (x,t),$$
(80)
with properties

$`_x\xi (x,t)`$ $`=`$ $`0,`$
$`_x\xi (x,t)_x^{}\xi (x^{},t^{})`$ $`=`$ $`_x^{}_x\delta (xx^{})\delta (tt^{}).`$ (81)
The diffusion constant takes the form

$$D=\frac{2\nu k_BT}{\rho }.$$
(82)
Other constants are

$`c^2={\displaystyle \frac{Ga}{\rho }},`$ $`p={\displaystyle \frac{2a^2AG}{\rho }},`$
$`h={\displaystyle \frac{a^3G}{12\rho }},`$ $`\rho =M/a.`$ (83)
### B.2 Noisy and damped KdV equation

We write Eq. (79) in the sound velocity frame and make a further approximation concerning variations in time. In this case we may use the reductive perturbation technique Taniuti . Therefore we rewrite Eq. (79) in the perturbation form

$`_t^2yc^2_x^2yp_xy_x^2yh_x^4y=`$
$`\kappa \left(\nu a^2_x^2_ty+a\sqrt{D}_x\xi (x,t)\right),`$ (84)
where we have introduced a small parameter $`\kappa `$ for convenience. Afterwards, we perform the following change of variables

$$s=\kappa \alpha (xct),\tau =\kappa ^3\beta t,u=\gamma _sy,$$
(85)
so that

$`_x=\kappa \alpha _s\text{and}_t=\kappa ^3\beta _\tau \kappa \alpha c_s,`$ (86)
where

$$\alpha =\frac{p}{\sqrt{6h}},\beta =\frac{p^3}{12c\sqrt{6h}},\gamma =\frac{1}{\sqrt{6h}}.$$
(87)
We have also expressed $`u`$ and $`\xi `$ in a perturbation series

$`u=\kappa u_1+\kappa ^2u_2+\mathrm{}`$ (88)
$`\xi =\kappa \xi _1+\kappa ^2\xi _2+\mathrm{}`$ (89)
Here the parameter $`\kappa `$ indicates the magnitude of the rate of change, the coefficients $`\kappa `$ and $`\kappa ^3`$ in (85) are chosen in order to balance the nonlinear term, and the dispersive term and the time derivative are of the same order in $`\kappa `$. The small noise expansion (89) is defined such that the lowest order terms of noise and damping are of the same order. Substituting Eqs. (86), (88) and (89) into (B.2), and keeping the lowest order terms, namely $`O(\kappa ^5)`$, we find that

$$_\tau u+6u_su+_s^3u=\nu _1_{ss}u\sqrt{D_1}_s\xi (s,\tau ),$$
(90)
where we have set $`u=u_1`$ and $`\xi =\xi _1`$.

$`\nu _1={\displaystyle \frac{\sqrt{6}\nu a^2c}{\sqrt{h}p}},D_1=D\alpha \beta \left({\displaystyle \frac{6a}{p^3}}\right)^2.`$ (91)
## Appendix C Perturbation analysis

In this Appendix we develop a perturbation approach to the equations (see Eqs (39) and (40))

$`dS(\tau )`$ $`=`$ $`4\eta ^2(\tau )d\tau +ϵ{\displaystyle \frac{5\sqrt{3}}{4\sqrt{7}}}\sqrt{{\displaystyle \frac{D_1}{\eta ^3(\tau )}}}dW_1(\tau )`$ (92)
$`d\eta (\tau )`$ $`=`$ $`{\displaystyle \frac{30\nu _1}{30+\pi ^2}}\eta ^3(\tau )d\tau +ϵ({\displaystyle \frac{225(231+8\pi ^2)D_1}{112(30+\pi ^2)^2}}d\tau +`$ (93)
$`{\displaystyle \frac{15\sqrt{21+\pi ^2}}{2\sqrt{7}(30+\pi ^2)}}\sqrt{D_1\eta (\tau )}dW_2(\tau )).`$
We interpret Eqs. (92) and (93) in the Ito sense where the Wiener process $`dW_i(\tau )=\xi _i(\tau )d\tau `$ with $`\xi _i(\tau )\xi _j(\tau ^{})=\delta _{ij}\delta (\tau \tau ^{})`$. We seek an asymptotic solution of the form

$`S(\tau )=s_0(\tau )+ϵs_1(\tau )+\mathrm{}`$
$`\eta (\tau )=\eta _0(\tau )+ϵ\eta _1(\tau )+\mathrm{}.`$ (94)
Inserting Eqs. (C) into Eqs. (92) and (93) and collecting powers of $`ϵ`$ we get
$`ϵ^0`$:

$`ds_0(\tau )`$ $`=`$ $`4\eta _0^2(\tau )d\tau `$ (95)
$`d\eta _0(\tau )`$ $`=`$ $`{\displaystyle \frac{30\nu _1}{30+\pi ^2}}\eta _0^3d\tau `$ (96)
$`ϵ^1`$:

$`ds_1(\tau )`$ $`=`$ $`8\eta _0(\tau )\eta _1(\tau )d\tau +{\displaystyle \frac{5\sqrt{3}}{4\sqrt{7}}}\sqrt{{\displaystyle \frac{D_1}{\eta _0^3(\tau )}}}dW_1(\tau )`$ (97)
$`d\eta _1(\tau )`$ $`=`$ $`{\displaystyle \frac{90\nu _1}{30+\pi ^2}}\eta _0^2(\tau )\eta _1(\tau )d\tau +`$ (98)
$`{\displaystyle \frac{15\sqrt{21+\pi ^2}}{2\sqrt{7}(30+\pi ^2)}}\sqrt{D_1\eta _0(\tau )}dW_2(\tau ).`$
Solving Eqs. (95) and (96) we obtain

$`s_0(\tau )`$ $`=`$ $`4{\displaystyle \frac{\eta _0^2(0)}{\lambda }}\mathrm{log}(1+\lambda \tau )`$ (99)
$`\eta _0(\tau )`$ $`=`$ $`{\displaystyle \frac{\eta _0(0)}{\sqrt{1+\lambda \tau }}}`$ (100)
with

$`\lambda ={\displaystyle \frac{60\nu _1\eta _0^2(0)}{30+\pi ^2}}.`$ (101)
Inserting Eq. (100) in (98) and solving, with the initial condition $`\eta _1(0)=0`$, we find

$`\eta _1(\tau )={\displaystyle \frac{45D_1(231+8\pi ^2)((1+\lambda \tau )^{5/2}1)}{56(30+\pi ^2)^2\lambda (1+\lambda \tau )^{3/2}}}+`$
$`{\displaystyle \frac{15\sqrt{21+\pi ^2}}{2\sqrt{7}(30+\pi ^2)}}\sqrt{D_1\eta _0(0)}{\displaystyle \underset{0}{\overset{\tau }{}}}(1+\lambda \tau ^{})^{5/4}𝑑W_2(\tau ^{}).`$
Then inserting Eqs. (C) and (100) in Eq (97) and integrating once we get

$`s_1(\tau )={\displaystyle \frac{15D_1(231+8\pi ^2)\eta _0(0)}{7(30+\pi ^2)^2\lambda ^2(1+\lambda \tau )}}(2(1+\lambda \tau )^{5/2}`$
$`5\lambda \tau 2)+{\displaystyle \frac{5\sqrt{3}\sqrt{D_1}}{4\sqrt{7}\eta _0^{3/2}(0)}}{\displaystyle }_0^\tau (1+\lambda \tau ^{})^{3/4}dW_1(\tau ^{})+`$
$`{\displaystyle \frac{60\sqrt{21+\pi ^2}}{\sqrt{7}(30+\pi ^2)}}\sqrt{D_1}\eta _0^{3/2}(0)\times `$
$`{\displaystyle \underset{0}{\overset{\tau }{}}}𝑑\tau ^{}{\displaystyle \frac{1}{\sqrt{1+\lambda \tau ^{}}}}{\displaystyle \underset{0}{\overset{\tau ^{}}{}}}(1+\lambda \tau ^{\prime \prime })^{5/4}𝑑W_2(\tau ^{\prime \prime }).`$ (103)
To this order we have $`Var(S(\tau ))=ϵ^2Var(s_1(\tau ))`$ and thus we finally obtain

$`Var(S(\tau ))`$ $`=`$ $`ϵ^2\left({\displaystyle \frac{75D_1}{112\eta _0^3(0)}}\tau +{\displaystyle \frac{225D_1\lambda }{448\eta _0^3(0)}}\tau ^2+O(\tau ^3)\right).`$
The full expressions of $`Var(S)`$ to this order of perturbation is given by Eqs. (III). Notice that the variance in the rest frame reads

$$Var(x(t))=\frac{1}{\alpha ^2}Var(S(\beta t)).$$
(105)
Concerning soliton velocity, up to first order perturbation it reads

$$4\eta ^2(\tau )=4\eta _0^2(\tau )+8ϵ\eta _0(\tau )\eta _1(\tau ).$$
(106)
Then, by substituting Eqs. (100) and (C) in Eq. (106), it is straightforward to see that

$`Var(4\eta ^2(\tau ))=64ϵ^2\eta _0(\tau )^2Var(\eta _1^2(\tau ))`$
$`=ϵ^2({\displaystyle \frac{3600(21+\pi ^2)D_1}{7(30+\pi ^2)^2}}\eta _0^3(0)\tau `$
$`{\displaystyle \frac{9900(21+\pi ^2)D_1\lambda }{7(30+\pi ^2)^2}}\eta _0^3(0)\tau ^2+O(\tau ^3)).`$ (107)
The full expressions of $`Var(4\eta ^2)`$ to this order of perturbation is given by Eqs. (III). In the rest frame this variance reads

$$Var(v(t))=\left(\frac{p^2}{12c}\right)^2Var(4\eta ^2(\beta t)).$$
(108)
In Fig. 6 we show some examples of $`Var(x)`$ and $`Var(v)`$ compared with results from a numerical solution of Eqs. (37) and (38) for which we have used the Heun method Heun . The variances have been obtained by averaging over 1000 realizations.

## Appendix D Thermalization process

From the generalized equipartition theorem Huang  we have that

$$\underset{n=1}{\overset{N}{}}Y_n(t)\frac{H(t)}{Y_n(t)}=Nk_BT$$
(109)
when the system is in thermal equilibrium with an external bath at temperature $`T`$. The relation (109) is strictly satisfied in the harmonic limit of the Hamiltonian (II), namely $`A=0`$. For finite values of $`A`$ the relation (109) is a rather good approximation to evaluate the temperature of our system when the relative displacements are sufficiently small, namely

$$Y_{n+1}(t)Y_n(t)<<\frac{3}{2A}.$$
(110)
Notice that $`A=1`$ in our simulations. The condition (110) can be obtained by comparing the harmonic term with the cubic anharmonicity in (II) and it is always satisfied in the present work.

So, in order to examine the thermalization process, one can define the ensemble average

$$H_{eqp}(t)=\underset{n=1}{\overset{N}{}}Y_n(t)\frac{H(t)}{Y_n(t)}.$$
(111)
For finite times $`H_{eqp}(t)`$ possesses two contributions, one due to the coupling to the thermal bath, $`H_{eqp}^{th}(t)`$, and the other one due to the soliton, $`H_{eqp}^{sol}(t)`$. So

$$H_{eqp}(t)=H_{eqp}^{th}(t)+H_{eqp}^{sol}(t),$$
(112)
Notice that for a very large system ($`N>>1500`$) the contribution of the soliton energy, $`H_{eqp}^{sol}(t)`$, can be neglected. However, in our system of 1500 sites this contribution is appreciable. In fact, the time evolution of $`H_{eqp}(t)`$ for different initial soliton energies presents different values due to the soliton contribution. Here, we remark that in thermal equilibrium, i.e. $`t\mathrm{}`$, these differences vanish. In Eq. (112) $`H_{eqp}^{th}(t)`$ gives us information about the time evolution of the temperature of the system in terms of the temperature of the thermal bath. The soliton contribution $`H_{eqp}^{sol}(t)`$ can be evaluated numerically using Eq. (111) when the soliton propagates in the presence of the damping but without noise. Then one can perform the numerical subtraction $`H_{eqp}(t)H_{eqp}^{sol}(t)`$ to obtain $`H_{eqp}^{th}(t)`$, which is shown in a normalized form in Fig. 7. We remark that we observe the same result here in systems with the double number of sites, namely 3000. The normalized thermalization process depends only on the damping constant which has the same value $`\nu =0.003`$ in all our simulations. Notice that for times $`t{}_{}{}^{<}\mathrm{\hspace{0.17em}3}/\nu =1000`$ there is a fast relaxation process, but for larger times, $`t{}_{}{}^{>}\mathrm{\hspace{0.17em}3}/\nu `$, the energy system approaches very slowly its stationary value. The thermal equilibrium of the system with the external bath corresponds to the case $`H_{eqp}^{th}(t)/k_BT=1`$.

## Appendix E Profiles

In Fig. 8 we compare snapshots of the system at $`t=5000`$ with and without noise and in the presence of damping. Figs. 8a and 8b correspond to both the kink shape (absolute displacements) and the pulse shape (relative displacements), respectively. Both shapes, with and without noise, in Fig. 8b were reconstructed from the shapes in Fig. 8a, respectively, by using the algorithm (45).


format: text

text:
# Thixotropy in macroscopic suspensions of spheres

## I Introduction

The rheology of suspensions has been studied for many years, mainly because of its obvious importance in a wide range of industrial applications (see for example , ).

A lot of work has been done especially for colloidal suspensions, i.e. suspensions where Brownian motion cannot be neglected , , , , , , , , , , , , , . In these colloidal suspensions shear-thinning can be observed, i.e. the viscosity of the fluid decreases with increasing shear rate , , , as well as shear-thickening for higher shear rates. It is known, that colloidal suspensions also show thixotropic behavior , , i.e. the continuous reduction in viscosity with time of shearing and the subsequent recovery under reduced shear , , , , . The idea is proposed, that the change in viscosity is caused by the formation of structures within a colloidal suspension under shear , , , . It is assumed that for colloidal suspensions the dynamics is dominated by Brownian diffusion , , .

A measure for the relative importance of shear forces compared to Brownian motion is given by the Péclet-number $`Pe=6\pi \eta \dot{\gamma }r^3/(k_BT)`$, where $`r`$ is the radius of the particles. Whereas colloidal suspensions have been investigated intensively up to a Péclet-number of 100, few experimental data exist for macroscopic suspensions, i.e. for suspensions of higher Péclet-number where Brownian motion can be neglected. It is assumed that for high Péclet numbers the shear-thinning region shifts rapidly to extremely low values and essentially becomes experimentally inaccessible , . The lack of Brownian motion for these large particles is believed to eliminate the reversibility in structural changes and should complicate experiments , . These difficulties explain the dearth of experimental data on the behavior of macroscopic, i.e. non-Brownian suspensions.

For this class interesting effects have been observed for non-Newtonian carrier fluids (see e.g. , , , ). Suspensions in Newtonian media appear to have been almost neglected. Some works showed the difficulties arising when trying to investigate the behavior of concentrated suspensions. It was demonstrated that shearing the macroscopic suspensions in a concentric cylinder-viscometer of the Searle-type leads to shear migration effects, i.e. the particles migrate in nonhomogeneous shear flows from high shear rates to low shear rates , ,  . For the geometry used in those experiments it meant that the particles even migrated out of the sheared volume and the particle concentration decreased. Since for high packing fractions the viscosity depends strongly on the concentration, this reduction in concentration was measured. So from these measurements no estimate of the real viscosity of the suspension could be made.

In another experiment  a macroscopic suspension was exposed to oscillatory shear and the decrease of the viscosity with time was measured. In this experiment there was a hint for an increase in the organization of the particles with shearing time. Nevertheless the suspension was not density matched, so it is not proved that the decrease in viscosity is not, at least partially, due to sedimentation.

The present paper presents experimental results on the behavior of macroscopic suspensions under steady shear at $`\dot{\gamma }=50/`$s of a Péclet-number $`Pe2.110^9`$ when using the viscosity of the suspension or $`Pe8.310^7`$ when using the viscosity of the carrier fluid, i.e. the Péclet number is seven orders of magnitude larger than the Péclet numbers which are covered by the experiments with colloidal suspensions. We find evidence for a thixotropic behavior of the suspension with extremely large relaxation times. The relaxation times are about 10000 times larger than the viscous time scale $`\tau _{visc}=a^2/\nu `$ and about 1000 times smaller than expected from the Péclet number $`\tau _{therm}=Pe/\dot{\gamma }`$. This makes it unlikely that a hydrodynamic or thermodynamic interaction is responsible for the effect. One can thus speculate that the time scale is dominated by ordering effects within the densely packed granular network provided by the beads . The decrease of the viscosity under shear coincides with an increase in order of the particles.

The paper is organized as follows: In the following section the experimental setup is described, in Section III the experimental results are discussed. A final section contains the conclusions.

## II Experimental setup

The suspension investigated consisted of nearly monodisperse glass beads (Würth Ballotini MGL) dispersed in an aqueous solution of zinc iodide $`ZnI_2`$. The glass beads are spherical and consist mainly of $`SiO_2`$ and have a material density of $`\rho =2.45`$ g/cm3. The particles are sieved in order to obtain a well defined particle size distribution which is then measured with a Multisizer 2 (Coulter). The particle size distribution (PSD) is shown in Fig. 1. This measurement of the PSD is based on the principle of monitoring the electrical current between two electrodes immersed in a conductive liquid on either side of a small aperture, through which the suspension of particles is forced to flow. As a particle passes through the aperture, it changes the impedance between the electrodes and produces an electrical pulse of short duration with an amplitude essentially proportional to the particle volume. The arithmetic mean of the particle size is 61.13 $`\mu `$m and the standard deviation is 18.71 $`\mu `$m. An alternative measurement of the particle size distribution was based on laser light diffraction experiments. These measurements confirm the results in Fig. 1. Before the glass beads were suspended in the solution they were washed 8 times in totally demineralized water, in order to purify them from any residues and dried afterwards.

To stabilize the suspension, i.e. to avoid sedimentation during the time-consuming rheological experiments, the density of the continuous phase has been matched to that of the glass beads. Hence the particles were suspended in almost saturated $`ZnI_2`$-solution, which had the same density $`\rho =2.45`$ g/cm3 at $`T=20.0^{}`$C as the glass beads. (The maximal density of the solution is about 2.6 g/cm3.) The density-matched, i.e. stabilized, suspension was prepared with the packing fraction $`\varphi =50.49`$%. It was checked at $`T=20.0^{}`$C for several hours that the suspension stayed homogeneous and the glass beads did not sediment to the bottom or rise to the top of the vessel. The carrier fluid, the pure $`ZnI_2`$-solution, is a Newtonian solution with a dynamic viscosity of 12 mPas at 20.0C.

The rheological measurements were performed with a rheometer SR-5000 (Rheometric Scientific, Inc.) with a rotating inner cylinder in a fixed outer cylinder. The inner cylinder has a cone at the lower end. Due to the relatively large particles a concentric cylinder geometry with a ratio $`\delta =R_i/R_o=0.922`$ of the inner cylinder radius $`R_i=14.75`$ mm to the outer cylinder radius $`R_o=16.0`$ mm and a gap width of $`a`$=1.25 mm has been used. The length of the inner cylinder is 44.3 mm. The shear stress in the small gap is constant within 15% for a Newtonian fluid. The thixotropic behavior was investigated in the controlled shear rate mode by a stepwise increase and subsequently stepwise decrease of the shear rate $`\dot{\gamma }`$ and registration of the shear stress $`\sigma `$ and the apparent viscosity $`\eta _{app}=\sigma /\dot{\gamma }`$. The temperature during the whole measurements was stabilized to 20.0 C.

The shear rates are below $`\dot{\gamma }=500/s`$, which corresponds to a Taylor number  $`Ta=\mathrm{\Omega }^2Ra^3\nu ^211.3`$, where $`\mathrm{\Omega }`$ is the angular velocity and $`R=(R_i+R_o)/2`$. Hence the experiments are performed far away from the onset of Taylor-vortices at $`Ta`$=1712.

In order to allow for an optical investigation of the structure-formation at the outer cylinder of the suspension, a shearing-apparatus was built which had the same gap width $`a`$=1.25 mm and almost the same ratio $`\delta `$=0.916 as the rheometer. The outer cylinder had a radius of 15.0 mm and was made out of glass to allow the view into the suspension layer at the outer cylinder. The inner cylinder did not have a conic but a flat lower end. Images were taken during shearing. In addition the torque of the motor was monitored, which is a measure for the viscosity of the suspension.

## III Experimental results

The following sequence of shear rate steps was used: $`\dot{\gamma }=`$ 0/s, 5/s, 50/s, 5/s, 0/s. The resulting shear stresses and the apparent viscosities $`\eta _{app}`$ were measured as a function of time. Fig. 2(a) shows the shear rate steps together with the system response, i.e. (b) the stress and (c) the dynamic viscosity $`\eta `$. When $`\dot{\gamma }`$ is set from 0 to 5/s an exponential decay of the viscosity to the final viscosity is observed. Then the shear rate is set to 50/s and another exponential decay to another apparent viscosity is observed. This final value of the viscosity is smaller than the final value at $`\dot{\gamma }=`$5/s. After setting $`\dot{\gamma }`$ back to 5/s the former final viscosity at this shear rate is gradually reached again. This is the behavior of a thixotropic material.

From this measurement the relaxation times $`\tau `$ are obtained from a fit to

$$\eta (t)=aexp(t/\tau )+\eta _{final}$$
(1)
Fig. 3 and Fig. 4 show the fits for $`\dot{\gamma }=`$50/s (from $`\dot{\gamma }=`$5/s) and for $`\dot{\gamma }=`$5/s (from $`\dot{\gamma }=`$50/s).

For the first shear rate $`\dot{\gamma }=`$5/s the relaxation time is $`\tau `$=1276 s. For the shear rate $`\dot{\gamma }=`$50/s the relaxation time is $`\tau `$=968 s and for the gradual rise in viscosity for $`\dot{\gamma }`$=5/s the relaxation time $`\tau `$=12182 s, i.e. it is one order of magnitude larger than for the decay in viscosity at the same shear rate.

The ramp for these shear rates was not applied in a repeated sequence with the same suspension, but was applied four times for freshly prepared suspensions and the thixotropic behavior of the suspensions was reproducible. For the following shear rates the ramp was applied once with alternating jumps from higher to lower and lower to higher shear rates.

This evaluation of the relaxation times was also done for other measurements for both possible cases: coming to the appropriate shear rate from a lower shear rate (step from $`\dot{\gamma }`$=0 to 5/s, from 5/s to 50/s and from 10/s to 100/s) or coming from a higher shear rate (step from $`\dot{\gamma }`$=100/s to 50/s and 500/s to 100/s). The duration of shearing was 3 hours when coming from a lower shear rate and 6 hours when coming from a higher shear rate, respectively. The relaxation times show a systematic decrease for increasing shear rates. For decreasing steps from higher shear rates to lower shear rates the relaxation times are larger than for increasing steps. All results are shown in Fig. 5.

Using the optical setup, the torque of the motor was monitored. Fig. 6 shows the decrease of the torque averaged over 100 s. An exponential fit according to equ. 1 gives the relaxation time of 1089 s, which is close to the one measured with the viscometer. Although the signal to noise ratio is worse than with the viscometer this allows the interpretation that the same mechanism is responsible for the torque reduction in both apparatuses.

Images with a Long Range Microscope (Questar QM100) which is connected to a CCD-camera using an exposure time of about 1 ms are taken from the suspension at the outer cylinder every 10 s during shearing at $`\dot{\gamma }=`$50/s. The starting time of the shearing is defined by $`t=0`$. The images show a cantle of 1.53 x 1.53 mm and are enhanced in contrast. Fig. 7(a) shows images of the suspension at the outer cylinder at $`t=0`$ and at $`t=62`$ min. There is almost no difference observable in real space but the 2D-Fourier spectra of the images (Fig. 7(b)) show a hexagonal structure after 62 min of shearing.

The increase of this hexagonal structure is quantified by calculating a measure for the contrast $`\alpha `$ as the root-mean-square normalized by the mean intensity of the spectral intensities for wave numbers $`k=\sqrt{k_x^2+k_y^2}`$ within a range of $`94/`$mm $`k106/`$mm, i.e.

$$\alpha =\frac{1}{\overline{i}}\sqrt{\frac{(\overline{i}i(k))^2}{n}},$$
(2)
where $`i`$ is the intensity in the spectrum and $`n`$ is the number of data points. This contrast is averaged over 20 images in order to reduce the noise. The unstructured spectra have a lower contrast $`\alpha `$ than the structured spectra, i.e. $`\alpha `$ is a measure for the structure of the suspension at the outer wall of the cylinder. Fig. 8 shows the temporal increase of that number.

## IV Conclusion

An experimental study of the viscosity of a macroscopic suspension of high packing fraction which is density-matched in a Newtonian carrier fluid under steady shear is presented. It is demonstrated that macroscopic suspensions can be considered as thixotropic materials.

Furthermore relaxation times at the order of hours were observed. They show a systematic decrease for increasing shear rates. For decreasing steps from higher shear rates to lower shear rates the relaxation times are larger than for increasing steps.

During the decrease of the viscosity under shear a structure-formation was observed. The Fourier spectra show that the system forms a hexagonal structure within time. Both effects, i.e. the decrease of viscosity and the increase of patterns, coincide.

We cannot provide an explanation for the physical origin of the thixotropic behavior. We note that the relaxation times are about 10000 times larger than the viscous time scale $`\tau _{visc}=a^2/\nu `$ and about 1000 times smaller than the thermodynamic time scale $`\tau _{therm}=Pe/\dot{\gamma }`$. This makes it unlikely that hydrodynamic or thermodynamic interactions are responsible for the effect.

One can thus speculate that the time scale is dominated by ordering effects within the densely packed granular network provided by the beads . In particular, we believe that the friction between the wetted surfaces of spheres is primarily responsible for the measured effective viscosity. In a densely packed granulate material that friction will be dominated by the fact that individual particles are subject to frustration if they feel more than one particle contact, which tends to increase the effective viscosity. Ordering under shear might mean that the particles reduce the mutual friction by reducing the number of frustrated contacts which in turn might cause a lower viscosity. This idea is related to the concept of dilatancy, which is fundamental in the description of fluid properties of granular matter. Additional investigations using three-dimensional imaging techniques and theoretical calculations are clearly necessary to substantiate these speculations.

## Acknowledgements

We like to thank Nuri Aksel for helpful discussions and Markus Horn for his help during the measurements. The experiment was supported by Deutsche Forschungsgemeinschaft through Re588/12.


text:
# Optical excitations of Si by time-dependent density-functional theory based on the exact-exchange Kohn-Sham band structure

## I Introduction

While the Kohn-Sham (KS) density-functional theory (DFT) Kohn99 ; Jones89 ; Koch00  has become the dominant method of first-principles investigations of energetic properties of solids, its scope in the treatment of electronic excitations has been rather limited due to deficiencies of the conventional local density approximation (LDA) and its semilocal extensions. One well-known case is the band gaps calculated within the LDA which are substantially smaller than experimental values Jones89 . Significant advancement has been made during the past decades in devising improved exchange-correlation functionals starting from the LDA Koch00 , but these semilocal approximations still fail to correct the so-called ‘band-gap problem’. Thus for the accurate treatment of band structures one at present usually leaves the domain of DFT and employs a many-body method based on the GW approximation (GWA) to quasiparticle energies Hybertsen85 ; Aulbur00 .

For the case of optical excitations in which bound exciton states can play a significant role, it is further required to go beyond effective one-particle approaches. In principle, this can be done within the time-dependent (TD) DFT Gross96 ; Casida95  or by solving the Bethe-Salpeter equation (BSE) for two-body Green’s functions Hanke79 . However, the lack of a proper exchange-correlation functional again remains as the hindrance to the TDDFT route, and the GWA-BSE currently represents the state-of-the-art method for such purpose Albrecht98 ; Benedict98 ; Rohlfing98 .

In this regard recent numerical realization of the exact-exchange (EXX) KS method for solids Bylander95 ; Gorling96 ; Stadele97  and molecules Grabo97 ; KLI99 ; Gorling99 ; Ivanov99 ; Veseth01 ; DellaSalla01  provides an interesting opportunity: According to theoretical and numerical evidence concerning the physical meaningfulness of the KS potential and KS eigenvalues Gorling96b ; Umrigar98 ; Al-Sharif98 , realistic band structures for semiconductors and eigenvalue spectra of molecules obtained within the EXX suggests a promising basis for the study of electronic excited states within DFT. In this work, we consider this point further by calculating the absorption spectrum of bulk silicon within TDDFT based on the EXX band structure and the adiabatic LDA (ALDA) exchange-correlation kernel. Silicon here serves as a representative semiconductor system, whose first absorption peak, the $`E_1`$ peak, is known to be strongly excitonic in its character. The large excitonic contributions to the $`E_1`$ peak has been proposed theoretically and experimentally, and a semiempirical calculation confirming this point has been performed two decades ago Hanke79 . However, truly first-principles calculations appeared only recently based on the GWA-BSE Albrecht98 ; Benedict98 ; Rohlfing98 . A different approach to investigate optical spectra within a DFT framework was recently explored in Ref. deBoeij01, , where the domain of conventional DFT was left and time-dependent current DFT was invoked. As a first step towards an accurate computation of optical spectra within TDDFT, we will examine the performance of the EXX-TDLDA approach in this paper.

## II Theory

To obtain the absorption spectrum of solids, we compute the imaginary part of the frequency-dependent macroscopic dielectric function $`ϵ(\omega )=ϵ_1(\omega )+iϵ_2(\omega )`$ defined as Pick70 ; Hybertsen87

$$ϵ(\omega )=[\underset{q0}{lim}ϵ^1(𝐪+𝐆,𝐪+𝐆^{};\omega )|_{𝐆=𝐆^{}=0}]^1$$
(1)
where

$$\begin{array}{cc}\hfill ϵ^1(𝐪+𝐆,𝐪+𝐆^{};\omega )& =\delta _{𝐆,𝐆^{}}+F_H(𝐪+𝐆)\hfill \\ & \times \chi (𝐪+𝐆,𝐪+𝐆^{};\omega ),\hfill \end{array}$$
(2)
with $`F_H(𝐪+𝐆)4\pi /|𝐪+𝐆|^2`$ (Hartree atomic units are used throughout the paper. The symbols $`𝐪+𝐆`$ refer to plane waves with wave vectors $`𝐪+𝐆`$, $`e^{i(𝐪+𝐆)𝐫}`$, i.e., all formulas are given in the plane wave representation). The linear response function $`\chi `$ that describes the full response of the first-order density change $`\delta n`$ for the given external perturbation $`\delta v_{ext}`$,

$$\delta n(𝐪+𝐆;\omega )=\underset{𝐆^{}}{}\chi (𝐪+𝐆,𝐪+𝐆^{};\omega )\delta v_{ext}(𝐪+𝐆^{};\omega ),$$
(3)
can be computed either via the BSE or TDDFT. Within TDDFT the linear response function $`\chi `$ is given in matrix representation as Gross96 ; Casida95

$$\chi (𝐪;\omega )=\left[1\chi _0(𝐪;\omega )\{F_H(𝐪)+F_{xc}(𝐪;\omega )\}\right]^1\chi _0(𝐪;\omega ).$$
(4)
where $`\chi _0`$ is the noninteracting KS linear response matrix. It is given in terms of the occupied and unoccupied KS orbitals $`\varphi _{v𝐤}`$ and $`\varphi _{c𝐤}`$ and their eigenvalues $`ϵ_{v𝐤}`$ and $`ϵ_{c𝐤}`$ through

$$\begin{array}{cc}& \chi _0(𝐪+𝐆,𝐪+𝐆^{};\omega )=\frac{2}{N\mathrm{\Omega }}\hfill \\ & \times \underset{v,c,𝐤}{}[\frac{v,𝐤|e^{i(𝐪+𝐆)𝐫}|c,𝐤+𝐪c,𝐤+𝐪|e^{i(𝐪+𝐆^{})𝐫^{}}|v,𝐤}{ϵ_{v𝐤}ϵ_{c(𝐤+𝐪)}+\omega +i\delta }\hfill \\ & +\frac{c,𝐤|e^{i(𝐪+𝐆)𝐫}|v,𝐤+𝐪v,𝐤+𝐪|e^{i(𝐪+𝐆^{})𝐫^{}}|c,𝐤}{ϵ_{v𝐤}ϵ_{c(𝐤+𝐪)}\omega i\delta }],\hfill \end{array}$$
(5)
where $`N\mathrm{\Omega }`$ is the crystal volume. In Eq. (4) $`F_{xc}`$ denotes the exchange-correlation kernel, the frequency-dependent functional derivative of the exchange-correlation potential.

Thus, the macroscopic dielectric function $`ϵ`$ can be computed by first calculating $`ϵ^1`$ in the reciprocal-space using Eqs. (4) and (5), and next taking the inverse of the $`q0`$ limit of its $`𝐆=𝐆^{}=0`$ element \[Eq. (1)\]. However, the ALDA kernel $`F_{xc}^{ALDA}=\delta ^2E_{xc}^{LDA}/(\delta n\delta n^{})`$ is local in real space and frequency-independent, which results in a reciprocal-space representation for $`F_{xc}`$ that is independent of $`q`$ and $`\omega `$, $`F_{xc}^{ALDA}(𝐪+𝐆,𝐪+𝐆^{};\omega )=F_{xc}^{ALDA}(𝐆𝐆^{})`$. This incorrect $`q0`$ behavior of the ALDA kernel then simplifies the calculation of $`ϵ(\omega )`$: Writing $`\chi `$, $`F_H`$, and $`F_{xc}`$ (from now on the superscript “ALDA” on the kernel is suppressed for notational simplicity) for $`q0`$ in block form in terms of the “head” \[$`𝐆=𝐆^{}=0`$ (00)\], “wings” \[$`𝐆=0;𝐆^{}0`$ (01) and $`𝐆0;𝐆^{}=0`$ (10)\], and the “body” \[$`𝐆0;𝐆^{}0`$ (11)\] as

$$\begin{array}{cc}\hfill \chi _0(q;\omega )& =Q\left(\begin{array}{cc}\chi _0^{00}(\omega )& \chi _0^{01}(\omega )\\ \chi _0^{10}(\omega )& \chi _0^{11}(\omega )\end{array}\right)Q;\hfill \\ \hfill F_H(q)& =Q^1\left(\begin{array}{cc}4\pi & 0\\ 0& F_H^{11}\end{array}\right)Q^1;\hfill \\ \hfill F_{xc}(q)& =Q^1\left(\begin{array}{cc}0& 0\\ 0& F_{xc}^{11}\end{array}\right)Q^1\hfill \end{array}$$
(6)
where $`Q^{00}=q;Q^{01}=Q^{10}=0;Q^{11}=\delta _{𝐆,𝐆^{}}`$, and performing matrix multiplications explicitly (with $`𝐪`$-singularities treated properly at each stage) yields a $`3\times 3`$ tensor Pick70 ; Hybertsen87 ,

$$\begin{array}{cc}\hfill ϵ(\omega )=14\pi & \chi _0^{00}(\omega )4\pi \chi _0^{01}(\omega )(F_H^{11}+F_{xc}^{11})\hfill \\ \hfill \times & [1\chi _0^{11}(\omega )(F_H^{11}+F_{xc}^{11})]^1\chi _0^{10}(\omega ),\hfill \end{array}$$
(7)
with which matrix multiplications including head and wings with explicit $`𝐪`$-dependence are avoided. Eq. (7) simply represents the independent-particle-level absorption spectrum if $`F_H^{11}`$ and $`F_{xc}^{11}`$ are ignored altogether, while it is reverted to the so-called RPA dielectric function formula if only $`F_{xc}^{11}`$ is neglected. For cubic systems, such as Si considered here, the macroscopic dielectric tensor is diagonal.

## III Computational methods

For the calculation of the dielectric function, one needs wavefunctions $`\varphi _{i𝐤}`$ and eigenvalues $`ϵ_{i𝐤}`$ at the chosen set of $`𝐤`$-points \[See Eq. (5)\]. To generate these, we first performed the LDA and EXX self-consistent ground-state calculations at the experimental lattice constant, $`a=5.43\AA `$ for Si, and obtained the KS potential $`v_{KS}`$ which serves as the only input to the TDDFT calculation of optical spectra. The EXX calculation was carried out using ten special-$`𝐤`$ points, orbital kinetic energy cutoffs of 12.5 Ha, and nonlocal pseudopotentials in the separable form of Kleinman-Bylander Kleinman82 . We refer the reader to Ref. Stadele97,  for further details of EXX calculations. In addition, because the ALDA kernel $`F_{xc}^{ALDA}`$ depends only on the density, we constructed it at the end of the ground-state calculations when fully converged densities were available.

With the converged $`v_{KS}^{EXX}`$, $`v_{KS}^{LDA}`$ and $`F_{xc}^{ALDA}`$ in hand, we moved on to the computation of the dielectric function and first solved the KS equations once again at a much larger number of $`𝐤`$-points but with the same energy cutoff 12.5 Ha to generate $`\chi _0`$. The calculation of the head and wings of $`\chi _0`$ needs special care since the $`q0`$ limit of $`i,𝐤|e^{i𝐪𝐫}|j,𝐤+𝐪`$ should be handled properly. Applying perturbation theory one obtains the following equation Hybertsen87 ,

$$\underset{q0}{lim}i,𝐤|e^{i𝐪𝐫}|j,𝐤+𝐪=𝐪\frac{i,𝐤|𝐩+[V_{NL},i𝐫]|j,𝐤}{ϵ_{j,𝐤}ϵ_{i,𝐤}},$$
(8)
which implies that, due to the extra powers of the energy denominator that appear in the head $`\chi _0^{00}(\widehat{𝐪};\omega )`$, computation parameters for the convergence are more or less determined by those necessary for the head element. We use the kinetic energy cutoff 7 Ha and 10 conduction bands for $`\chi _0`$ which was sufficient to reach the convergence of $`ϵ_2(\omega )`$.

For the $`𝐤`$-point sampling, we chose to adopt a regular uniform grid, with which one could evaluate $`lim_{q0}i,𝐤|e^{i𝐪𝐫}|j,𝐤+𝐪`$ by a finite difference approach without resorting to the analytic replacement of Eq. (8). In this case one had to pay the price of performing three additional diagonalizations at $`𝐤+q𝐞_\alpha (\alpha =x,y,z)`$ for each $`𝐤`$-point. Therefore, in this paper we report results obtained by treating the $`q0`$ limit according to Eq. (8) with ignoring the nonlocal pseudopotential contributions due to the commutator $`[\widehat{V}_{NL},i𝐫]`$ term. In test calculations with a small number of $`𝐤`$-points we found that the neglect of the nonlocal pseudopotential contributions has little effect on the positions of absorption peaks while the peak heights are somewhat reduced. We slightly shift the uniform $`𝐤`$-mesh in order to break any present symmetry thus achieve a more effect $`𝐤`$-point sampling Rohlfing98 . After checking the convergence behavior by increasing the number of $`𝐤`$-points and changing the shifting vector, we used a shifted $`17\times 17\times 17`$ regular grid which gives the convergence of $`ϵ_2`$ up to the frequency range of at least 6 eV. The use of a shifted uniform $`𝐤`$-mesh results in slightly non-identical diagonal elements of the dielectric tensor whose differences have been also used to check the convergence with respect to the number of $`𝐤`$-points.

## IV Results and discussions

We first consider the absorption spectrum calculated at the LDA and EXX independent particle levels. In Fig. 1, LDA and EXX absorption spectra are shown together with experimental Herzinger98  and GWA Albrecht98 ; Benedict98  curves. Note that the two GWA results are rather different due to operational ambiguities resulting from the post-DFT nature of the GWA Aulbur00 , while the LDA and EXX which strictly reside in the KS DFT scheme do not suffer from such problems. Within the LDA, the absorption spectra is shifted to the lower frequency region by about 1 eV due to its well-known underestimation of bands gaps. In addition, the $`E_1`$ peak height is significantly underestimated, while that of the $`E_2`$ peak is much overestimated. In contrast, the EXX absorption spectrum shows a much better overall agreement reflecting realistic descriptions of the band structure within the EXX Stadele97 : Positions of the absorption edge and the $`E_2`$ peak are in excellent agreement with the experiment, and the $`E_2`$ peak position is even slightly better than the GWA one. Also the height of the $`E_2`$ peak is smaller than the LDA value and is much closer to the experimental data. One noticeable disagreement is the absence of the $`E_1`$ peak as in the LDA spectrum. This, however, does not reflect the deficiency of the EXX approach, but the excitonic nature of the $`E_1`$ peak Albrecht98 ; Benedict98 ; Rohlfing98  which cannot be described at the independent particle level as can be also seen in the GWA curves.

We feel that the excellent features of the EXX absorption spectrum warrant more discussions. Perturbation theory along the adiabatic connection shows that KS eigenvalue differences are well-defined approximations of excitation energies Gorling96b . Indeed Umrigar et al. recently found for several atoms that quasi-exact KS eigenvalue differences are in surprisingly good agreement with experimental excitations Umrigar98 . In addition, applying the quasi-exact KS potential to the calculation of the quantum defects for the Ne atom, Al-Sharif et al. argued that the KS potential of a neutral atom (with $`N`$ electrons) corresponds to the Green-function “potential” of its daughter cation ($`N1`$ electrons), or equivalently the KS equation for $`N`$ electrons corresponds to the Dyson equation where the reference ground state is chosen with $`N1`$ electrons Al-Sharif98 . The appealing description of excitations by KS eigenvalue differences, in contrast to the treatment of ionizations and electron affinities in the many-body (GWA) approach, has not drawn much attention due to the lack of realistic exchange-correlation approximations which are suitable for practical calculations. Exchange-correlation potentials derived from the LDA and other semilocal exchange-correlation energy functionals do not cancel Coulomb self-interactions, decay exponentially rather than as $`1/r`$ for localized systems, and inherently lack the term that results from the functional derivative of the pair-correlation function Gritsenko94 . The self-interaction-free EXX method corrects such problems and provides accurate KS potentials and eigenvalues. Then, we propose that the good absorption spectrum from such EXX method represents further evidence for the $`N`$-electron excitations picture within the KS scheme.

Now we consider the effects of the inclusion of two-particle interactions via TDLDA on top of the independent-particle EXX result. The Hartree-only or RPA (EXX+TDH) and full TDLDA (EXX+TDH+ALDA) spectra are shown together with the EXX independent-particle spectrum in Fig. 2. Granted the overall good feature of the EXX spectrum, the most interesting question is that whether the ALDA can generate the excitonic $`E_1`$ peak. Fig. 2 shows that this is not the case. TDH alone suppresses the EXX spectrum and shifts the $`E_2`$ peak to a slightly higher frequency region, while adding the ALDA counteracts such effects and more or less restores the EXX independent-particle spectrum. In Fig. 2 we also display the TDLDA absorption spectrum based on the LDA band structure (LDA+TDH+ALDA) to show that the LDA independent-particle result (see Fig. 1) is again insignificantly changed by TDLDA effects.

## V Conclusions

In this work, by calculating the frequency-dependent dielectric function of bulk Si we studied the quality of EXX one-particle band structure and the description of excitonic effects within the TDLDA. We showed that, at the EXX independent-particle level, the position of the $`E_2`$ peak is in good agreement with experimental data and is even slightly better than that obtained in the GWA approach. We argued that this is because the EXX approach, unlike procedures employing the LDA and semilocal density-functionals, yields realistic KS orbitals and eigenvalues which accurately model electronic excitations at a fixed particle number. Employing TDDFT within the ALDA on top of the EXX band structure however results in only a minor modification of the EXX spectrum and fails to produce the excitonic $`E_1`$ peak. Given the fact that proper inclusion of excitonic effects is mandatory for the complete description of optical excitations, devising a TDDFT scheme that overcomes the deficiency of the ALDA kernel will be an interesting task. Work along this line using the nonadiabatic and nonlocal EXX kernel is in progress TDEXX .

###### Acknowledgements.

We thank Professor W. Domcke for his support including computational resources. This work was supported by the Deutsche Forschungsgemeinschaft and the Fonds der Chemischen Industrie. Y.-H. Kim also acknowledges receipt of a Humboldt Foundation award.