Wave propagation properties in oscillatory chains with cubic nonlinearities via nonlinear map approach

Free wave propagation properties in one-dimensional chains of nonlinear oscillators are investigated by means of nonlinear maps. In this realm, the governing diﬀerence equations are regarded as symplectic nonlinear transformations relating the amplitudes in adjacent chain sites ( n , n + 1) thereby considering a dynamical system where the location index n plays the role of the discrete time. Thus, wave propagation becomes synonymous of stability: ﬁnding regions of propagating wave solutions is equivalent to ﬁnding regions of linearly stable map solutions. Mechanical models of chains of linearly coupled nonlinear oscillators are investigated. Pass-and stop-band regions of the mono-coupled periodic system are analytically determined for period-q orbits as they are governed by the eigenvalues of the linearized 2D map arising from linear stability analysis of periodic orbits. Then, equivalent chains of nonlinear oscillators in complex domain are tackled. Also in this case, where a 4D real map governs the wave transmission, the nonlinear pass-and stop-bands for periodic orbits are analytically determined by extending the 2D map analysis. The analytical ﬁndings concerning the propagation properties are then compared with numerical results obtained through nonlinear map iteration.


Introduction
One-dimensional chains of coupled oscillators are suitable for modeling a number of physical systems arising in different scientific contexts such as condensed-matter physics, optics, chemistry and mechanics.In the latter, periodically reinforced structures as well as trusses, are frequently used in aerospace, naval, and civil engineering and can be conveniently modeled as periodic systems.It is well known that the dynamics of linear multi-coupled periodic systems is governed by frequency intervals or bands where disturbances propagate harmonically (pass-bands), decay (stop-bands) or propagate harmonically with attenuation (complex-bands).Moreover, the natural frequencies fall within pass-bands.In contrast to the vast literature concerning linear periodic structures, few works have so far appeared in the mechanical context on the dynamics of nonlinear periodic systems.Monocoupled periodic systems of infinite extent with material nonlinearities have been addressed in [8].Two different asymptotic approaches have been devised for studying standing (stop-band) and traveling (pass-bands) waves; amplitude dependent frequencies bounding nonlinear propagation and attenuation zones have been found.In [9] an array of elastic oscillators coupled through buckling sensitive elastica has been addressed both numerically and experimentally.The existence of transitions from soliton-like motions to spatially and temporally disordered motions due to a sudden excitation has been shown relying on a modified Toda lattice model.Chains of oscillators with cubic nonlinearities have been studied in [7,3] through numerical and asymptotic approaches, respectively.
Being interested in dynamical phenomena with the same length scale of the interoscillator distance, discrete nonlinear systems are considered in this work, whose main goal is to analytically investigate the modification of the boundary of the linear propagation/attenuation zones due to the nonlinearities.In particular, one-dimensional chains of linearly coupled nonlinear oscillators are addressed by means of nonlinear maps [1,2].In this realm, the governing difference equations are regarded as symplectic nonlinear transformations relating the amplitudes in adjacent chain sites (n, n + 1), thereby considering a dynamical system where the location index n plays the role of the discrete time.By doing so, wave propagation becomes synonymous of stability: finding regions of propagating wave solutions (passbands) is equivalent to finding regions of linearly stable map solutions.
The paper is organised as follows.In Section 2 monocoupled chains of oscillators with cubic nonlinearity are studied.Pass-and stop-band regions are analytically determined for period-q orbits as they are governed by the eigenvalues of the linearized 2D map arising from linear stability analysis of periodic orbits.In Section 3 chains of nonlinear oscillators in complex domain are tackled.The interest in this model lies in its mathematical representation as it is equivalent to a fictitious bi-coupled chain of nonlinear oscillators.In this case, where a 4D real map governs the wave transmission, the nonlinear pass-and stop-bands for periodic orbits are analytically determined by extending the 2D map analysis.In essence, the nonlinear map, embedding the linear one whose propagation zones boundaries have been derived in closed form for bi-coupled periodic systems in [4], allows to determine the corrections of such boundaries due to the nonlinearities.In both 2D and 4D cases the nonlinearity causes amplitude dependent pass-and stopbands, and the occurrence of vibration modes in the linear attenuation zones observed in the literature [5] can thereby be explained.
Numerical investigations carried out by nonlinear map iteration are presented in Section 4. They are mainly aimed at characterizing the bounded solutions occurring within the passing band, where, besides periodic orbits, quasiperiodic and chaotic orbits do exist.Good agreement between the approximate analytical prediction of the propagation regions and the numerical evidence is found.

