Nonlinear waves in a model for silicate layers

CHAOS 28, 083119 (2018) Nonlinear waves in a model for silicate layers Juan F. R. Archilla,1,a)Yaroslav Zolotaryuk,2Yuriy A. Kosevich,3and Yusuke Doi4 1Grupo de Física No Lineal, Universidad de Sevilla, ETSI Informática, Avda Reina Mercedes s/n, 41012 Sevilla, Spain 2Bogolyubov Institute for Theoretical Physics, National Academy of Sciences of Ukraine, vul. Metrologichna 14-B, 03680 Kiev, Ukraine 3Semenov Institute of Chemical Physics, Russian Academy of Sciences, Kosygin street 4, 119991 Moscow, Russia and Plekhanov Russian University of Economics, 36 Stremyanny per., 117997 Moscow, Russia. 4Department of Adaptive Machine Systems, Graduate School of Engineering, Osaka University, 2-1 Yamadaoka, Suita, Osaka 565-0871, Japan (Received 21 March 2018; accepted 8 August 2018; published online 24 August 2018) Some layered silicates are composed of positive ions, surrounded by layers of ions with opposite sign. Mica muscovite is a particularly interesting material, because there exist fossil and experimental evidence for nonlinear excitations transporting localized energy and charge along the cation rows within the potassium layers. This evidence suggests that there are different kinds of excitations with different energies and properties. Some of the authors proposed recently a one-dimensional model based on physical principles and the silicate structure. The main characteristic of the model is that it has a hard substrate potential and two different repulsion terms, between ions and nuclei. In a previous work with this model, it was found the propagation of crowdions, i.e., lattice kinks in a lattice with substrate potential that transport mass and charge. They have a single specific velocity and energy coherent with the experimental data. In the present work, we perform a much more thorough search for nonlinear excitations in the same model using the pseudospectral method to obtain exact nanopteron solutions, which are single kinks with tails, crowdions, and bi-crowdions. We analyze their velocities, energies, and stability or instability and the possible reasons for the latter. We relate the different excitations with their possible origin from recoils from different beta decays and with the fossil tracks. We explore the consequences of some variation of the physical parameters because their values are not perfectly known. Through a different method, we also have found stationary and moving breathers, that is, localized nonlinear excitations with an internal vibration. Moving breathers have small amplitude and energy, which is also coherent with the fossil evidence. Published by AIP Publishing. https://doi.org/10.1063/1.5030376 Some fossil dark tracks in a layered silicate, mica muscovite, have been identified with localized excitations that travel through the rows of atoms within layers. Also, experiments show that these excitations are able to expel atoms from the surface of the mineral and transport charge. The different tracks that can be observed suggest that different kinds of excitations exist. In this paper, we use a model where the potentials have been obtained from the physical properties and the structure of the mineral. We make a systematic search for the possible nonlinear excitations and find different ones. Some imply the movement of charge and mass, and others have smaller energy and are only vibrational. Some are localized, but others are coupled to low amplitude plane waves. We explore the energies, velocities, and stability of these solutions and explain how they can be produced by the recoil after different radioactive decays of the isotope 40Kandrelate them with the different tracks. I. INTRODUCTION Muscovite is a layered silicate that has been from long ago of interest for different reasons. It was the second a)[email protected] material1for which nuclear tracks from fission fragments were observed, and the first material for which fossil nuclear tracks were discovered.2Also, fossil tracks from different elementary swift particles as positrons and muons could be observed.3These tracks were formed by the precipitation of excess impurity of Fe growing by accretion into the mineral magnetite and can be viewed without instruments. Many other tracks were in the closest-packed directions of the K+plane and were attributed to some kind of quasi one-dimensional lattice excitations, whose exact nature is not well known. See Ref. 4for a recent review. The existence of these excitations was also demonstrated experimentally by bombarding with αparticles a muscovite sample and detecting the ejection of an atom at the direction of the lattice closest-packed lines.5Recently, it was deduced that the lattice-excitations called quodons carry electric charge,4could be of many different types with different charge states,6and demonstrated experimentally.7Trying to understand the basic properties of the mineral that led to the existence of nonlinear excitations, a one-dimensional model was introduced with potentials deduced from basic principles and physical properties that led to the numerical finding of the existence of a moving kink or crowdion, that is, a moving interstitial, with a single energy which was the appropriate for the experimental and fossil observations.8,9 1054-1500/2018/28(8)/083119/15/$30.00 28, 083119-1 Published by AIP Publishing. 083119-2 Archilla et al .Chaos 28, 083119 (2018) In this article, we use powerful numerical techniques to perform a more systematic search of nonlinear waves in the same model. We have found other solitary waves as kinks, double kinks, and nanopterons with different energies and velocities. However, among them, only the faster entities as one single kink, two double kinks, and fast nanopterons are stable. Their approximate velocity is about 5 times the sound velocity in the system. We have also found stationary breathers that cannot be moved, and also small amplitude moving breathers that propagate smoothly through the lattice. The physical consequences of these findings are analyzed. The generic solutions found are nanopterons, i.e., nonlocal solutions that can be described as the coupling of a kink and a plane wave. Nanopterons can be considered as an approximation to a physical solution, but can also be of interest by themselves.10–14 Nonlocal breather solutions are also a subject of current research.15–17 It should be stressed that finding localized non-radiating travelling-wave solutions in lattices with an on-site potential was and remains an important problem in the physics of nonlinear waves. Discreteness of the media seriously obstructs existence of exact moving topological solitons.18 However, it has been shown in numerous cases both analytically19–22 and numerically23–26 that these solutions exist, also in twodimensional lattices,27–29 despite the obvious resonance with the plane waves of the system. Therefore, they fall into the family of the embedded solitons that are literally embedded in the linear spectrum of the system as defined in Ref. 30. Contrary to the conventional solitons that usually exist as a one-parametric family of solutions, where the parameter may be their velocity, frequency, or some other characteristic, the embedded solitons form a discrete set of solutions. In this paper, we are going to demonstrate the existence of nonradiating moving topological solitons in a model for mica muscovite with realistic potentials. The outline of the article is as follows: we describe the mineral and previous results in Sec. II and the details of the model in Sec. III. Section IV describes the peculiarity of the phonon spectrum and the phonon different velocities. Section Vexplains the pseudo-spectral method for calculating travelling solitary waves with precision. In Sec. VI, we describe the characteristics of the nanopterons and kinks obtained, and their velocities and energies. Their stability is analyzed in Sec. VII. In Sec. VIII, we explore the consequences of changing the strength of the interaction or the on-site potential. In Sec. IX, the properties of some double kinks or 2-crowdions are presented. Breathers are studied in Sec. Xand the article finishes with the conclusions in Sec. XI. II. MUSCOVITE AND NONLINEAR EXCITATIONS In this section, we describe the silicate structure of muscovite, the physics of the model, previous results, and some considerations about the possible existence of different nonlinear excitations in the material and the model. A. Silicate structure Most silicates are formed by the basic building block of a tetrahedron formed by the cation silicon surrounded by oxygen anions at the vertices. Often those oxygen anions are shared with other tetrahedra forming sheets of several layers, sometimes tetrahedra sheets, joined with other ones or with octahedra sheets, are sandwiched between tetrahedra sheets bringing about thicker layers. These sheets can be joined together by van der Waals forces. Often, substitutions of silicon by other cations like aluminium with different charges are common and, therefore, the sheets have a net electrical charge. Ions of different signs are, therefore, brought in to neutralize the charge and the forces between sheets become stronger and predominantly electrostatic. The result is a layer of ions that are coupled to the silicate sheets but less stronger than the silicate structure of the sheets. The inter-layer ions are kept in place at specific places where the attraction of the silicate layer is stronger. If we consider the ion lattice, it can be viewed as subjected to an on-site potential where the ions interact between them through electrical repulsion. It is a natural and important example where on-site potentials become essential in determining the properties of a lattice. The layered silicate mica muscovite is of particular interest in this article. It consists of a silicate layer with the structure of tetrahedraoctahedra-tetrahedra. Substitution of Si+4by Al+3results in a net negative electrical charge, which is compensated by the introduction of a cation layer made out of K+.Often,there are also different substitutional impurities which have different properties; however, we will refer in this article to the ideal structure. B. A realistic model for muscovite and the existence of acrowdion A model with realistic potentials was constructed starting from the simplest model with only electrostatic repulsion between K+ions31 and later augmented with a Ziegler- Biersack-Littmark (ZBL) short-range repulsive potential between K+ions,32 which was important for describing the interaction at very short distances bringing about reasonable minimal distances for energetic excitations.8,9In the latter references, it was also introduced an on-site potential calculated using empirical potentials33 and the interaction with the two closest two-dimensional layers of ions above and below the K+layers with some adjustment to take into account further layers. The resulting periodic potential, which we will use in this article, could be expressed as a cosine Fourier series with five terms and was able to reproduce approximately the absorption infrared frequency 110cm−1obtained experimentally34 and the 20 eV potential well obtained with molecular dynamics.35 This model led to the numerical finding of a crowdion, that is, a kink in a lattice with a periodic on-site potential, where a particle is moved to another site pushing the particle that sits there and becoming a replacement cascade, which can also be described as a moving interstitial.8,9That crowdion or lattice kink has several important characteristics from the mathematical and physical point of view: (a) it has a unique supersonic velocity and energy and travels long distances; (b) it survives to perturbations of the media as nonlinear waves and thermalization;9(c) its energy of about 26 eV is in the right order of magnitude for mica muscovite crystal, smaller than most of the recoils after βdecay of 083119-3 Archilla et al .Chaos 28, 083119 (2018) 40K,36,37 which has a maximum recoil of 42 eV for β−,9and larger than the typical ejection energies of 4–8eV of atoms in silicates,38 making a good candidate for the experimental results.5The property that supersonic crowdions have a discrete set of velocities is well known.39,40 The key point for the existence of supersonic crowdions is the presence of the on-site potential, which can be provided by an external system or simply by the rest of the crystal, if we consider, for example, a row of nearest-neighbour atoms as a subsystem subjected to the lattice potential. In a layered silicate crystal, the on-site potential is provided by the silicate sheets above and below the layer of potassium atoms. The high bond stiffness of the silicate layers make them dynamically decoupled from the K+layer. C. Other nonlinear localized waves in muscovite In this paper, we expand the search to other nonlinear localized waves in the same model for mica muscovite. The reason for that is partly theoretical, as a peculiar characteristic of the on-site potential is that it is hard for relatively small amplitudes and later becomes soft. This hard-soft characteristic does not appear as a mathematical curiosity but as a consequence of constructing a potential based on physical laws of a particular system. The same can be said about the interaction potential between ions, which is obtained simply by electrostatic repulsion and short range universal ZBL repulsion, which describes the electrostatic repulsion between the nuclei screened by the electron shell. On the other hand, the dark tracks in muscovite that are attributed to nonlinear lattice excitations are not equal; for example, some are strong and continuous and others are faint and discontinuous and are scattered at 60◦from the strong ones.41 It can also be deduced from the analysis of the first collisions after 40K recoil.6These secondary tracks cannot have the same energy as the primary one, which means there may be different excitations at play. We describe three types of lattice excitations for this model, kinks (without tails), nanopterons, which are kinks with tails, and breathers, which are vibrational localized oscillations. The terminology is not uniform partly due to different characteristics of a excitation depending on the variable used. The word soliton is used also for kinks, because the bond compressions and the velocities are solitons, i.e., their amplitudes tend to zero at ±∞, and breathers can be seen as envelope solitons with an internal vibration. Kinks and crowdions are the same entity in a system with an external potential as the one studied here. A method for finding kinks was developed in Ref. 23, where the terminology bound solitons was used because the profile of the bond compressions looks like the union of several solitons. Nonlocal kinks or nanopterons that are composed of a localized kink core coupled with a plane wave, which are solutions of the system, are constructed. Among them we search for solutions that have zero amplitude for the plane wave, which are proper kinks. Nanopterons are not proper kinks because they have oscillating plane wave asymptotics. They are, however, long-lived and appear in many different physical systems10 and are, therefore, interesting by themselves. Stationary breathers are found by continuation from the continuous limit. They are site-centered ones that are very stable and have a wide range of frequencies. They are, however, not movable but they can have relevance in other phenomena as the experimentally observed fast reaction kinetics in muscovite.42–44 However, small amplitude moving breathers are found with simple initial conditions. D. Nonlinear waves in two-dimensional lattices There is also interest in the study of the possibility of excitation and propagation of solitons, discrete breathers, and kinks in two-dimensional lattices. Different interatomic potentials, including discrete sine-Gordon and Lennard- Jones27,28,45,46 and Morse29,47 potentials, were used in these studies which confirm the existence and long lifetime of the kink solutions in two-dimensional lattices, including a twodimensional model proposed for mica muscovite.27,28 In the present work, we consider a one-dimensional lattice model which allows us to perform more rigorous study of the existence and stability of nonlinear excitations in a system with realistic interatomic potentials. This also allows us to obtain the values of magnitudes as frequency, velocity, and energy in physical units. III. A MODEL WITH A HARD SUBSTRATE POTENTIAL AND TWO REPULSIVE INTERACTIONS As stated above, the model and parameters correspond to the system described in detail in Refs. 8and 9. We consider a one-dimensional system of particles corresponding to ions K+separated by a distance rwhich in equilibrium is the lattice constant a=5.19Å which is set to 1 in scaled form with distance unit uL=a. The interaction between particles is given in scaled form as U(r)=1 r+Be−r/ρ r,(1) where the first term corresponds to the Coulomb repulsion between ions and the second Yukawa-type potential is the ZBL potential which describes the Coulomb repulsion between nuclei screened by the electron shell.32 The actual ZBL potential has four Yukawa terms that are important for different energy ranges from eV to MeV, but for the range of energies below 200 KeV one term is enough. The scaled unit of energy is the Coulomb energy at the equilibrium distance uE=Kce2/a≃2.77eV, and the constant B=184.1 is of the order of Z2,beingZ=19 the atomic number of K. The screening range is ρ=0.0569, which corresponds to ≃0.3 Å in physical units. The on-site potential has the periodicity of the lattice and is given in scaled units by V(u)= 4  m=0 Vmcos(2πmu),(2) with {Vm}=γ[2.4473,−3.3490, 1.0997, −0.2302, 0.0321]. It has been obtained using empirical potentials33 for the interaction with the two-dimensional silicate layers above and below the K+layers.8The parameter γ=1 will be used later to check the effect of a change on the strength of the on-site potential. 083119-4 Archilla et al .Chaos 28, 083119 (2018) Each particle with index nis described by a coordinate xn, which in equilibrium has the value x0 n=n. We perform a change of variables to un=xn−nwhich represents the deviation from the equilibrium position. In the variables un the interaction potential is given by U(un+1−un)=1 1+un+1−un +B 1+un+1−un e−(un+1−un+1)/ρ (3) without change for V(un)due to its periodicity. The kinetic energy of each particle is given by Kn=1 2mK˙u2 n,(4) with mK=39.1amu being the mass of a potassium atom. It will be taken as the unit of mass uM=mKin scaled variables, in which the kinetic energy is simply 1 2˙un2. The scaled unit of time uτcan now be obtained from uE=uM(uL/uτ)2,being uτ=(mka3/Kce2)1/2≃0.2ps. Therefore, the total energy is given in scaled units by H= n En= n1 2˙u2 n+V(un)+U(un−un−1).(5) The equations of motion have the following form: ¨un=−V(un)+U(un+1−un)−U(un−un−1),(6) with n=1,2, ...,N. The periodic boundary conditions un+N=un+Qare applied with Qbeing the topological charge that will be explained below. IV. PHONONS AND THEIR PHASE AND GROUP VELOCITIES The system described by Eq. (6) can be linearized for small amplitudes of unand the properties of phonons can be obtained easily.9They are typical of systems with an on-site potential but some of their properties and the specific values of some magnitudes are worth being remarked. The phonon dispersion is given by ω(q)=ω2 0+4c2 ssin2q 21/2,(7) with ω2 0=V(0)=γ 4  n=1 (2πn)2Vn, where qis the phonon wavevector. The frequency ω0=4.48 corresponds to linear oscillations of the isolated particle and cs=√2 is the sound velocity for the system without substrate, i.e., it is both the group and phase velocity for q→0, and it is also the maximum phonon velocity in the system without substrate. The wavevectors qcorrespond to plane waves cos[qn −ω(q)t]andsin[qn −ω(q)t] and due to the periodic boundary conditions they become discrete variables qn=2πm/N, for any integer m. As the wavenumbers qand q+2πnare equivalent, it is usually reduced to the first Brillouin zone q∈[−π,π] or the positive part of it for symmetry. However, the wavevectors in the extended zone are important in this system. The phonon spectrum is bounded from below by ω0and, therefore, the phase velocity Vph =ω(q)/q tends to ∞when q→0. This is an important property for Klein-Gordon system (systems with substrate potential) very different from FPU systems (without substrate). Note also that Vph is different for the wavevector qand q+2πin spite of representing indistinguishable phonons. This brings about that the phase velocities although bounded from below in the first Brillouin zone by Vph(π) ≃1.68, they are not in the extended zone. They have, however, a representative in the first Brillouin zone as can be seen in Fig. 1(top). However, the group velocity Vg=∂ω ∂q=c2 ssin(q) ω=c2 ssin(q) [ω2 0+4c2 ssin2(q/2)]1/2(8) FIG. 1. Top: Phonon properties: dispersion relation ω(q)(upper oscillating curve), phase velocity Vph (monotonic function), and group velocity Vg (lower oscillating curve). It can be seen that for any value of the velocity of a solitary wave there is a plane wave with the same phase velocity for which the wavenumber qand the frequency can be immediately deduced. The wavenumber qcan be considered in the extended zone but it also has a representative in the first Brillouin Zone. The group velocity Vgis zero at q=0andq=π, having the maximum value Vg≃0.4091 at q=0.4733π. Bottom: Plane wave associated with the sliding kink with velocity s=0.825 found in Sec. VI where it can be seen the plane wave with phase velocity V=sand q=1.77π, in the extended zone (smaller wavelength), and its representative in the 1st Brillouin zone, with q=q−2π=−0.23π(larger wavelength) and velocity V=−6.37. The minimum velocity in the first Brillouin zone is 1.68. 083119-5 Archilla et al .Chaos 28, 083119 (2018) is zero at q=0andq=±πand has a maximum value which is easy to calculate Vg,max ≃0.4091 for q=0.4733π. The phonon dispersion law and the phase and group velocities can be seen in Fig. 1. Therefore, for any velocity of the kinks there is always a plane wave with the same velocity, which can be coupled with it to form a nanopteron, in some cases the amplitude of the plane wave becomes zero, which corresponds truly to kinks. This is a different situation with FPU systems where the phonon velocities are bounded from above and supersonic kinks have velocities above the phonon ones. However, plane waves cannot transport energy as they are extended in the whole space, for that wave groups are necessary. In systems with substrate, “supersonic” means above the group velocities. In this work, we will find that all kinks and nanopterons are supersonic. V. PSEUDO-SPECTRAL METHOD FOR CALCULATING SOLITARY WAVES We look for the travelling-wave solutions of Eq. (6) with velocity sgiven by un(t)=u(n−st)≡u(z).Forthe unknown function u(z), the set of ordinary ODEs transform to the following differential-delay equation with the advanced and retarded terms: s2u(z)=−V[u(z)]+U[u(z+1)−u(z)] −U[u(z)−u(z−1)]. (9) Specific boundary conditions for this model will be discussed later in this section. Unfortunately, there is no clear analytical technique for finding solutions of such equations, so that we have to use numerical methods. A very accurate and simple method for solving these equations and finding soliton solutions has been developed in Refs. 48–50. We look for the solution u(z)in the form u(z)≃u(0)(z)+ k  m=1 cmsin2πmz L, (10) where the second term is a Fourier series of plane waves with wavenumbers qm=2πm/Land represents any periodic function with period Lin zwith coefficients cmthat have to be determined. The period Lis chosen very large so that the influence of boundaries would be negligible. The integer part of L will effectively play the role of the chain length in the original system Eq. (6).Thefirsttermu(0)(z)is an initial approximation which is used in order to enhance the initial guess and/or to enforce the appropriate boundary conditions. For this function, we specify the boundary conditions as follows: u(−∞)=a1,u(∞)=a2,a1,a2∈Z.(11) The topological charge is defined as Q=a2−a1,where Q=−1 is called an antikink and corresponds to a shift in position of a lattice site after the kink has passed, i.e., to a moving interstitial, while Q=1 corresponds to the opposite phenomena, i.e., a moving vacancy. For simplicity, we will use the term kink for both entities clarifying the topological charge when necessary. In order to extend the method to multi-soliton solutions24 the initial approximation u(0)(z)for a general P-soliton solution is chosen in the form u(0)(z)=a1+a2 2+a2−a1 2P l tanhμz+lδ 2, (12) where l=0,±2, ...,±(P−1)if Pis odd and l=±1, ±3,...,±(P−1)if Pis even. Here μand δare trial parameters which will be determined below for each iteration step by using various minimization techniques. The parameter μ describes the inverse width of one kink while the parameter δ can be considered as the distance between the adjacent kinks. We substitute the ansatz Eq. (10) into the original differential-delay-advance Eq. (9) and divide the interval [−L/2, L/2] into kcollocation points [due to Eq. (10)]. The solutions we look for are antisymmetric, so we can restrict ourselves to the interval [0, L/2]). As a result, we obtain a system of knonlinear algebraic equations on unknown coefficients cmwhich can be solved by various modifications of the Newton method. After finding the coefficients cm,we reconstruct the shape of the function u(z). VI. OBTENTION OF TRAVELLING NANOPTERONS AND KINKS WITH Q =−1 In this section, we apply the pseudospectral method defined above to the differential-delay equation Eq. (9) and discuss the main properties of its solutions. We investigate the simplest solution, namely, the compression kink (strictly antikink) with Q=−1. To monitor the amplitude of the oscillating tails of the nanopteron,10 we calculate the tail amplitude defined by A=max z∈[L/2−ξ,L/2] |u(z)−a2|, (13) which measures the amplitude of the tail oscillations in the right end interval ξ. The length of this interval was taken to be ξ=L/8. For non-oscillating (monotonic) solutions the function Eq. (13) should be very close to zero and in the limit L→∞this function tends to zero. The tails of the kink decay quite slowly after the slope and in order to be sure that u(z)is horizontal enough we choose Lto be approximately 10u(0)|a2−a1|. We find numerically that for the non-oscillating kink solution A(s)will be almost undistinguishable from the machine zero. Minimizing it for each value of system parameters, we find an appropriate value of the velocity sand the corresponding monotonic profile while the tail amplitude function Eq. (13) reaches a minimum which is almost equal to zero. A. Description of solutions: profiles, amplitudes, velocities, and energies The typical solution of Eq. (9) appears to be a non-local solution, a bound state of a kink with the plane wave that moves with the same velocity s, also known as a nanopteron.10 For the fixed value of Lwe trace the nanopteron solutions and monitor the amplitude in their tails A(s)defined in Eq. (13). This dependence A(s), which can be seen as a equation of state, is shown in Fig. 2and appears to be a set of pairs of branches that come from large values of Aand that join 083119-6 Archilla et al .Chaos 28, 083119 (2018) together at some local minima. It appears that by gradually increasing seach branch can be continued to an arbitrarily large value of A. The branches amplitude-velocity in each pair as seen in Fig. 2can be described as an upper one and a lower one. Each of these branches is a single-valued function and these branches join each other at certain bifurcation points where dA/ds=∞. Alternatively, we can consider the dependence given in Fig. 2as a function s=s(A). In that case, the pairs of these branches can be classified as two singlevalued functions of Athat join together at the bifurcation point ds/dA=∞. In this case, we should call them as faster and slower branches since for the fixed amplitude their velocities differ. Because the kink velocity is the variable parameter in our problem, we choose the former classification. Figure 2shows that for a fixed value of the kink velocity sthere exists a countable set of solutions as Ais increased. While sdecreases, the branches are packed more densely and each individual branch A(s)becomes more and more steep. If the velocity sis fixed, then two different solutions that belong to one pair of curves differ by the amplitude and the wavelength of the oscillating tail. Depending on the pair of branches, they may also differ by the behaviour near the solutions center. The typical nanopteron profiles are given in Fig. 3, and it can be observed that if a solution lies on the upper curve of the pair, we may have u(0)<0, otherwise u(0)>0, as shown in panel (b). However, for larger velocities the behaviour in the core remains very similar with both profiles having u(0)>0. The wavelength of the nanopteron tail approximately equals to 2π/q∗, where the wavenumber q∗is the root of the equation sq∗=ω(q∗)and ω(q)is the dispersion law given in Eq. (7). For the two nanopterons with s=3.9 given in Fig. 3, the respective wavelength equals 4.47 for the one marked by the thin line and 4.70 for the one marked by the thick line. The solution of the equation sq∗=ω(q∗)yields 2π/q∗=5.15. Close to the sliding velocity the correspondence is much better. Indeed, for the nanopteron from the lowest branch of Fig. 2at s=2.75, we obtain the wavelength 3.45 and the root of sq∗=ω(q∗)gives the same value up to the third digit. We should keep in mind that the asymptotics of the nanopteron tails are not exactly phonons, but rather nonlinear periodic waves. Certainly, for smaller amplitudes they become closer to phonons. The A(s)behaviour for this lattice is different from the analogous figure for the discrete Klein-Gordon lattice with harmonic interactions,51 where the dependence A(s)is generally single valued. It appears that within the range of velocities 0 <s <4, the dependence A(s)hits the A=0 axis at several points, namely, s1=2.73025, s2=1.51338, s3=0.82501, s4=0.67820, s5=0.42719, .... As far as velocities s>4 are concerned, we did not detect any zeroes of the A(s) dependence up to s=12. We did not continue into smaller velocity values because of the bad convergence of the numerical scheme, which is very likely due to entering into the phonon group velocities zone. Examples of non-oscillating solitons for two values of the sliding velocity are given in Fig. 4. Both solutions can FIG. 2. Left: Oscillation amplitude Aof the nonlocal solution of the nanopteron with topological charge Q=−1 (compression) as a function or its velocity s and for L=40. Note that there are some branches that touch zero, but not others. Also, for some branches there are two solutions that correspond to the same branch. Right: Zoom of two branches that touch zero giving rise to a sliding solution or lattice kink. 083119-7 Archilla et al .Chaos 28, 083119 (2018) FIG. 3. Profiles of the non-local moving solitons (nanopterons) that belong to the same pair of branches in Fig. 2(thick line: upper branch, thin line: lower branch) and have the same velocity s=3.9 (a) and s=1.5 (b). The inset shows the details in their core. be treated as a bound state of two acoustic solitons, and especially this can be seen well from the u(z)profiles. The difference between these solutions is the behaviour in the centre (z=0). As sdecreases, the separation between the global minima of u(z)becomes more significant. This remains to be true for the rest of the non-oscillating solutions, those with s<1. If the length Lis changed, the dependence A(s)also changes, but insignificantly. The sliding velocities [those at which A(s)=0] do not change. This is logical, because when A(s)=0 the boundary conditions become irrelevant. It has not been possible to find other non-oscillating bound states in spite of trying different δin the initial ansatz Eq. (12), with one, three, or other number of kinks in one bound state. FIG. 4. Profiles of the sliding kink solutions with Q=−1 (compression). [u(z)thick line, compression −u(z)thin line] for s=1.51338 (a) and s=2.73025 (b). The length of the chain is L=40. Note that the number of activated bonds, i.e., the length of the kink is approximately 4 at (a) and 3 at (b), the two largest velocities. B. Energies We set the zero of the energy at the equilibrium distances, that is, the energy is given by E=H−U0,whereU0is the sum of the potential interactions at the equilibrium distances, U0=[nU(un−un−1)]u=0. Note that the on-site potential V(un)is already zero for un=0. The values of scaled and physical velocities are given in Table I. Note that in spite of the different velocities of the sliding solitons, their energies appear to vary by less than 1%. This is logical because the ions need to have the same energy to overcome the same potential barrier and they need the same energy increase to be two of them fairly close in the same potential well. VII. STABILITY OF KINKS AND NANOPTERONS WITH Q =−1 Stability of the moving sliding kinks with Q=−1has been checked by integrating numerically the dynamical equations of motion where the exact or perturbed solution of Eq. (9) was used as an initial condition. The perturbation was introduced by distorting the velocity field ˙un(0)in several ways, for example, ˙un(0)→(1+)˙un(0), (14) where is the perturbation parameter. The displacement field un(0)was kept unchanged. Alternatively, the displacement field uncan be excited in the similar manner while the velocity field was undistorted. Although the two methods seem quite different, they lead to similar results for the stability as can TABLE I. Scaled and physical velocities and energies of the localized kinks with Q=−1. Note that in spite of the different velocities, their energies are very similar. Notation s5s4s3s2s1 V(scaled) 0.4271903 0.678198 0.825005 1.513382 2.73025 E(scaled) 9.177956 9.230015 9.241154 9.10141 9.427962 V(km/s) 1.1171 1.7734 2.1573 3.9574 7.1394 E(eV) 25.4648 25.6093 25.6402 25.2525 26.1585 083119-8 Archilla et al .Chaos 28, 083119 (2018) FIG. 5. Center of mass (17) evolution for: (a) the kink with the sliding velocity s=2.73025 with the randomly perturbed initial velocities according to Eq. (15) and ≃0.01 of the order of thermal perturbation at room temperature; (b) nanopterons with =0.001 in Eq. (15),s=2.8 and A=0.0079; (c) nanopteron with s=3andA=0.064. Panels (d)–(f) correspond to the same as (a)–(c) but with the initially perturbed displacement field as in Eq. (16). be seen in Fig. 6. The perturbation parameter has been varied from −0.1 to 0.1. The boundary conditions were periodic. Only the kink with fastest velocity s1=2.73025 is stable for positive values of =0.1 but gets pinned for negative =−0.01. This is due to the fact that being the perturbation proportional to the particle velocities in the kink core, which are large, the diminution of the velocities takes too much energy away from the kink. A different perturbation can be done with a random variation of the initial displacements or velocities: un(0)→(ξn−1/2)+un(0),n=1,2,...,N, (15) ˙un(0)→(ξn−1/2)+˙un(0),n=1,2, ...,N, (16) where {ξn}Nl n=1∈[0,1] are random numbers. This perturbation is of the order of the thermal velocities at room temperature for ≃0.1. The kink also shows significant robustness. The center of mass Xc(t)=N n=1nEn N n=1En , (17) dynamics is given in Fig. 5(a). Here, the energy density Enis defined in Eq. (5). The stability of nanopterons depends on many details such as speed, amplitude, and also the lattice size. Generically speaking, nanopterons with fast velocities s2.5 −3 can travel 10–25 times the whole lattice before getting pinned, while slower nanopterons s2 travel much shorter distances or get pinned in the lattice from the very beginning. If the nanopteron speed is close to the sliding velocity, it can be significantly robust, travelling successfully around the lattice with N=40 more than 60 times [see Fig. 5(b)]. Another nanopteron solution with the velocity s=3 and same perturbed initial conditions has managed to travel much shorter distance as shown in Fig. 5(c). In all three cases of Fig. 5, the periodic boundary conditions un+N=un+Qwere used. From the physical point of view, we can conclude that nanopterons can play a role as a transient state before the sliding kink velocity is established. It is, of course, impossible, to produce an exact nanopteron in a physical process as it should include a very large part of the lattice, but approximations to them might happen. To get better understanding of the nanopteron and kink stability, we have performed the global stability studies. For that, we need to follow each of the branches that appear in Fig. 2. We have used the homogeneous velocity perturbation, described in Eq. (14). In this way, we can exclude the dependence on the particular realization of disorder. We also have used the random perturbation scheme (16) for one fixed disorder realization. To characterize the robustness of the solitary wave, we introduce the survival time tsfor which the progressive propagation with the constant velocity is sustained: X(t)=X0+st,t<ts, pinned, t>ts . (18) For example, in the case shown by Fig. 5(a) ts∼800. The logarithmic dependence of this critical time is given in Fig. 6. We numerate the pairs of curves in Fig. 2from right to the left. Thus, the first pair is positioned in the top right corner of that figure. The pair number 4 contains the non-oscillating kink with s1=2.73025. The second non-oscillating kink belongs to the pair 15. Then we take each branch (upper and lower) from each pair, and compute the time tswhile the velocity is varied. For the first pair we have used three values of the perturbation =0.1, =0.01 and =0.001. This case is depicted in Fig. 6(a). As one can see, for the value =0.1 the original solution gets destroyed too quickly and not much is seen on the logarithmic scale. Therefore, we drop it for the next cases. The main results of the global stability studies are the following: • Only the nanopterons from the lower branches can be truly stable, and the survival time for the upper branch solutions never exceeds 102. That means it can hardly travel more 083119-9 Archilla et al .Chaos 28, 083119 (2018) FIG. 6. Survival time ts(see text) as a function of the velocity for different branches of the A(s)dependence represented in Fig. 2. The panels (a)–(f) correspond to the pairs of curves 1–6 starting from the right, i.e., from higher velocities. The lower curve is generally more stable than the upper one for the same branch. Lower branches: thick continuous line for the perturbation parameter =0.001, dashed for =0.01; upper branches: dash-dotted line for =0.001, dotted line for =0.01. Also for panel (a) there are two thinner lines at the bottom for =0.1, continuous for the upper branch and dotted for the lower branch. Circles (◦) mark results for the randomly perturbed chain described in Eq. (16) with =0.01. Vertical thinner lines mark the position of the local minimum of the respective A(s)dependence. than 100 lattice sites. On the contrary, the lower branch solutions can be extremely robust with ts>2×105(the computation was stopped at those times). • The intervals of high nanopteron stability are positioned at the lower parts of the A(s)dependencies (for small amplitudes). • High robustness of the nanopterons is seen for the highspeed branches, especially for the branches 1, 2, and 4. The branch 4 contains a non-radiating kink; however, the most stable solutions are nanopterons with velocities s∈[2.805,2.807] that slightly exceed the non-radiating kink velocity s1. This can be understood easily as the faster nanopteron can slow down because of disorder keeping its stability. The intervals of the extremely high kink robustness are very close to the respective minima of the A(s)curve for the pairs 1 and 2. Thus, the most stable nanopterons have the smallest tail amplitude. However, as commented above, this is not true for the branches that