Analytical solutions and collision stability analysis of solitary waves in a pre-compressed one-dimensional granular crystal

In this paper, the governing equation in a pre-compressed one-dimensional granular crystal, which was previously discussed by Nesterenko [J. Appl. Mech. Phys. 24, 733 (1983)], is solved analytically. Multiple solitary wave solutions are obtained by using the homogeneous balance principle and Hirota’s bilinear method. We analyze the difference between the original system and the KdV system and examine the collision of solitary waves in some special parameters. The dynamic behavior and stability of the double solitary waves are also studied. We find that the opposite collision between single solitary waves may be stable and thus generate a stable double solitary wave. It is concluded that the collision is a special stable double solitary wave solution. We further propose a possible way to determine the stability of multiple solitary waves qualitatively.


introduction
One distinguishing feature of wave propagation in a linear periodic structure (also termed phononic crystal) is the appearance of "stop bands" [1][2], which can be used to control the propagation of waves. These stop bands allow phononic crystals to serve as mechanical filters [3], waveguides [4][5], diodes [6][7][8][9] and resonators [10]. Although we have found a lot of interesting phenomena and extensive applications in linear phononic crystals, nonlinear phononic crystals have attracted much attention in recent years due to their tunable phononic characteristics [11]. Moreover, nonlinear phononic crystals have some phenomena that linear phononic crystals do not have, such as solitons [12][13][14], second harmonic [15], etc.
The simplest case of nonlinear phononic crystals is one-dimensional granular chain of mass [27]. The Fermi-Pasta-Ulam (FPU) problem was proposed from such nonlinear systems by Fermi et al. [28]. Later, the Korteweg-de Vries (KdV) equation [29] was derived from the FPU lattice. The soliton solution for the KdV equation explains the FPU problem under certain conditions. There are not only single solitons but also multiple solitons in KdV equation, which makes it very important in describing the dynamic behavior of weakly nonlinear systems. It should be noted that in the long wavelength approximation, the discrete system (e.g. the granular chains) can be described by KdV (or m-KdV) equation [27]. The FPU chain is a model for nearest-neighbour coupled with harmonic plus cubic (or quartic) nonlinear interactions. The long wavelength approximation gives rise to a Boussinesq equation or a nonlinear Schrodinger (NLS) equation which admit soliton solutions. Certain studies have also employed this approximation [30].
Granular crystals consisting of tightly packed aggregates of particles in a periodic arrangement have evinced special interest due to their tunable nonlinear [27]. It is shown that the interactions between spheres are highly nonlinear under weak precompression but are weakly nonlinear under higher precompression. Recently, these nonlinear systems have been studied extensively as they can be used as an ideal experimental platform for the study of interplay between discrete and continuum systems. There are a lot of mathematical methods for studying wave phenomena in nonlinear continuous systems, such as the perturbation techniques [31][32]. However, discrete nonlinear systems have received little attention. Due to the complexity of discrete systems, the research methods are mostly focused on numerical simulation. And the study of discrete systems based on analytical solutions is very few.
In this paper, the continuous equation corresponding to the discrete granular system, namely Boussinesq equation [30], is studied by solving the equation analytically. Then, the analytical solution is used as the initial signal to study the dynamic behavior and stability of the solution of the original discrete granular system. In Ref. [33], a continuous equation (Eq. (2.3) in Ref. [33], which belongs to Boussinesq equation [30]) in the original system without coordinate transformation was obtained in an initially compressed chain of granular spheres. In Ref. [34], we obtained the exact analytical single solitary wave solution from this continuous equation (Eq. (6) in Ref. [34]). However, only approximate analytical solutions of multiple-solitary waves were constructed without rigorous mathematical derivation. Moreover, there are some restrictions on the solutions, e.g. they do not hold when the wavenumbers of arbitrary two single solitary waves are opposite. In order to resolve the above problems, the bilinear method and homogeneous balance principle are suggested to solve the continuous equation directly. Based on the obtained analytical solutions, the dynamic behavior and stability of the double-solitary wave are studied in detail.