Chain of nonlinear oscillators in real domain
A mechanical model for a chain of linearly coupled nonlinear oscillators, schematically depicted in Fig. 1, has been chosen in the form: The equation of motion ( 1) is characterized by on-site cubic nonlinearity, as e.g. in [3,7], describing a chain of Hamiltonian oscillators.Moreover it is worth noticing that Eq. ( 1) is formally equivalent to the discrete nonlinear Schro ¨dinger equation (DNLS) in real domain considered in [1].Periodic solutions of Eq. ( 1) are sought for by assuming the time harmonic solution q n = a n cos(xt) (harmonic balance with only the first harmonic), and setting the following second-order difference equation for the stationary amplitude is obtained Eq. ( 3), relating the amplitudes a in adjacent chain sites n À 1, n and n + 1, can be rewritten in matrix form as or a n+1 = T(a n )a n , where T(a n ) is the nonlinear transfer matrix.In the linear case, T represents a symplectic linear transformation and its reciprocal eigenvalues k satisfy k 1 k 2 = 1.As well known [6], such eigenvalues govern the stationary wave transmission properties: if the eigenvalues lie on the unit circle, then free waves propagate harmonically without attenuation (pass-band, P); if the eigenvalues are real, then free waves decay without oscillations (stop-band, S).In the more general nonlinear case, T(a n ) belongs to the class of area preserving maps [1] such that det(DT(a n )) = 1, where DT is the Jacobian or tangent map with reciprocal eigenvalues.Therefore, in order to study the stationary wave transmission properties of the one-dimensional nonlinear chain (1) of length N, it is convenient to rely on the eigenvalues of the linearized map equations in the neighborhood of an orbit ranging from (a 0 , a 1 ) to (a NÀ1 , a N ).Indeed, according to [2], the transformation (4) can be considered as a dynamical system where the chain index n plays the role of discrete time n, so that the analysis of the transmission properties is equivalent to the stability analysis of the orbits.
The linear stability of a given orbit is investigated by introducing a small perturbation u n ; linearizing the map equations, a second-order difference equation for the perturbation is obtained Eq. ( 5) can be put in matrix form leading to the two-dimensional Jacobian The interest lies in the linear stability of spatially periodic orbits a n+q = a n with cycle length q.The eigenvalues of DT are determined by its trace, so the stability of period-q orbits is described by tr(DT q ), where DT q is given by the matrix product If the eigenvalues lie on the unit circle, then stable elliptic periodic cycles or oscillating solutions (pass-band, P) occur; if the eigenvalues are real, then unstable hyperbolic periodic cycles or exponentially increasing solutions (stop-band, S) occur [2], see Fig. 2.
By setting a nþ1 ¼ xnþ1 and a n ¼ ỹnþ1 , we can express the map (4) as For period-1 orbits the curves bounding the propagation regions, where the eigenvalues lie on the unit circle, can be determined by the condition jtr(DT)j and are shown in Fig. 3a in the x-positive half-plane.The curves r and s represent hyperbolic (k 1 = k 2 = 1) and reflection hyperbolic (k 1 = k 2 = À1) boundaries, respectively.In Fig. 3b-d further curves t i , lying inside the pass-band region, are depicted; they are determined by satisfying the condition jtr(DT q )j = 2 for q = 2, 3, 4 and their number depends on the periodicity of the orbit, i.e. t i with i = 1,. . ., q À 1.While the curves r are always hyperbolic boundaries, the curves s are either hyperbolic with reflection boundaries, for q odd, or hyperbolic boundaries, for q even.Whenever a period-q orbit crosses a curve t i it temporarily loses its stability or, equivalently, does not propagate, through either a saddlenode or a period-doubling bifurcation for i even or i odd, respectively.Such alternating pattern of the internal curves t i implies that, while the elliptic even-period orbits become hyperbolic at the upper pass-band boundary, the odd-period ones, at the same boundary, become hyperbolic with reflection.It must be emphasized that the periodicity q governs the new born internal curves without affecting the overall propagation region.As expected, the nonlinearity (b 5 0) implies a propagation region depending on the amplitude of oscillations , in contrast to the linear case, where the propagation region is given by jaj 6 2: this is better shown in terms of the dimensional amplitude x in Fig. 4.
In order to find periodic orbits of the map T defined by Eq. ( 4), it is convenient to exploit its reversibility, therefore it can be cast into the product of two involutions T = T 1 T 2 , given by The involutions (10) satisfy T 2 1 ¼ T 2 2 ¼ I so that the inverse of T can be expressed as The sets of fixed points of T 1 and T 2 form two curves in the phase plane, x ¼ ỹ and ỹ ¼ 1=2xðbx 2 þ aÞ respectively, called symmetry lines.The latter curves are very useful to determine periodic orbits since a periodic orbit initially on the symmetry line stays on the line as the parameters are changed [1].For instance, period-1 and period-4 orbits, representing the fixed points of T and T 4 , respectively, born on the symmetry line of T 1 , are determined as Fig. 3. Propagation region of period-q orbits: (a) q = 1, (b) q = 2, (c) q = 3, (d) q = 4.
Period-1 : and are superimposed in Fig. 3a and d, respectively.The propagation properties of the periodic orbits (11) will be numerically investigated in Section 4. It can be noticed that, regardless of the periodicity q, the bounded orbits region coincides with that of the period-1 case whose boundaries are given by the curves of Eq. ( 9).However, the loss of stability through the upper bound involves different bifurcations if even-period or odd-period orbits are considered.Namely, while odd-period ones lose stability via period-doubling bifurcation, such as point A at a = À4 in the following Fig. 7 for the period-1 case, the even-period orbits lose stability via saddle-node bifurcation, such as point B at a = À1 for the period-4 case.

