# Solitons and thermal fluctuations in strongly non-linear solids

###### Abstract

We study a chain of anharmonic springs with tunable power law interactions as a minimal model to explore the propagation of strongly non-linear solitary wave excitations in a background of thermal fluctuations. By treating the solitary waves as quasi-particles, we derive an effective Langevin equation and obtain their damping rate and thermal diffusion. These analytical findings compare favorably against numerical results from a Langevin dynamic simulation. In our chains composed of two sided non-linear springs, we report the existence of an expansion solitary wave (anti-soliton) in addition to the compressive solitary waves observed for non-cohesive macroscopic particles.

###### pacs:

45.70.-n, 61.43.Fs, 65.60.+a, 83.80.FgIn linear elastic solids, phonons are the basic mechanical excitations responsible for energy propagation. By contrast, an aggregate of macroscopic grains just in contact with their nearest neighbours constitute a novel elastic material where solitary waves or shocks replace phonons as the basic excitations Nesterenko_1984 (); Rosenau_1986 (); Gomez_2012 (). The origin of these strongly non-linear waves can be traced to the fact that, unlike the case of harmonic springs, the repulsive force between two grains in contact does not depend linearly on the relative compression. So far, little effort has been directed to determine the fate of these strongly non-linear excitations in a background of thermal fluctuations because temperature is clearly not a parameter relevant to the elastic response of macroscopic grains.

However, granular aggregates at zero pressure are just one example of a broader class of materials that can be prepared in a unique mechanical state called sonic vacuum Nesterenko_1984 (). This term originally coined by Nesterenko in the context of strongly non-linear granular chains designate a material characterized by a vanishing elastic moduli and linear speed of sound Nesterenko_1984 (); Nesterenko_Book (); Coste2 (); Daraio2005 (); Job (); Sen (). Grafted colloidal particles Anand_2007 () and ultra-cold atoms in optical lattices Chang_2013 () are microscopic systems that allow for tunable non-linear interactions, while being naturally coupled to a source of fluctuation (thermal or quantum). These fluctuations restore rigidity and generate long wavelength phonon modes Chaikin_2000 (); Zhirov_2011 (). However, the physics of very high amplitude strain propagation is still predominantly non-linear and resembles the state of sonic vacuum perturbed by background fluctuations, even if the interaction potentials are typically two sided, unlike granular ones. The non-linear regime is particularly relevant for some biological systems, where energy transport occurs through localized non-linear excitations with energy significantly higher than the thermal energy Davydov_1977 (); Peyrard_Book ().

Moreover, systems such as polymer networks and colloidal glasses undergoing an unjamming transition are also characterized by vanishing elastic moduli as the coordination number or packing fraction are lowered towards the critical point OHern (). The effect of thermal fluctuations on the non-linear response of materials undergoing an unjamming transition is relatively unexplored, despite they are obvious examples of a sonic vacuum state at zero temperature Mackintosh_2012 (); berthier (); Ning (). Note, that in the case of jamming the linear elastic moduli can be lowered towards zero even if the microscopic interactions are harmonic, simply because there are not enough forces to prevent floppy motions.

In this article, we focus on strongly non-linear mechanical waves propagating in a background of small thermal fluctuations, a non-equilibrium problem that lies outside the realm of perturbation theory. The starting point of conventional perturbation methods is a linear elastic solid, possibly at finite temperature, perturbed by small anharmonic terms. By contrast, we adopt as a starting point the fully non-linear state of sonic vacuum whose elementary excitations are long-lived solitary waves remark1 (). Subsequently we switch on temperature as a small perturbation that creates a background of thermal fluctuations.