Continuous equation and analytical solutions
The governing equation for a monatomic granular chain of spherical particles under the pre-compression 0  % can be written as [33] 3/2 3/2 where n u % indicates the displacements of the n-th spheres; and the superscripted dot denotes the derivation w.r.t time t % ; and the Hertzian constant [33] with R % , m, E and m  being the radius, mass, Young's modulus and Poisson's ratio of the spheres, respectively. For convenience of discussion, we deal with the governing equation (1) in the dimensionless form. To this end, we set L to be a given characteristic length scale, and introduce the following dimensionless variables: and three dimensionless parameters Then Eq. (1) can be transformed into the following dimensionless equation in terms of the dimensionless variables and parameters: It should be indicated that the superscripted dot in the above equation represents the derivation w.r.t the dimensionless time t.

 
; and x is the dimensionless special coordinate.
Equation (5) belongs to Boussinesq equation [30] and can be solved by using the homogeneous balance principle [36][37][38][39]. To this end, we can assume that the derivation of the solution to Eq. (5) w,r,t x is 0 11 Substituting Eq. (6) into Eq. (5), we can obtain 22 1 , To reduce the calculation, we integrate Eq. (7) w.r.t x by setting the integral constant to be zero, and then have 22 Considering the following relations: we can rewrite Eq. (8) as   2  2  2  2  , , 0 , 0 , , , where D is the Hirota bilinear operator [35]. The above equation can be solved by using Hirota's bilinear method. Expand f into a power series of a small parameter ε:  (13) Setting the coefficient of the term for each order of ε in Eq. (14) to be zero, we obtain 2 2 2 4 00 of which the first four equations are Without of loss of generality, we assume 0 f to be a constant and Eq. (16) has an exponential eigen-solution: where k and  may be interpreted as the wavenumber angular frequency, respectively; and  is an arbitrary constant. Substitution of (19) into Eq. (16) yields from which we obtain the dispersion relation: The sum of the eigen-solution, Eq. (19), w.r.t all or partial possible values of j (j = 1L N) also satisfies Eq. (16). Therefore, we have the following eigen-solution: 1 1 , from which we obtain 20 1 Following the same process, we can obtain All the other terms, Here as examples, we present the detailed expressions of the specific solutions with ε=1 for N=1,2,3,4.

1) The bright single-solitary wave solution
When N = 1, we set 2 0 f  and 0 1 f  , and then have . Finally, we can obtain, from Eq.