Chain of nonlinear oscillators in complex domain
In this section the analysis is focused on a chain of nonlinear oscillators in complex domain.The dynamics is governed by the following differential equation which describes a number of physical models including molecular crystal chains, coupled nonlinear waveguides in optics and nonlinear electrical lattices.By setting and by substituting q n = A n e Àixt , with A complex, the following equation for the stationary solutions is obtained Eq. ( 14) can be rewritten in the form Eq. ( 15) can be put in the compact form A n+1 = T(A n )A n , where T(A n ) is still an area preserving map enabling to analyse the transmission properties along the line followed in the previous section for the real domain case.However in this case A is a complex amplitude leading to a four-dimensional real map.Thus, by introducing a small complexvalued perturbation u n and setting A n = B n + iC n , u n = x n + iy n , a 4D real map can be cast: Fig. 4. Amplitude dependent propagation band vs increasing nonlinearity. where By rewriting Eq. ( 16) as the stability of spatially periodic orbits A n+q = A n is governed by the reciprocal eigenvalues of the corresponding linearized map DT q which transfers the state through the complete cycle of length q.If the eigenvalues of DT q lie on the unit circle, then stable elliptic periodic cycles or oscillating solutions (pass-band) occur.If the eigenvalues are real, then unstable hyperbolic periodic cycles or exponentially increasing solutions (stop-band) occur.Since each member of a periodic orbit family exhibits the same stability type, it is sufficient to consider only one of the periodic points of each family [2], e.g.period-1 orbits and the associated map DT(A n ).The characteristic polynomial of the matrix DT(A n ) has the form with Although the original problem is mono-coupled and implies two-dimensional maps, a fourth order characteristic equation is obtained due to the adopted complex-valued representation.Such instance prompts the mathematical analogy with [4], where linear bi-coupled periodic structures have been analysed on the basis of the invariants I j (j = 1, 2) of transfer matrix and analytical results on the single unit free waves (characteristic waves) propagation properties have been derived.According to that paper, an exhaustive geometrical representation of such propagation properties is achieved by identifying in the two-dimensional space {I j } the domains in which the four eigenvalues k are of the same type.In particular, perfect transmission is allowed in the passing domain where the following condition holds The invariants I j (j = 1, 2) in Eqs. ( 20) and ( 21) of the four-dimensional map in tangent space DT(A n ) are amplitudedependent, therefore, from (22) it follows that transmission is possible within the following boundaries Eq. ( 23) represent analytical approximations to the pass-band boundaries modified by the presence of nonlinearities (see Figs. 5 and 6).The curves r and s represent hyperbolic and reflection hyperbolic boundaries, respectively; regardless of the level of nonlinearity b, their intersection occurs at a ¼ ffiffi ffi 6 p .It must be noted that the same condition for the corresponding linear model gives j2 À a 2 j 6 2, in agreement with the linear dispersion relation that can be obtained by considering a periodic traveling wave solution of the linearized equation (12) in the form q n (t) = Ae i(jsnÀxt) , leading to 2 À a 2 = 2 cos(js), where j and s represent the wave number and the spacing between adjacent masses, respectively.According to the actual mono-coupled nature of the system, complex conjugate eigenvalues cannot occur as the condition for their existence, 8 þ I 2 1 À 4I 2 < 0, leading to 4b 2 jA n j 4 < 0, cannot be satisfied.
For period-2 orbits, Eq. ( 22) must be satisfied by the invariants I j (j = 1, 2) of the map DT 2 (A n ).In this case the external boundaries of the pass-band are still given by Eq. ( 23), while two new born curves, shown in Fig. 5a, arise inside the P region, namely t 11 :¼ ðjA n j; aÞ j jA n j ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ða 2 À 2Þ=b p and t 12 :¼ ðjA n j; aÞ j jA n j ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ða 2 À 2Þ=3b p .Moreover, the boundaries r and s become both hyperbolic, while the new internal ones are both of reflection hyperbolic type.Whenever a period-1 orbit crosses one of these curves, it temporarily loses its stability through a period-doubling bifurcation.