As a minimal model that is analytically tractable, we study impulse propagation in a one dimensional lattice of non-linear springs with a tune-able power law interaction. By coupling the lattice to a heat bath, we then study the effects of the thermal fluctuations on the leading solitary wave generated in response to an impulse of energy much higher than the background thermal energy. Our approach in a nutshell is to treat the solitary wave as a quasi-particle and derive an effective Langevin equation that describes its stochastic dynamics. We corroborate our analytical predictions for the damping rate and thermal diffusion of the solitary waves with Langevin dynamic simulations. The sonic vacuum is usually studied with a chain of non-cohesive beads that only interact upon compression (one-sided interaction). The system with springs has some of the same properties as the sonic vacuum, since the sonic vacuum is a property of the non-linear power law interaction (without a harmonic term). In addition to the compressive solitary waves seen in a lattice of macroscopic grains with one-sided repulsive interaction, we report an accompanying anti-solitary wave solution for the lattice of non-linear springs with two sided interactions.

## I The impulse response of non-linear springs

In Fig. (1), we demonstrate that the compressional solitary wave (SW) excitation discovered by Nesterenko in a chain of non-cohesive beads is also seen in a lattice of springs with two sided interactions. However, unlike the case of a one sided potential, each compressive SW generated in response to an impulse is accompanied by a corresponding expansion solitary wave (formed by local stretching of springs) of the same magnitude but moving in the opposite direction. This anti solitary wave (ASW) is not sustained by beads interacting with purely repulsive potentials – the beads would merely loose contact.

In Fig. (1) left panel, we show the SW/ASW excitations for (a)beads , (b-c) particles connected by springs. In Fig. (1), right panel, we plot, , the total energy carried by the soliton (after summing over all particles envolved) versus , the total momentum, for the leading SW/ASW. Note that , where the constant can be viewed as the effective mass of the solitary or antisolitary waves. Inspection of Fig. (1) demonstrates that SW excitations in a lattice of repulsive beads (black circles) have the same effective mass as a SW and ASW in two-sided springs (red squares).

As shown in Fig. (1) left panel (b-c), the leading SW-ASW generated in response to an impulse imparted to one of the particles towards the right (direction of arrow) is followed by a train of alternating SW-ASW’s excitations, of progressively smaller magnitudes. The smaller SW/ASW’s are generated as the particle that initially imparted the impulse, recoils with its left-over energy. This process is repeated several times, leading to the generation of the train of smaller excitations. Since the speed of propagation depends upon the amplitude, the SW and ASW that start propagating together initially (appearing bounded), eventually separate and become clearly distinguishable.

## Ii Langevin Equation

The classical energy-momentum relation satisfied by the SW motivates the interpretation of the solitary wave as a quasi-particle Nesterenko_1984 (); Job (). For small perturbations, the SW can still be treated as a quasi-particle provided the effects of the perturbations accrue gradually such that the SW retains its functional form. We now apply this adiabatic approximation to derive an effective Langevin equation for the SW quasi-particle when the lattice of springs is coupled to a heat bath. Recall first, the Langevin equation for a particle of mass undergoing Brownian motion in one dimension is

(1) |

Here, are the total and kinetic energies respectively, are the dissipation and diffusion coefficients related via the fluctuation dissipation theorem , where is the Boltzmann constant. is a normal random variable with mean 0 and variance 1, and encapsulates the effects of random fluctuations during the time interval . For a free particle of unit mass moving with speed , and upon substituting in Eq. (1), we recover the Langevin’s equation conventionally expressed as the rate of change of momentum of the particle Kampen_Book ().

We now derive an equation analogous to Eq. (1) for the compressive solitary wave quasi-particle. Let the displacement of a particle from its initial equilibrium position in the continuum limit be . If we identify the lattice spacing as a characteristic length scale and as an inverse time scale, then the equation of motion for the compressive displacement field in dimensionless units reads

(2) |

where subscripts denote partial derivatives with respect to space and time . Eq. (2) is a simplified form of the Nesterenko equation Nesterenko_1984 (); Rosenau_1986 (), see appendix A for details. The first two terms express the rate of change of momentum while the third term represents the force. Although the solitary wave solution to Eq. (2) is not exact (lacking compact support), Eq. (2) provides a good approximation while being analytically more tractable especially since we are interested in keeping the non-linear exponent general Rosenau_1986 (); Gomez_2012 (). Note that the equation for the ASW (stretching) is obtained by modifying the third term in Eq. (2).

