Propagation of transition front in bi-stable nondegenerate chains: model dependence and universality
1 Propagation of transition front in bi-stable nondegenerate chains: model dependence and universality
I.B. Shiroky and O.V. Gendelman* Faculty of Mechanical Engineering, Technion * - contacting author, [email protected]
We consider a propagation of transition fronts in one-dimensional chains with bi-stable nondegenerate on-site potential. If one adopts linear coupling in the chain and piecewise linear on-site force, then it is possible to develop well-known exact solutions for the front and accompanying oscillatory tail. We demonstrate that these solutions are essentially non-robust. Various approximations for the on-site potential with the same basic parameters (height and coordinate of the potential barrier, energy effect and distance between the potential wells) lead to substantially different front velocities. Besides, inclusion of even weak nearest neighbor nonlinearity drastically modifies the front structure and parameters. The energy concentration in the front zone leads to a dominance of the nonlinear term. It turns out that the dynamics can be efficiently studied in terms of an equivalent model with a single degree of freedom. This estimation leads to accurate prediction of the front velocity and parameters of the oscillatory tail. Moreover, it turns out that the solution is robust - exact shape of the on-site potential weakly effects the front parameters. This finding also conforms to the simplified model, since the latter invokes only the general shape characteristics of the on-site potential. Introduction
In condensed matter, many processes accompanied by energy release, occur through propagation of transition fronts in the bulk at a nanosecond scale [1, 2, 3, 4]. Estimation of the characteristics of such transitions, especially of the front velocity and structure of the oscillatory tail, are of considerable interest in theoretical and experimental studies. It is beneficial in many cases to invoke discrete models for description of these phenomena. In such models, the transition front (or other structural defect) is to overcome the potential barrier caused by the discrete lattice. For this sake, it requires either external forcing or energy release in the lattice, related to the front propagation. This micro- scale effect disappears in continuous models, where the solutions, which represent defects, can move freely without drag. The discrete models, on the other hand, in principle allow derivation of a relation between the external force or energy gain and the front velocity, commonly known as the “kinetic relation”. At the micro level, when no dissipation is assumed, the effect, which describes the removal of the released energy, is known in literature as “radiative damping” [5]. The lattice waves that accompany the transition front have to remove at least the major portion of the gained energy, while the rest of the energy may be spent for creating new surfaces as it happens in the case of crack initiation [6, 7]. Analysis of the discrete models with multi-equilibria potentials of interaction goes back to early work of Dehlinger [8] and dynamical model of Frenkel and Kontorova [9]. Various derivatives of these models appeared throughout the years to describe different phenomena. Classic examples of this sort are dislocations in metals [3, 4], lattice distortions around twin boundaries [10], domain walls in ferro-electrics [1], and crack propagation [11]. Truskinovsky and Vainchtein in [2] address martensitic phase transitions by presenting a discrete model with long- range interactions. The model allows the derivation of a macroscopic dissipation law specified as a relation between force and velocity. The dissipation is due to radiation of lattice waves that carry energy from the front. Dynamics of crowdions in anisotropic crystals was studied in [12, 13]. In [14] a damped and externally driven FK chain is studied and the threshold forcing amplitude is found. This model may be beneficial in describing a chain of ions trapped in a metallic surface with an external AC electric field as the drive, as well as in the study of Josephson junctions. Other applications can be found in the comprehensive review of Braun and Kivshar [4]. A common configuration that has been adopted in many works is that of a bi-stable on-site potential [15, 16, 5]. The most interesting case is the non-degenerate on-site potential with certain energetic difference between the minima; this difference dictates the direction of the front propagation. The linearly coupled chain with a bi-stable on-site potential was studied in previous works. The simplest model where both wells have same curvature was proposed in the seminal work of Atkinson and Cabrera [3] and further explored in the early works of Ishioka [17] and Celli [18] . The advantage of this potential is that an analytical solution can be obtained through direct Fourier transform. The technique was further implemented in [2, 5, 19]. The case with different well curvatures requires application of the Wiener-Hopf method. It was done for systems with asymmetric parabolic double well in [20, 21, 22]. Similar models with nonlinear coupling received less attention. Recent studies [23, 24] address the case of generic coupling with on-site dissipation and suggest a law that connects the transported energy with the velocity and dissipation ratio. The studies mentioned above widely explored the models with linear and piecewise linear interactions. The reason is obvious – a possibility of exact solutions and comprehensive exploration of the transition front and the oscillatory tail. However, even if one restricts himself by the one-dimensional setting and considers a model with minimal possible set of interactions (for instance, only the on-site and nearest – neighbor), still a robustness of the linear approximation is rather questionable. In this work, we are going to demonstrate that the linear approximation is indeed not robust in a sense that the addition of small nonlinearity leads to qualitative changes in the structure of the propagating front and the oscillatory tail. Besides, even if the realistic system is reduced to the one-dimensional model, the exact shapes of the interaction potentials are hardly known and one can operate only with some particular characteristics. Commonly, only such basic parameters as the energy barrier, energy effect of the transition, the distance between the minima and, sometimes, the frequencies of oscillations in the potential wells can be measured or simulated. Any specific set of these parameters can correspond to 3 uncountable infinity of possible approximations of the potential functions. We are going to explore how the specific choice of such approximation affects the transition front characteristics. It turns out that the effect is noticeable in the case of the linear nearest-neighbor interaction and almost absent if the nonlinearity is taken into account. Description of the model and the case of the linear nearest-neighbor interaction . The basic model is adopted from [15, 5] where the problem of a chain with the linear inter-particle interaction and an on-site bi-parabolic potential is addressed. Here we expand this model in two ways. First, we include a cubic nonlinear inter-particle coupling. Second, we don’t specify the exact shape of the on-site potential, but only fix its main characteristics. The bi-parabolic shape will be a particular case that we will refer to later. The Hamiltonian of the system is given by Eq. (1). n n n n n nn p cH U (1) n is the displacement of the th n particle from the initial equilibrium state (meta-stable). c and are respectively the stiffnesses of the linear and nonlinear nearest-neighbor springs, n n p , the mass of each particle is set to unity. The shape of the nondegenerate bi-stable on-site potential U is defined by the energetic effect Q , the height of the potential barrier B , coordinate of the barrier b and coordinate of the stable state * . Three examples that conform to these requirements are shown in Figure 1. Position of the meta-stable state is set to zero. In this setting, the energetic effect Q is the driving force which determines the favorable direction of the reaction. Figure 1 - On-site nondegenerate potential ( ) U . Tree possible approximations with the same general shape parameters are presented: solid-blue-bi-parabolic potential, line-dotted green – 4 th order polynomial, dashed red – 6 th order polynomial. Without loss of generality the linear coupling coefficient is set to unity; c . The special case of system with Hamiltonian (1) with linear inter-particle interaction ( ) and a symmetric bi-parabolic potential with equal curvatures is studied in [5]. Linear dispersion relation for this model is given by: k (2) Here is a frequency of particle oscillations in each of the two potential wells. The typical dynamic response in this case is shown in Figure 2 for certain time instance. The only nonzero initial condition is the velocity of particle - . From here on, this condition is denoted in figures as “impulse 10”. It is seen that about 5 particles are within the transition area simultaneously. Figure 2 – Particle displacements in the chain with linear coupling. at t . Detailed structure of the front zone is presented in the inset. Parameters: / 6, 0.5 Q B , initial conditions: impulse 10
For stationary propagation of the transition it is necessary that the front velocity will be equal to the phase velocity of the accompanying oscillatory tail, f ph
V V .It is shown in [5] that the phase velocity can be expressed analytically as a function of the driving force Q through the following relationship:
22 2 00 k M k Qk k B (3) Here, M are the real and complex roots of dispersion relation (2) of the chain that satisfy : 0, Im 0 : 0, Im 0, 0 k M k L k k L k kL The phase velocity is implicitly determined for any set of the system parameters, as Eq. (3) can be satisfied only by a unique set of roots M that are found from the dispersion relation for a single value of frequency * . This frequency corresponds to a single value of real wavevector * k k Then, one obtains: 5 * * / f ph V V k (4) To examine robustness of solution (3)-(4) to variations in the shape of the on-site potential, we numerically integrate the evolution equations obtained from Hamiltonian (1), with different U that have common general shape ( * , , , B Q b ), as explained in Figure 1. The front velocity as a function of * is presented in Figure 3. For three different potentials: bi-parabolic, 4 th order polynomial and 6 th order polynomial, we obtain different results for the front velocity. For the 6 th order polynomial, the front propagation is not observed in some range of parameters, for which it is observed for other approximations for the potentials. This finding itself indicates a significant non-robustness of the linear chain to variations in on-site potential shape. Figure 3 – Phase velocity as a function of * ; Solid blue – bi-parabolic potential, ‘o’ red – 4 th order polynomial potential, ‘x’ green – 6 th order polynomial potential; parameters: / 1 Q B The nonlinear nearest-neighbor coupling.
Let us now consider the chain with Hamiltonian (1) and nonlinear nearest-neighbor stiffness . A typical dynamical response is shown in the n plane for fixed time instance (Figure 4), and in the t plane for one of the chain particles (Figure 5). All parameters and initial conditions besides β are similar to those used for the simulation presented in Figure 2. Figure 4 - Dynamic response of the chain with cubic coupling; n at t . Parameters: / 6, 0.5, 0.2 Q B , initial conditions: impulse 10
Figure 5 - Dynamic response of the nonlinear chain; t ; Parameters: / 6, 0.5, 0.2 Q B , initial conditions: impulse 10
One immediately observes that the structure of the solution strongly differs from the case of linear interparticle interaction, even for a small nonlinearity. The transition front is considerably accelerated. The gradient at the transition area is extremely steep, with at most 3 particles in the front zone (see zoom in Figure 4). Besides, the oscillatory tail has very large wavelength. To gain further insight in the front structure, we numerically evaluate the average contributions of the quadratic and quartic terms of strain energy in the chain with the help of the following equations: t tn n n nl nlt t t t t te n dt e n dt (5) Here ph V is the characteristic time of transition, l e and nl e are the average distributions of linear and nonlinear portions of interaction energy respectively. The numeric results are presented in Figure 6. For both quadratic and quartic components of the coupling energy, the concentration is extremely high in the narrow transition area compared to the rest of the chain. In addition, one notes that the 7 quartic term is responsible for about 6 times more energy concentrated in the front than the quadratic term. The characteristic time that is needed for the solution to acquire its steady state front velocity is rather large compared to the case of the linear chain (Figure 7). This can be explained by high amount of energy concentrated in the transition zone; a lot of time is required to accumulate this amount of energy in the transition zone and to react the steady-state propagation. In conclusion, it is clear that the nonlinear interaction term qualitatively modifies the transition front, and the latter cannot be considered as perturbation of the exact solution available for the linearly coupled chain. Figure 6 – Nearest-neighbor interaction energy of the chain with cubic coupling term. a) linear coupling l l e e n b) nonlinear coupling nl nl e e n . Common parameters: / 6, 0.5, 0.2, 950 Q B t , initial conditions: impulse 10
Figure 7 - Transient time as a function of ; parameters: / 1, 0.5 Q B
Simplified model of the transition front.
The simplified model of the front zone for the case of nonlinear coupling is based on three observations. The first is the dominance of the nonlinear term in the transition area compared to contributions of the linear coupling and the on-site potential. Thus, in equations of motion for the particles belonging to the transition zone we can neglect in the first approximation all terms besides the nonlinear coupling. The on-site potential will affect only boundary conditions of the obtained solution. Thus, the following approximate equations of motion are considered for the particles in the transition area: i i i i i (6) Then, the transition front is extremely narrow, and only very few particles take part in the transition simultaneously (see Figure 4). Therefore, it is necessary to consider only very few particles in Eqs. (6) Visually, one encounters 2-3 particles inside the transition zone, but sufficient description may be obtained by considering the rapid jump of a single particle from the meta-stable state to the vicinity of the stable state. Moreover, the gradient in the transitional region is extremely steep when compared to adjacent layout within the two wells. Thus, one can adopt that the single particle inside the front is attached to a fixed particle that still has not left the metastable position . For the steady state, the front velocity must be equal to the phase velocity of the oscillatory tail. As it follows from dispersion relation (2), very large phase velocities / ph V k are possible only if the wavenumber k is small, i.e. close to the left bandgap (see Figure 8). The energy in the oscillatory tail is very small compared to those concentrated in the front; therefore, one can admit that the transiting particle is attached to an immobile point with coordinate , defined as the first maximum of the oscillatory tail behind the transition front (see Figure 9). The resulting approximate single-DOF equation for the particle inside the transition front is presented in (7). Figure 8- A typical dispersion relation plot
Δ 9
Figure 9- Definition of (7) By substituting z , integrating over the entire motion range z and by using t t dt , the following solution is obtained for the SDOF model: dzz z (8) Evaluation of the integral and substitution V further yields: V (9) Here K is a complete elliptic integral of the first kind. Expression (9) can be represented decimally as V . Of course, the numeric coefficient in (9) should not be taken too seriously because of numerous approximations. The qualitative prediction to check is that the front velocity is related to the system parameters in the following way: V (10) However, it was found from numerical simulations that in order to obtain very good quantitative coincidence, it is enough to multiply (9) by constant coefficient . Necessity of this correction stems from the simplification and assumptions taken; in the same time, with good accuracy, it is a constant, which does not depend on any of system parameters. With this correction, one can also obtain an approximate equation for the time history of the particle inside the propagating front z z t : z dzt zz z (11) Here the values of the actual time and can be obtained from: , t t z . Eq. (11) describes the movement of the considered single particle inside the front. The absence of system parameters in the rescaled representation implies that the basic structure of the kink does not depend on the system parameters and any specific solution is a rescaling of this basic shape. Comparison of this rescaled shape of the transition region to the one obtained numerically in the complete system is presented in Figure 10. The coincidence is reasonable, albeit not perfect. 0 Figure 10 – Comparison between analytic prediction and numeric result for the rescaled front shape.
Kinetic relation (10) is even more important. It implies that the following relations should be satisfied with sufficient accuracy:
2) ln 1 ln ln 2 K 41 2) ln ln ln 2 K2 4 a Vb V (12) Direct comparison between the numerical and the analytical results for the front velocity are presented in Figure 11 and Figure 12 as functions of and respectively. Figure 11 – Validation of approximation of the nonlinear chain ln ln ph V f ; Line-dotted blue – Numerical, Solid red – Analytical; parameters: / 1, 0.5 Q B
One of the main assumptions that led to derivation of the simplified model was the possibility to neglect the particularities of the on-site potential. The characteristics of the potential effected only the fixed value of the “amplitude of transition” . Each on-site potential yields different value of ; however, the scaling law (13) that relates between the front velocity and Δ is valid for broad range of 11 the on-site potentials with similar general characteristics ( * , , , B Q b ). To check this claim, we examine three different potentials: bi-parabolic, 4 th and 6 th order polynomials with similar characteristics (see Figure 1). The comparison between these simulations and the analytic model is also presented in Figure 12. It is seen that all potentials provide same ratio between and ph V . This supports the assumption that the phase velocity depends only on the coordinate , rather than on the exact potential shape. This weak dependence leads to an important conclusion, that the model (10) is robust for modified shapes of the on-site potential. Also, the approximate model provides good estimation of the numerical results with accuracy that asymptotically improves at higher velocities. Figure 12 – Front velocity as a function of potential parameter ; ‘o’ blue – parabolic potential, ‘*’ red – 4 th order polynomial, ‘x’ green – 6 th order polynomial; solid black – Analytical; Parameters: As it was demonstrated above, the front velocity linearly depends on . The latter parameter, however, is not known for a generic on-site potential, as it depends on the specific oscillation regime within the stable well. In the special case of bi-parabolic potential it is possible to evaluate this parameter explicitly. For this sake, we derive a nonlinear dispersion relation for the oscillatory tail. The piecewise parabolic potential has the following shape: n nn n n BbU b (13) Here, * 0 0
Q B B . Let us introduce complex variable: n n n i (14) The derivatives in terms of are written as follows: * * * , ,2 2 2 n n n nn n n n n n ii (15) 2 Substitution into the equations of motion for the stable branch and substitution of the modulated harmonic function: , i t i t i tn n n n e e i e yields:
2* * * * *0 1 1 1 13 3* * * *1 1 1 13 i t i t i t i t i t i t i tn n n n n n n n n n ni t i t i t i tn n n n n n n n ii i ie e e e e e ei e e e e (16) Division by i t e and subsequent averaging leads to the following slow-flow equations: n n n nn n n n n n n n n n i (17) By substituting iknn Ae one obtains the desired nonlinear dispersion relation:
22 2 2 2 20 2 k A kk (18) We see that the dispersion relation is indeed a small correction for the dispersion relation for linear chain (2). As it was mentioned above, the phase velocity of the oscillatory tail is large, consequently, the wavenumbers are near the left bandgap and the group velocity is small. Therefore, in the basic approximation the energy transport through the oscillatory tail can be neglected and the energy released due to the front propagation is almost not transferred towards or from the front. This situation strongly differs from one observed for the case of the linear coupling. [15]. The energy balance for arbitrary particle n in the oscillatory tail can be simply expressed as:
22 *2 0 nnn E (19) Due to the monochromatic structure of the tail, its form can be approximated as follows: * cos n A kn t (20) By substitution of (20) into (19) one obtains: sin2 n AE kn t (21) We plug in the nonlinear dispersion relation (18) into (21). In the limit k one obtains: QA (22) In Figure 13 the amplitude is shown as a function of nonlinear coupling parameter . The amplitude converges very quickly to the asymptotic value. At value of the convergence is within of the limit value. 13 Figure 13 - Amplitude of oscillation as a function of ; parameters:
1, 1, 0.5
Q B
Thus, for the case of bi-parabolic on-site potential one obtains the explicit expression: * 0 Q (23) Concluding remarks.
First, the results presented above demonstrate the lack of robustness of the models with linear coupling. The addition of even small nonlinearity to the nearest-neighbor interaction drastically modifies the dynamic response and causes extremely high velocity and strong localization of the transition front. This effect can be attributed to combination of focusing properties of the gradient nonlinearity considered here, and the non-degenerate on-site potential. The latter provides constant supply of energy, which counterweights effects of pinning and radiation. As a result, extremely narrow kinks can propagate through the chain with far supersonic velocities. Obviously, a continuous model of the system cannot explain such process. In the same time, namely the extreme localization and significant energy concentration in the front lead to considerable simplifications and reduction of the system to the effective single DOF model. Then, the dominance of the nonlinear term allows further simplifications and explicit prediction of the scaling relationships between the front velocity and other system properties: nonlinear coupling coefficient and maximal displacement of the particles inside the oscillatory front. These scaling relationships turn out to be universal for the on-site potentials with the same shape characteristics. The authors are very grateful to Israel Science Foundation (grant 838/13) for financial support.
References [1]
D. A. Bruce, "Scattering properties of ferroelectric domain walls,"
Journal of Physics C, vol. 14, pp. 5195-5214, 1981. [2] L. Truskinovsky and A. Vainchtein, "Kinetics of martensitic transitions: lattice model,"
SIAM Journal of Applied Mathematics, vol. 66, no. 2, pp. 533-553, 2005. [3] W. Atkinson and N. Cabrera, "Motion of a Frenkel–Kontorova dislocation in a one-dimensional crystal,"
Phys. Rev. A, vol. 138, no. 3, p. A763–A766., 1965. [4] O. M. Braun and Y. S. Kivshar, "Nonlinear dynamics of the Frenkel-Kontorova model,"
Physics Reports, vol. 306, pp. 1-108, 1998. [5] O. Kresse and L. Truskinovsky, "Mobility of lattice defects: discrete and continuum approaches,"
Journal of the Mechanics and Physics of Solids, vol. 51, pp. 1305-1332, 2003. [6] A. M. Balk, A. V. Cherkaev and L. I. Slepyan, "Dynamics of chains with non-monotone stress–strain relations. I. Model and numerical experiments,"
Journal of the Mechanics and Physics of Solids, vol. 49, no. 1, pp. 131-148, 2001. [7] A. M. Balk, A. V. Cherkaev and L. I. Slepyan, "Dynamics of chains with non-monotone stress–strain relations. II. Nonlinear waves and waves of phase transition,"
Journal of the Mechanics and Physics of Solids, vol. 49, no. 1, pp. 149-171, 2001. [8] U. Dehlinger, "Zur Theorie der Rekristallisation reiner Metalle,"
Annalen der Physik, vol. 394, no. 7, p. 749–793, 1929. [9] Y. I. Frenkel and T. Kontorova, "On the theory of plastic deformation and twinning,"
Phys. Z. Sowietunion, vol. 13, pp. 1-10, 1938. [10] M. Suezawa and K. Sumino, "Lattice distortion and the effective shear modulus along a coherent twin boundary,"
Physica Status Solidi (A) Applied Research, vol. 36, pp. 263-268, 1976. [11] M. Marder and S. Gross, "Origin of crack tip instabilities,"
Journal of the Mechanics and Physics of Solids, vol. 43, no. 1, pp. 1-48, 1995. [12] A. I. Landau, A. S. Kovalev and A. D. Kondratyuk, "Model of Interacting Atomic Chains and Its Application to the Description of the Crowdion in an Anisotropic Crystal," physica status solidi (b), vol. 179, no. 2, p. 373–381, 1993. [13] A. S. Kovalev, A. D. Kondratyuk, A. M. Kosevich and A. I. Landau, "Theoretical Description of the Crowdion in an Anisotropic Crystal Based on the Frenkel-Kontorova Model Including and Elastic Three-Dimensional Medium," physica status solidi (b), vol. 177, pp. 117-127, 1993. [14] L. L. Bonilla and B. Malomed, "Motion of kinks in the ac-driven damped Frenkel-Kontorova chain,"
Phys. Rev. B, vol. 43, no. 13, pp. 11539-11541, 1991. [15] V. V. Smirnov, O. V. Gendelman and L. I. Manevitch, "Front propagation in a bistable system: How the energy is released,"
Physical Review E, vol. 89, 2014. [16] L. I. Slepyan, "Dynamic factor in impact, phase transition and fracture,"
J. Mech. Phys. Solids, vol. 48, pp. 927-960, 2000. [17] S. Ishioka, "Steady motion of a dislocation in a lattice,"
J. Phys. Soc. Jpn., vol. 34, pp. 462-469, 1973. [18] V. Celli and N. Flytzanis, "Motion of a Screw Dislocation in a Crysta,"
J. Appl. Phys., vol. 41, no. 11, pp. 4443-4447, 1970. [19] A. Vainchtein, "The role of spinodal region in the kinetics of lattice phase transitions,"
Journal of the Mechanics and Physics of Solids, vol. 58, pp. 227-240, 2010. [20] L. I. Slepyan and L. V. Troyankina, "Fracture wave in a chain structure,"
J. Appl. Mech. Tech. Phys, vol. 25, no. 6, p. 921–927, 1984. [21] L. I. Slepyan, "Feeding and dissipative waves in fracture and phase transition: II. Phase-transition waves,"
Journal of the Mechanics and Physics of Solids, vol. 49, no. 3, p. 513–550, 2001. [22] O. Kresse and L. Truskinovsky, "Lattice friction for crystalline defects: from dislocations to 15 cracks,"
Journal of the Mechanics and Physics of Solids, vol. 52, p. 2521 – 2543, 2004. [23] N. Nadkarni, A. F. Arrieta, C. Chong, D. Kochmann and C. Daraio, "Unidirectional Transition Waves in Bistable Lattices,"
Physical Review Letters, vol. 118, p. 244501, 2016. [24] N. Nadkarni, C. Daraio, R. Abeyarante and M. Kochmann, "Universal energy traansport law for dissipative and diffusive phase transitions,"