Real domain
Periodic orbits of map T in Eq. ( 8) are numerically investigated in this section.In the numerical procedure only the orbits with amplitude smaller than a selected upper limit are considered bounded; it has been verified that the stability or boundedness zone boundary is not affected by the upper limit value as long as the number of iterations exceeds the order of 10 2 .In Fig. 7 bifurcation diagrams for period-1 and period-4 orbits obtained by iterating the map T starting with initial conditions along the curves of Eq. ( 11), are shown.According to the analytical predictions, bounded periodic orbits are found within the P zone as indicated by the continuous lines, whereas after crossing the upper boundary unbounded orbits take place, as shown by the occurrence of empty intervals along the curves.The eigenvalues obtained by the numerical iteration of the maps confirm for odd-and even-period orbits the different type of loss of stability analytically predicted.Yet, it is worth noticing that, within the propagation zone, not only regular solutions but also a rich variety of nonregular bounded orbits can occur [1].In Fig. 8 several orbits occurring numerically in the neighbourhood of the intersections of period-1 and period-4 orbits with internal critical thresholds, when varying initial conditions, are shown on the phase plane of T. Temporary loss of stability of period-1 orbits according to a tripling bifurcation at a = À3.5 and a quadrupling bifurcation at a = À3 are shown in Fig. 8a and b, respectively.These bifurcated solutions occur where the period-1 orbits cross the curves t 2 of T 3 (point A in Fig. 8a) and T 4 (point B in Fig. 8b), respectively, where saddle-node bifurcations are analytically predicted.In particular, Fig. 8a shows the fixed point of the map T (period-1 orbits) on the symmetry line at x = y = 1.224 splitting into a period-3 orbit, as well as the ensuing quasiperiodic orbits diverging along three paths when varying the initial conditions.In Fig. 8b the fixed point at x = y = 1 splits into a period-4 orbit and then quasiperiodic orbits and chaotic layers occur before the onset of unbounded solutions.In turn, Fig. 8c, besides a number of neighboring quasiperiodic orbits, shows two pairs of doubled numerical solutions at a ¼ À ffiffi ffi 2 p =2, consistent with the analytical crossing of the reference period-4 orbit with the period doubling bifurcation curve t 3 of T 4 (point C in Fig. 8c).It must be noticed that since period-2 orbits do not exist, period doubling bifurcations of period-1 orbits at the intersection with t 1 of either T 2 or T 4 cannot occur.