In analogy with the Langevin equation for a particle, we model the coupling to a heat bath as the sum of two contributions- an external drag and a random fluctuating force, phenomenologically introduced into the equation of motion as :

(3) |

where is the dimensionless drag coefficient that couples to the momentum . It is useful to define a coupling constant as the ratio of potential to thermal energy in terms of which, the dimensionless diffusion coefficient is . The last (noise) term on the right of Eq. (3) in conjunction with , satisfies the fluctuation dissipation theorem Bishop_2008 (). Here, is a Gaussian random noise during the time interval with the moments and respectively, where angular brackets denotes ensemble averaging.

To study the propagation of the SW in a background of thermal fluctuations, we now make a working assumption based on the quasi-particle approximation to the SW: whenever the energy of the SW, , the SW functional form is unaltered and only its amplitude becomes time-dependent. The amplitude is the collective variable for the SW quasi-particle and other properties of the solitary wave, such as its energy and momentum may be determined from it. Note, the width of the SW is independent of its amplitude and therefore we do not consider its time dependence Bishop_2008 ().

From Eq. (2), the conserved energy is

(4) |

and the energy of the SW may be obtained by integrating Eq. (4) over the width of the SW of order . (This avoids including the energy of small SW that separate from the main wave). Using Eq. (2), the rate of change of energy is,

(5) |

where, is the kinetic part of the energy,

(6) |

The last two terms on the right of Eq. (5) describe the possible mechanisms of decay of the SW by “friction” from the heat bath () and the ”phonon drag” induced by the thermal motion of the chain . The first term is the fluctuating part of the energy. In the following, we make the assumption (verified numerically) that the coupling to the heat bath is more important and therefore, ignore (see SI for details).

Solving for the SW solution from Eq. (2), we find , where

(7) |

is the functional form of the SW with amplitude , speed and width in units of the lattice spacing, see appendix A, subsection 1, for details. The SW energy, kinetic energy, and momentum may now be expressed in terms of the collective variable :

(8) |

and from the the virial theorem,

(9) |

Additionally, the solitary wave momentum is

(10) |

Here,

(11) | |||

(12) |

are constants obtained by integrating over all space Rosenau_1986 ().

Substituting for and in terms of , we cast Eq. (5) into the form of an ordinary Langevin equation (with additive noise) for the collective variable ,

(13) |

where, . Eq. (13) is the central result of our work whose analytical predictions we derive and test numerically in the next sections. The first term can be written as , where is a white noise signal; that is, its correlations are given by . Using the fact that the correlations of are described by delta-functions, the correlations of can be related to the kinetic energy Eq. (6), which can be replaced by .

## Iii Time dependence of mean and variance

Taking the expectation value (ensemble average) of Eq. (13), we find

(14) |

where, owing to the noise term (which acts between times ) and (which is a solution at time ) being statistically independent, the expectation value . Consequently, the solitary wave amplitude decays as

(15) |

where, is the initial solitary wave amplitude. Note, the effective damping rate is independent of inverse temperature but rescales with the exponent of the non-linear potential .

Similarly, we solve for the variance of the solitary wave amplitude or equivalently, the variance in the square root of energy. Re-defining, , we solve for the variance in the solitary wave amplitude by first evaluating the differential , by substituting from Eq. (7) and retaining terms to order Gillespie_Book (); Lemons_Book (),

(16) |

Taking the expectation value, the second term on the right vanishes (as discussed for the mean) and using the property that the noise term is delta-correlated in space, we obtain

The last term when expanded gives twice the solitary wave kinetic energy , see Eq. (6), plus an integral , that vanishes by symmetry for the SW solution. Moreover, the SW kinetic energy is related to its total energy via the virial relation . Hence, we obtain the ordinary differential equation correct to order -

(17) |

Solving, the differential equation subject to the initial condition and substituting for , we obtain,

(18) |

Using Eq. (15), this may be expressed as

(19) |

Using the relation in Eq. ( 8), we rewrite the above equation as