2) The bright double-solitary wave solution
When N = 2, we set 3 0 f  and 0 1 f  , and then have  i  ); and and  A and 12 A is observed in Fig. 1(a) near the point 21 ( 1) kk becomes far away from this point, 12 A and 12 A are close to each orther. In Fig.  1 Although it seems that the corresponding dark multiple-solitary wave solutions can be obtained by integrating the bright ones (Eqs. (28), (30) and (31)) based on the relation / v u x    , it is not easy to obtain the closed-form expressions. However, as we mentioned above, the derivation of the multiple-solitary wave solutions in Ref. [34] (Eqs. (31)-(33)) w.r.t x is the same as the solutions, Eqs. (28), (30) and (31), in form. Therefore, we can directly write the dark multiplesolitary wave solutions corresponding to the bright ones (Eqs. (28), (30) and (31) where C is a constant; and 12 A is given by Eq. (29). The other dark multiplesolitary wave solutions can be obtained similarly.
In addition, the multiple-solitary wave solutions obtained here still have the dynamic characteristics of the solutions to the KdV equation described in Ref. [ (ii) 1   A bright single solitary wave evolving in the space-time domain is observed in Fig. 2(a). Its dispersion curve shown Fig. 2   In fact, Fig. 2 and Fig. 3 describe the same physical phenomenon, but in different ways. However, the dispersion curve in Fig. 3(b) is relatively simple. It seems clearer and more intuitive to describe the physical phenomena of solitary waves with the bright solitary waves. Thus, the following analysis will be mainly focused on the bright solitary waves.
For the complicated double-solitary wave, we consider the following two different cases: (I) 12 [41]. However, the Type-I double-solitary wave does not present a phase shift, but exhibits a sharp peak at the center of the collision between two single-solitary waves. Figure 4(b) illustrates the dispersion curves of the Type-I bright doublesolitary wave corresponding to Fig. 4(a). It is shown that the dispersion curves are concentrated in a cross-like region near two diagonal lines passing through the zero point. The energy transfer during the dynamic collision is shown Figs. 4(c) in the space-frequency domain and in Fig. 4(d) in the time-wavenumber domain. Before the collision, the energy flows from high frequencies/wavenumbers to low frequencies/wavenumbers; during the collision process, the energy concentrates to the low frequency/wavenumber region very quickly; and after the collision, the energy flows from low frequencies/wavenumbers to high frequencies/wavenumbers. For the elastic collision, the total energy is conserved. Figure 5 shows the evolution of the Type-I dark double-solitary wave [Eq. (35)] in the space-time domain corresponding to Fig. 4. It can be seen that the waveform is very similar to a group of steps with obvious changes in amplitude.  show that before the collision, the energy flows from the high frequencies/wavenumbers to low frequencies/wavenumbers. During the collision process, the energy is rapidly concentrated to the low frequency/wavenumber region. After the collision, the energy flows from the low frequencies/wavenumbers to high frequencies/wavenumbers. The total energy is conserved. The Type-II dark doublesolitary wave [Eq. (35)] corresponding to Fig. 6(a) is shown in Fig.7. Its waveform, although exhibits steps with changes of amplitude, is obviously different from that of the Type-I dark double-solitary wave shown in Fig. 5.
Through the above analysis, we can summarize the differences between the two types of double-solitary waves: Type-I can be seen as two single-solitary waves colliding in the opposite direction (cf. Figs. 4 and 5) without phase shift; and Type-II can be regarded as two single-solitary waves colliding in the same direction (cf. Figs. 6 and 7) with a phase shift after collision.
The same method can be used to analyze the triple-and quadruple-solitary wave solutions. The results are similar and will not be presented here.

III. Numerical analysis of stability of double-solitary waves
The stability of double-solitary wave will be studied numerically by following the idea proposed in Ref. [34]. Here we consider the collision of the following two dark single-solitary waves running in the opposite directions: which, according to the previous analysis, should generate the following Type-I dark double-solitary wave: x t x t To examine the evolution stability of the dark double-solitary wave, we use Runge-Kutta method to solve Eq. (4) numerically by assuming an initial excitation which is the displacement field from Eq. (38) at time t = 0 with a uniformly distributed random perturbation of amplitude 10 −4 [41][42]. The evolution of the dark double-solitary wave solution obtained from the Runge-Kutta method is demonstrated in Fig. 8(a) Fig. 8(b) with the same parameters as in Fig. 8(a). The similarity of Figs. 8(a) and 8(b) implies that the solitary wave in Fig. 8(a) is the Type-I doublesolitary wave. x t x t x Following the same process as above analysis, we solve Eq. (4) numerically using Runge-Kutta method with the initial excitation from the perturbed initial values of Eq. (40). The results are shown in Fig. 9(a). Before the collision, the waveforms of the two single-solitary waves are well preserved, and the perturbation is not amplified. At 700 t  , the two solitary waves collide, and the perturbation increases exponentially. After the collision, the waveform of the solitary wave becomes distorted obviously. Therefore, the collision is unstable, and the excited double-solitary wave are also unstable. The space-time evolution of Eq. (41) is shown in Fig. 9(b) with the same parameters as in Fig. 9(a). Comparing these two figures, we find that the solitary wave in Fig. 9(a) is the Type-I dark double-solitary wave, but it is not convergent and unstable. Further numerical tests show that the generated double-solitary wave can propagate stably for a longer time (i.e. the life-span becomes longer [43]) if the velocity difference of two single-solitary waves are smaller. This is an interesting phenomenon in granular phononic crystals. It seems that only the opposite collision between two single-solitary waves with the same velocities, as shown in Fig. 8 with Runge-Kutta method. The data from the Type-II double-solitary wave shown in Fig. 7 at a certain time is used as the input signal (see Fig. 10(a)). The perturbation magnitude of 10 -4 is added to the signal. Figure  10(b) shows that a double-solitary wave is excited in the original system before the time is less than 1000, and the waveform keeps good. But after t=1000, the perturbation is amplified exponentially, and the waveform loses the characteristics of the solitary wave. It is obvious that the Type-II double-solitary wave is unstable. Further numerical calculations support this conclusion although we cannot exhaust all cases. Collision of moving solitary waves is an interesting phenomenon which is closely related to the stability of the single solitary waves. However, our analytical solutions and corresponding numerical simulations show that the stability of collision is related to the stability of double-solitary waves. In other words, only when there is a stable double-solitary wave, is the collision of two single-solitary waves stable.
Next, we will present qualitative analysis about the stability of multiple-solitary waves from the point of view of geometry.

IV. Qualitative analysis of stability of multiplesolitary waves
According to the analysis in Section III, the multiple-solitary waves can be regarded as (or equivalently be excited by) the dynamic interaction (collision) of single-solitary waves. We first consider the case of a double-solitary wave. The two kinds of double-solitary waves (Type-I and Type-II) can be schematically illustrated in Fig. 11 where the strips represent the peak regions of the dark singlesolitary waves or the rapid-changing regions of the bright single-solitary waves. In the other regions, (I)-(IV), the amplitude tends to zero/constant for the bright/dark double-solitary wave. So, except in the strongly interacting regions marked by area S, the double-solitary wave can be regarded approximately as the superposition of two single-solitary waves (also refer to Ref. [34]). In the strongly interacting region, the two solitary waves interact sufficiently and generate intense energy exchange and flow. We find that bigger the area is, the longer the running time of the simulation program is. Thus, it is reasonable to infer that the stability of the doublesolitary wave is relevant to the area of the strongly interacting region, and the smaller the area is, the more stable the system is.
where  ( 0 90     ) is the angle between the two single-solitary waves; and d1 and d2 are the effective widths of the single-solitary waves. It is easy to see that S is a monotone decreasing function of  with the minimum value of min   . If the double-solitary wave is stable at an angle smaller than  , then it is also stable at the angle  . The bigger the angle  is, the smaller the area S is, and therefore the more stable the double-solitary wave is. 90   yields the minimum S. The Type-I double-solitary wave includes the case of 90   , which means it may be stable. It is obvious that the angle  of Type-II doublesolitary wave is smaller than that of Type-I and is generally much smaller than 90. Therefore, the Type-II double-solitary wave is generally unstable.
It is worth noting that the width of solitary waves is the diameter of five spheres [33]. Since the space width of solitary wave is also determined after the sphere chain is selected, we ignore the influence of the width of solitary wave on the stability when we consider its geometric structure. Next, we consider the case of a triple-solitary wave which can be treated as the dynamic interaction of three single-solitary waves. There are two situations: one single-solitary wave propagates against the other two as shown in Fig. 12(a), and all three single-solitary waves propagate in the same direction as shown in Fig.  12(b). In both situations, we will definitely have two single-solitary waves propagating in the same direction. And the above analysis demonstrates that the collision of these two single-solitary waves generates an unstable Type-II doublesolitary wave. Therefore, the triple-solitary wave is generally unstable. Indeed, we could not find a stable triple-solitary wave through many numerical tests.
Similarly, a multiple-solitary wave higher than the third order is generally unstable.

V. Conclusions
In this paper, the analytical solutions of multiple-solitary waves in a precompressed spherical chain are derived based the homogeneous balance principle and bilinear method. The dynamic behavior and stability of the double-solitary waves are studied in detail. The main results and conclusions may be summarized as follows: 1) As in Ref. [34], the present solutions also show that the multiple-solitary waves and be regarded as the dynamic interaction (collision) of multiple single solitary waves. The obtained solutions have the same form those in Ref. [34] which were constructed through the inspiration of the solutions to the KdV system. However, a great difference between the present solutions and those in Ref. [34] appears when the wavenumbers of arbitrary two single solitary waves are opposite (i.e. when ij kk  in Eq. (34)). This implies that the KdV system cannot exactly capture the nature of the original granular system.
2 and Type-II. Numerical simulation shows a stable Type-I double-solitary wave.
Other high-order multiple-solitary waves can be analyzed in the similar way.
3) The strongly nonlinear interaction between single solitary waves takes place in a small spatio-temporal region. We argue that the smaller this region is, the more likely the system is to be stable. This provides a possible way to determine the stability of multiple-solitary waves qualitatively.

VI. Acknowledgements
This piece of work is supported by the National Natural Science Foundation of China under grant numbers 11532001 and 11991031. The third author (Y. S. Wang) is also grateful for the support of National Natural Science Foundation of China under grant number 12021002.

VII. Declarations
We would like to declare on behalf of our co-authors that the work described in this paper is original, has not been published previously, and is not considered for publication elsewhere, in whole or in part. All the authors listed have approved the submission of the manuscript. No conflict of interest exits in the submission of the manuscript.