Complex domain
To the aim of numerically investigating the periodic orbits of Eq. ( 12), the stationary equation ( 14) can be recast using polar coordinates by substituting A n ¼ r n e i#n and separating real and imaginary parts, thus leading to where D# n+1 = # n+1 À # n .Next, by introducing the following change of variables, representing a bilinear combination of the wave amplitudes on each segment of the chain [2], and defining J = r n r nÀ1 sin D# n , after rescaling the quantities x n ¼ 2Jx n , z n ¼ 2Jz n and b ¼ J b, Eq. ( 24) can be written, omitting the bar, as a two-dimensional real map T : The numerical analysis will be focused on period-1 and period-2 orbits of the map (26), which are determined by [2] Period-1 : . In order to analyse the results of the numerical investigations against the propagation properties of the original chain of complex nonlinear oscillators (12), the orbits in Eq. ( 27) are mapped onto the a À jAj plane where the stability boundaries (23) have already been determined.By considering that r n = jA 0 j and setting r 1 = r 0 = jA 0 j, the mapping can be obtained by the following relation ensuing from Eq. ( 25) on account of the J-rescaled variables.The rescaling parameter J, with meaning of a conserved physical quantity, is determined as follows.The period-1 orbit equation (27 1 ) for b > 0 has either no roots or two real roots corresponding to one hyperbolic and one elliptic fixed point.Therefore, it can be inferred that it crosses the hyperbolic boundary given by Eq. (23 1 ) with b ¼ b=J , where the two real roots are born.By labeling the coordinates of the intersection point with ða Ã ; x Ã 0 Þ, for a given strength of nonlinearity b, J can be obtained by matching Eqs. ( 28) and (29) thus leading to With a path following procedure, the coordinates of each intersection point, corresponding to different values of the conserved quantity J according to (30), are then used as starting points to obtain bifurcation diagrams of period-1 orbits through numerical iteration of the map (27 1 ).In Fig. 9a the bounded period-1 orbits, represented by dotted lines, are shown for five different values of J and b = 0.2.Within the parameter range investigated, 0 < J < 2, only bounded orbits are found inside the P zone, consistent with the analytical prediction.The bifurcation diagram for period-2 orbits, obtained through iteration of the map (27 2 ), is also shown in Fig. 9b.For each value of J, these orbits, represented by the thicker dotted lines, are born at the period doubling bifurcations of period-1 orbits; according to the analytical predictions, such bifurcation points are located along the internal line t 11 ðjA n j ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ða 2 À 2Þ=b p Þ, whereas the line t 12 seems to play no role.
In Fig. 10a a number of orbits obtained numerically from the map (26) when varying either a control parameter or the initial conditions are shown on the phase plane for b = 0.2 and J = 1.Fig. 9a shows at a = 1.111 period-1 orbits and then, as a increases, quasiperiodic orbits; afterwards, at a = 1.68, period-2 orbits surrounded by quasiperiodic and chaotic orbits take place.A similar involved picture of regular and complex orbits is found by varying initial conditions at a = 1.63 as shown in Fig. 10b where period-2 orbits, surrounding quasiperiodic orbits, and a chaotic separatrix layer are observed.Finally, Fig. 10c shows several period-1 orbits (represented by small circles) and a chaotic sea developed after the period doubling bifurcation occurring at a ¼ ffiffiffiffiffiffiffiffiffiffiffi 2 þ b p .

Conclusion
Propagation properties of periodic orbits in nonlinear oscillatory chains have been studied based on the analogy with the linear stability analysis of the relevant maps (space/time, i.e. a nonlinear dynamical approach to wave propagation properties).The boundaries of pass-band regions have been analytically determined for both real and complex oscillatory chains with cubic nonlinearities.Within pass-band regions bifurcation thresholds of various order periodic orbits have also been determined.Regardless of the periodicity q, the bounded orbits region coincides with that of the period-1 case.However, the loss of stability, or orbitsÕ unboundedness, occurring at the upper bound involves different bifurcations if even-period or odd-period orbits are considered, namely, saddle-node and period-doubling bifurcation, respectively.The analytical predictions have been validated through comparison with numerical solutions of the maps highlighting the occurrence of a variety of bounded, also nonperiodic, orbits.

= 2 .
Fig. 2. Stability-propagation equivalence based on the eigenvalues of the Jacobian DT. r

Fig. 10 .
Fig. 10.Numerical phase portraits (right), for b = 0.2 and J = 1, with a varying control parameter (a, c) or varying initial conditions (b), within the analytical propagation band (left).