(20) |

The coefficient reduces to when the energy is not measured in units of , so this expression is analogous to the velocity variance of a Brownian particle. Note, for large , the SW is effectively one particle wide and thus Eq. (20) captures the correct thermal equilibration of the particle energy with the heat bath. However, for the dynamics of the SW, Eq. (20) is only useful as long as the SW is identifiable against the background thermal energy, that is .

## Iv Simulations

We consider a one dimensional chain consisting of particles each having a mass placed regularly on a lattice with spacing (spring rest length) interacting pair-wise with a nearest neighbour interaction , where is the compression/stretching induced during the dynamics. We model the coupling to the heat bath by numerically integrating Eq. ( 1) for each particle using the velocity-verlet algorithm Allen_Book (). In thermal equilibrium the mean kinetic energy is KE and potential energy is PE, where their ratio satisfies the virial relation, see Fig. (3), inset for . In the following, all numerical data is presented in dimensionless units, ensemble averaged over 1000 samples.

### iv.1 Fluctuation induced rigidity

To extract the equilibrium properties in the thermalized state, we define the longitudinal current density of particles as , and its Fourier transform , where is the longitudinal collective mode along the -direction. Thus, the corresponding longitudinal current density auto-correlation functions is where the angular brackets denote ensemble averaging over the initial time. The longitudinal power spectral density is then obtained as the Fourier transform of the respective current density auto-correlation functions as, . The Fourier transforms defined above are evaluated using fast Fourier transform from simulation data. The sound speeds in Fig. (3) correspond to the linear part of the dispersion curves, obtained by projecting the power spectral densities on the frequency ()- wavenumber () plane.

In Fig. 3, we plot the sound speed from the slope of the dispersion curves for for a range of . At thermal equilibrium, the mean kinetic energy and hence the temperature satisfy the virial relation , where is the average displacement of the particles induced by thermal fluctuations. Defining the sound speed as the second derivative of the induced potential energy leads to the relation, Zhirov_2011 (). For , we find , closely matching the linear fit in SI Fig. ( 1). Thus, coupling the lattice of non-linear springs that is initially in its state of sonic vacuum (implying the absence of linear sound) to a heat bath, leads to hydrodynamical sound modes with a linear sound speed that scales with the temperature of the heat bath Zhirov_2011 (). Note, setting (harmonic springs) yields a sound speed that is independent of temperature while the limit yields , a result in agreement with the entropic elasticity for hard sphere colloidal crystals Chaikin_2000 ().

### iv.2 Comparison with analytics

Once the lattice reaches thermal equilibrium, we excite a solitary wave (SW) by imparting one of the particles an initial energy of order in dimensionless units. In Fig. (2) left panel, we show a snapshot of two SWs at the same time, propagating in a background of thermal fluctuations for (red) and (black). We see that the SW with lower is wider and moves faster for the given amplitude, in qualitative agreement with the analytic widths and speeds .

In Fig. (4) left panel, we plot the numerical data (symbols) for the attenuation of the SW amplitude as a function of time for various values of and and we find a very good match to the analytic expression in Eq. (15) (solid curves). For the range of explored, we find the damping rate is independent of temperature () but depends on the environmental drag and .

In Fig. (4) right panel, we show the increase in the variance of SW amplitude (or the square root of its energy) as a function of time for multiple values of , and obtained numerically (symbols) and compare them with the complete analytical solution Eq.(20) finding good agreement. Notice, the final value of the variance correctly approaches the thermal energy, as expected for a Brownian particle. However, since the solitary wave is a dynamical object that decays under the influence of the external drag, once the solitary wave energy becomes comparable to the background thermal energy, it is no longer meaningful to consider it as a Brownian particle.

## Appendix A Continuum approximation

In this section, we will review the Rosenau approximation to the Nesterenko solitary wave solution, that is valid for any general non-linear potential Rosenau_1986 (); Gomez_2012 ().

Here, we adopt as our starting point the Lagrangian for a one dimensional chain of identical spheres that are just touching each other, i.e., in the limit -

(21) |

where, is the displacement of the th sphere from its equilibrium position, is the equilibirum lattice spacing and is the spring constant. In order to avoid doing a Binomial expansion in powers of , we will define the continuum field variable as

(22) |

where, primes denote derivative with respect to . We now take the continuum limit, i.e., and Taylor expand the right hand side about :

(23) |

Integrating both sides once with respect to , we obtain

(24) | |||||

(25) |

Inverting the differential operator, we obtain

(26) |

Thus, in the contiuum limit, the Lagrangian becomes

(27) | |||||

(28) |

By using the Euler-Lagrange equation, we obtain the equation of motion as

(29) |

Note, here corresponds to the continuum displacement field. The corresponding equation in the strain field reads

(30) |

Upon substituting for the compressive SW or for the expansive ASW, we find the same functional forms for the solitary wave solutions in both cases. Here, represents the compression of two adjacent particles i.e., the strain field.

### a.1 Solitary wave solution

The solitary wave solution of Eq. 30 can be obtained by looking for propagating solutions of the form :

(31) |

which can be expressed in the form of Newtons’s-like equation

(32) |

Multiplying both sides by and integrating,

(33) | |||||

(34) | |||||

(35) |

Integrating again, we find

(36) |

Substituting, , and writing for brevity and , we need to integrate

(37) |

Making the change of variables, , we find

(38) |

Therefore,

(39) |

that yields,

(40) |

or

(41) |

Squaring and substituting ,

(42) |

that yields,

(43) |

Therefore, the solitary wave solution is

(44) |

## Appendix B Solution from discrete equations of motion

The energy and the fluctuations of the solitary wave can also be derived for a discrete chain, without making the continuum approximation. From Eq. (21), the equation of motion is

(45) |

A solitary wave solution to the above equation of motion has the form of a wave moving at a constant speed :

(46) |

where is the amplitude of the SW and is a function that describes the shape of the SW (Since we define the amplitude as the maximum speed of the particles enveloped by the solitary wave, the displacement will be proportial to .)

In analogy with Eq. (1), we now couple the discrete equation of motion Eq. (45) to a source of Gaussian noise and drag :

(47) |

Here, the noise satisfies . The coefficient of the noise ensures that in equilibrium the particles satisfy the fluctuation-dissipation theorem such that their average kinetic energy is .

From Eq. (21), the energy of the chain is . If we restrict the summation to particles enveloped by the solitary wave such that the sum does not include particles that are far away from the SW, then corresponds to the energy of the SW. Multiplying Eq. (47) by , we obtain the rate of change of energy

(48) |

where, the contribution of terms of the form (i.e., with index and ) cancel out upon summation over .

Consistent with the quasi-particle interpretation of the SW, we again assume that the SW in a background of noise and drag, is still given approximately by Eq. (46). Therefore, we write the terms in Eq. (48) in terms of the SW amplitude. The SW energy is , and the first term on the right hand side is proportional to the kinetic energy, which is by the virial theorem. Thus the energy decays at a rate times the decay rate for the velocity of a single particle (). The SW amplitude therefore decays as:

(49) |

A term that compensates for the diffusion of is omitted here (this term is derived in Ito calculus), but it is only a small amount that is negligible if . Now, the different contributions in the last term just add up to a new noise term , which is also not correlated in time. (Assume that is the amplitude of the noise while has a noise with one unit of strength.) The original noise term has correlations because depends on the noise at an earlier time. But if the noise mostly just changes the amplitude of the SW by a random amount, and does not cause the particles in it to be displaced randomly relative to each other, then these correlations are not too important, and that is why they cancel out in the equation for . The variance of the noise is therefore

(50) |

so . This can be determined by the virial theorem since it is proportional to the kinetic energy of the SW, hence . Therefore the amplitude of the soliton satisfies

(51) |

Upon taking the mean and variance of the amplitude from Eq. (51), we recover the solutions derived using continuum approximation.

Acknowledgments NU acknowledges financial support from FOM.

## References

- (1)