Yu V Mikhlin

T V Shmatko

G V Manucharyan

LYAPUNOV DEFINITION AND STABILITY OF REGULAR OR CHAOTIC VIBRATION MODES IN SYSTEMS WITH SEVERAL EQUILIBRIUM POSITIONS

This paper deals with forced vibrations of two-DOF systems with more than one equilibrium positions. Such systems may be obtained by digitization of elastic post-buckling systems. A vibration mode, which is periodic at small force amplitudes and becomes chaotic as the force amplitudes are slowly increased, is selected. It is possible to formulate and solve the problem of stability of a periodic or chaotic vibration mode in a space with greater dimension using the classical Lyapunov stability definition and some calculating procedures. Instability of phase trajectories is used as a criterion of the chaotic behavior in dynamical systems. Trajectories with very close initial values are compared. Use of the Lyapunov stability definition shows mutual stability/instability of the trajectories. Calculations permit to observe an appearance and enlargement of the chaotic behavior regions. Specific results are obtained for the nonautonomous Duffing equation and pendulum system.

Introduction

Models under consideration may be obtained by digitization of some elastic systems that have lost stability under a constant compressive force. If displacements of the nonlinear elastic system are approximated by a single harmonic of the Fourier series expansion for space coordinates, a system having a single degree of freedom is obtained. Behavior of the model described by the nonautonomous Duffing equation was examined in papers [1,2 et al.]. At small force amplitudes, periodic motions close to one of two stable equilibrium positions are observed but irregular, apparently chaotic motions begin as the force amplitudes are slowly increased.

If the first two harmonics of the Fourier series for space coordinates are used, one obtains a set of two-second order ODEs, coupled in nonlinear terms only. In this case an energy 1 "pumping" from one vibration mode to another mode is possible. Thus one can formulate a problem of the stability of periodic or chaotic vibration mode in a space with a greater dimension. There are many publications on the general and special problems of stability of periodical solutions, e.g., fundamental works by Lyapunov [START_REF] Liapunov | Stability of Motion[END_REF]; Ince [START_REF] Ince | Ordinary Differential Equations[END_REF]; Andronov, Vitt and Khaikin [START_REF] Andronov | Theory of Oscillations[END_REF]; Minorsky [START_REF] Minorsky | Nonlinear Oscillations[END_REF] etc. But a problem of stability of chaotic modes in a space with a greater dimension has no analytical solutions in a general case.

We make use of the classical Lyapunov stability criterion for a case when initial variations are not arbitrary small that is they are limited below (Section 3). As a result we obtain some test suitable for a computer realization. It is discussed a connection of this test with an approach used by Schiehlen [START_REF] Schiehlen | Nonlinear oscillations in multibody systems[END_REF], and with so-called Lyapunov exponents(Section 3). Calculations, corresponding to the proposed test, are realized as long as boundaries of stability and instability regions in the space of system parameters are variable. Stability of regular or chaotic modes of nonlinear beams and shells (Section 4) are considered. Note that first results obtained for this problem were published in papers [START_REF] Yu | Stability of periodical and chaotic vibrations in systems with more than one equilibrium positions[END_REF][START_REF] Yu | Stability of periodic and chaotic modes in systems with several equilibrium positions[END_REF].

The chaotic behavior of nonlinear dynamical systems is an object of recent and current investigations. The first ideas on the chaotic behavior of dynamical systems can be found in works by Poincare' [START_REF] Poincare | Les Methods Nouvelles de la Mecanique Celeste[END_REF] and Birkhoff [11 et al.].

Criteria for the onset of chaos beginning in dynamical systems are very diverse. A detailed survey of various such criteria of the chaos beginning was published by Moon [START_REF] Moon | Chaotic Vibrations[END_REF]. The most important analytical or numerical-analytical criterions of the chaos beginning: criteria are based on the configuration of the homoclinic trajectories, period doubling, overlapping of resonances, analysis of spectra, calculation of Lyapunov exponents, "Smale horseshoe" et al [START_REF] Lorenz | Deterministtic non-periodic flow[END_REF][START_REF] Ruelle | On the nature of turbulence[END_REF][START_REF] Melnikov | On the stability of the centre for time periodic perturbations[END_REF][START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems,and Bifurcations of Vector Fields[END_REF]. All trajectories contained in the "strange attractor" region are unstable, that is trajectories which are close to each other at some point of time then diverge exponentially [START_REF] Moon | Chaotic Vibrations[END_REF][START_REF] Ueda | Randomly transitional phenomena in the system governed by Duffing's equation[END_REF][START_REF] Wolf | Quantifying chaos with Lyapunov exponents[END_REF].

We suggest here an approach to determine the onset of chaos based on the Lyapunov definition and the proposed calculating test (Section 5). Mutual instability of phase trajectories is used as criterion of chaotic behavior in dynamical systems. One compares trajectories that are initially very close. The proposed test shows the mutual instability of the trajectories. Calculations permit to observe a process of appearance and enlargement of the chaotic behavior regions if some selected parameters of the dynamical system under consideration are changing. Specific results were obtained for the nonautonomous Duffing equation and pendulum system when an amplitude of the external periodic action increases.

Principal Model

Consider the following system that can be obtained by discretization of equations of nonlinear dynamics of some elastic systems:

ÿ1 + δ ẏ1 -α y 1 + β y 3 1 + c y 1 y 2 2 = f cos ω t, ÿ2 + δ ẏ2 + ay 2 + b y 3 2 + c y 2 y 2 1 = 0, (1)
where y 1 (t) and y 2 (t) are unknown functions; δ is the coefficient determining friction; all coefficients are positive (excepting may be the coefficient α).

The system (1) can be obtained in the following problem. Consider nonlinear bending vibrations of a beam with a length equal to l. Within framework of the Kirchhoff hypothesis the vibrations are described by the following equation [START_REF] Kauderer | Nichtlinear Mechanik[END_REF]:

µ ∂ 2 w ∂ t 2 + EI ∂ 4 w ∂ x 4 - ES 2l l 0   ∂w ∂ x 2 dx   ∂ 2 w ∂ x 2 + Γ ∂ 2 w ∂ x 2 = F (x, t), ( 2
)
where x is the space coordinate; w(x, t) is the bending deflection; µ is the density of the beam material; E is the Young's modulus; S is the area of the transversal cross-section; I is the corresponding moment of inertia; Γ is the axial compressive force; friction will be introduced later; F (x, t) is the distributed external effect. Let w = w xx = 0 if x = 0 and x = l. Let the function F (x, t) be presented as F (x, t) = f 1 cos ωt sin π x l . Representing the deflection as truncated Fourier series (by sines) and making use of the Bubnov-Galerkin procedure for [START_REF] Moon | Chaotic Vibrations[END_REF], one comes to a system of nonlinear ODEs connecting coefficients of the series. Taking into account two harmonics of the Fourier series expansion,

w = y 1 (t) sin π x l + y 2 (t) sin nπ x l , (3)
one obtains the system (1) for determination of y 1 (t) and y 2 (t), where δ is the coefficient determining friction;

α = -1 µ EI π l 4 -π l 2 Γ ; β = ES 4µ l π l 4 ; c = ES 4µ l π l 4 n 2 ; a = 1 µ EI nπ l 4 -nπ l 2 Γ ; b = ES 4µ l nπ l 4 ; f = 1 µ f 1 .
In the post-buckling dynamics we have Γ > π 2 EI l 2 , and the coefficient α > 0 too. If y 2 = 0, one obtains a nonautonomous Duffing equation with two potential wells (for a case α > 0). It was shown [START_REF] Holmes | A nonlinear oscillator with a strange attractor[END_REF][START_REF] Moon | Chaotic Vibrations[END_REF] that as the force amplitude f increases, the behavior of the system changes. A sequential doubling of the vibration period results, and beginning with some values of f it is possible to observe irregular, apparently chaotic motions (for example, if α = 10, β = 100, δ = 1, ω = 3.76, the chaotic vibrations appear if the force amplitude f is a little greater than 1). Additional analyses (spectral, and also based on Poincare' mapping) show that a region of the phase space exhibits the properties of a "strange attractor".

Equations which are similar to the system (1) can be obtained also in problem of prebuckling or post-buckling nonlinear dynamics of plates or shells.

Lyapunov stability definition and its computer realization

One consider a problem of stability of the periodic or chaotic vibration mode y 2 = 0 of the system (1). Instability of the mode means energy "pumping" from one harmonic of the Fourier series to another. The variables y 2 , ẏ2 emerge in the form of variations (or perturbations). Assume first that values of the variations are significantly smaller than values of the variable y 1 (in the region of stability of the mode y 2 = 0). This corresponds to a conventional formulation of the problem of local stability in configuration space.

Consider the well-known Lyapunov stability criterion [START_REF] Liapunov | Stability of Motion[END_REF][START_REF] Minorsky | Nonlinear Oscillations[END_REF][START_REF] Perko | Differential Equations and Dynamical Systems[END_REF] stating that: the solution y = 0 is stable if for all ε > 0 there exists a δ > 0 such that for all y 0 ∈ N δ (0) and t ≥ 0 we have y (t) ∈ N ε (0). Here N (0) δ and N (0) ε mean δ-and ε-neighborhoods of the mode y = 0. For example, N γ = {y ∈ R n |y| ≤ γ } (the neighborhoods could be calculated differently).

Introduce a relation between the quantity ε and the initial value of the variable y. Let

ε = ρ |y 0 | ≤ ρδ ( ρ = const) ( 4 )
This means that a value of δ is not arbitrarily small because δ ≥ ε ρ . Note that this assumption does not contradict to the Lyapunov definition meaning because in this definition the initial values of δ can not tend to zero. One has from (4) that ρ ≥ ε δ , that is ρ is a high limit of the fraction ε δ . Besides, one obtains from the Lyapunov stability criterion taking into account the inequality (4) that |y (t)| ≤ ρ |y 0 |.

One has from the preceding the following stability test for the system (1): Instability of the mode

y 2 = 0 is established if max t∈[0,T ] |y 2 (t)| ≥ ρ |y 2 (0)| (5)
The inequality [START_REF] Andronov | Theory of Oscillations[END_REF] gives us a possibility to use a computer in the next stability analysis.

It is necessary to choose values of ρ and T . Here a value ρ -1 is a measure of smallness of initial variations with respect to maximum admissible variations. An increase of the value ρ means that feasible initial values of variation decrease. There is some arbitrariness in a choice of the value ρ. In engineering one says that some quantity α is greater than another quantity β "in order of magnitude" if α β ≈ 10. We can choose, for example, ρ ≈ 10 which constitutes as an "engineering test". The discussion on the choice of the constant ρ is continued in the Subsection 3.2.

Calculations are conducted as long as boundaries of stability and instability regions (in a chosen scale) on the system parameter space are variable. This is the principal criterion for the choice of the calculation time T . Note that all calculations are realized by using the standard Runge-Kutta procedure. The calculations are made at points on some chosen mesh of the system parameter space. It is clear that if the mesh widths decrease and the number of mesh points increases infinitely, the interval of time T tends to infinity.

We now discuss the dependence (or independence) of the stability analysis on the variations initial conditions. The linear stability results are not dependent on initial conditions, this is a property of linear variational equations. But it is known [START_REF] Siegel | Lectures on Celestical Mechanics[END_REF][START_REF] Pecelli | Normal modes, incoupling, and stability for a class of nonlinear oscillators[END_REF] that additional nonlinear instability regions (obtained if we take into account nonlinear terms) have a smaller dimension in parameter space than instability regions obtained by the linearized stability analysis. Examples verify that the stability analysis is independent of initial variations if the initial variations are small.

Schiehlen's analysis of stability

In 1993 W. Schiehlen [START_REF] Schiehlen | Nonlinear oscillations in multibody systems[END_REF] proposed an approach for the stability assessment of nonlinear oscillations "from an engineering point of view". His approach was based on the Lyapunov stability definition. Schiehlen introduced the "inverse stability measure"

IS1(x 0 , t 0 ) =        x(t) ∞ x 0 for x 0 = 0 1 for x 0 = 0 (6)
where

x(t) = max 1≤i≤n |x i (t)| , x(t) ∞ = lim T →∞ x(t) T = lim T →∞ max t∈[t 0 ,T ]
x(t) .

The measure IS1(x 0 , t 0 ) "characterizes the ration between a given initial state and the corresponding maximal displacement of the trajectory". Here x(t) is the time-dependent n-state vector.

The inverse stability measure

IS2(r, t 0 ) = max x 0 ∈{x: x =r} IS1(x 0 , t 0 )
"characterizes the maximal displacement of all trajectories starting out of the initial conditions' subspace which is by definition a hyper-cube with respect to the equilibrium point x = 0". These measures were used in a numerical calculations where the integration interval is limited, 0 ≤ t ≤ T . So, the measure (4) was replaced by

IS1 T (x 0 , t 0 ) = x(t) T x 0 .
There is a relation between the stability definition and the inverse stability measure. In the unstable case it is possible to choose the following series of initial conditions: x 0 1 > x 0 2 > ... > x 0n . One has lim n→∞ x 0n = 0 and, as an result, the IS1 ∞ (x 0n , t 0 ) → ∞. But if the inverse stability measure is limited for all x 0 , then x(t) ∞ ≤ s x 0 . This means that the system is stable by Lyapunov.

The approach was used by Schiehlen to analyze a stability of the single and double pendulums. It is possible to select both stability and instability regions in place of the systems initial values. But it should be noticed that the choice of the calculation time T is not discussed in paper [START_REF] Schiehlen | Nonlinear oscillations in multibody systems[END_REF].

The stability test (4) and Lyapunov exponents

"Lyapunov characteristic exponents" are introduced in some problems of stability or analysis of the "strange attractors" [START_REF] Moon | Chaotic Vibrations[END_REF][START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems,and Bifurcations of Vector Fields[END_REF][START_REF] Wolf | Quantifying chaos with Lyapunov exponents[END_REF][START_REF] Parker | Practical numerical algorithms for chaotic systems[END_REF][START_REF] Rugonyi | An evaluation of the Lyapunov characteristic exponent of chaotic continuous system[END_REF]. Positiveness of the exponents means the global mutual instability of phase trajectories; this is accepted as a criterion of the chaotic behavior onset in the nonlinear system under consideration. The largest Lyapunov exponent can be defined, for instance, as the following:

λ = lim t→∞ 1 t ln |y(t)| |y(0)| . ( 7
)
One has from here

|y(t)| |y(0)| = e (λ+∆)t , where ∆ → 0 for t → ∞ (8)
Comparing ( 8) and ( 5), one obtains the next inequalities

ρ ≤ e (λ+∆)T , T ≥ ln ρ λ + ∆ (λ > 0), ( 9
)
where T is the maximal calculation time. Thus the constants ρ and T presented in the criterion (5) are linked by the inequality [START_REF] Yu | Stability of periodic and chaotic modes in systems with several equilibrium positions[END_REF]. It is clear that a decrease of the value ρ gives the corresponding decrease of the calculation time T (it is confirmed by concrete computations too). Near the boundaries of the stability/instability regions values of the Lyapunov exponent λ decrease, and the calculation time increase. But the value ρ can not be chosen very small because near the boundaries values of λ are small, and the calculation time is depended on the uncertain value ∆. The significant arbitrariness in a choice of the value ρ will be showed in the next section, too. We note that in the work by Rugonyi and Bathe [START_REF] Rugonyi | An evaluation of the Lyapunov characteristic exponent of chaotic continuous system[END_REF] the concept of Lyapunov stability is connected with the Lyapunov characteristic exponents calculation. The effect of velocities (see the next Remark here) in the exponents calculation is analyzed too. Besides, the procedure proposed in [START_REF] Rugonyi | An evaluation of the Lyapunov characteristic exponent of chaotic continuous system[END_REF] is used in the stability investigation of some continuous nonlinear systems, namely, the buckled beam and fluidstructure interaction systems. Additional discussion on the proposed approach and the Lyapunov exponents calculation approach is contained in the section Discussion. REMARK. We can use in place of the test (5) the following "energetic" test. Namely, instability is fixed if

y 2 2 (t) + ẏ2 2 (t) ≤ ρ y 2 2 (0) + ẏ2 2 (0) (0 ≤ t ≤ T ) . ( 10
)
It is clear, and checking calculations confirm it, that regions of instability in the system parameter space, obtained if we use [START_REF] Poincare | Les Methods Nouvelles de la Mecanique Celeste[END_REF] are greater that those obtained if we use test [START_REF] Andronov | Theory of Oscillations[END_REF]. This follows from the fact that there are regions where the variations y 2 (t) are small but the corresponding velocities ẏ2 (t) are large.

Example

We show now the significant arbitrariness in a choice of the value ρ. For definiteness we consider the following system which can be obtained from (1) if we use some values of the physical and geometrical characteristics of the corresponding elastic system and use some scale transformations:

..

y 1 + δ . y 1 + [ 1 -q + y 2 1 + 4ε y 2 2 ] y 1 = f cos ωt, .. y 2 + δ . y 2 + 4 [4 -q + ε y 2 1 + 4y 2 2 ] y 2 = 0. ( 11
)
Consider the stability of the mode y 2 = 0, which can be both regular and chaotic. One selects the following domain of plane of parameters (f, ε): 0 ≤ f ≤ 10, 0 ≤ ε ≤ 20. Here δ= 0.1. Mesh width in the domain was chosen as: ∆f = 1/3, ∆ε = 1.

In the chosen scale of mesh, boundaries of stability and instability regions for the chosen region of the parameters plane stabilized in the interval 0 ≤ t ≤ T (T ≈ 250).

Increasing a value of the constant ρ presented in the test (5), we consider results of the calculations for ρ = 10; 5; 2.5 . It is important that the stability regions are the same for the different values of ρ and q. Corresponding results are presented in Fig. 1, where q = 0.9. A smoothing function is used in tracing of the graphic. The results show that for the system under consideration the localization of the vibration mode is the same as the stability of the mode. 11), obtained using the stability test [START_REF] Minorsky | Nonlinear Oscillations[END_REF]. Here q=0.9. Boundaries of stability and instability regions for the chosen region of parameters plane are stabilized on the interval of calculations, 0 ≤ t ≤ T (T ≈ 250). Results are the same for ρ = 10; 5; 2.5 .

Stability of Forced Vibration Modes in Nonlinear Shells

Consider equations of the nonlinear shallow-shell theory connecting the normal deflection w(x, y, t) and the stress function Φ(x, y, t) [START_REF] Th | The buckling of thin cylindrical shells under axial compression[END_REF][START_REF] Dowell | Modal equations for the nonlinear flexural vibrations of a cylindrical shell[END_REF][START_REF]Thin-Shell Structures. Theory, Experiment and Design[END_REF][START_REF] Volmir | Nonlinear Dynamics of Plates and Shells[END_REF]:

D h ∇ 4 w = L (w, Φ) + 1 R ∂ 2 Φ ∂ x 2 - ρ g ∂ 2 w ∂ t 2 + q, 1 E ∇ 4 Φ = -(1/2) L (w, w) - 1 R ∂ 2 w ∂ x 2 , ( 12
)
where x, y are respectively the longitudinal and circumferential coordinates; E is the Joung's modulus; the cylindrical stiffness D = Eh 3 12(1-ν 2 ) ; ν is the Poisson ratio; h is the shell thickness; R is the radius of the shell; ρ is the density of the shell; q(x, y, t) is the transversal loading; the operator L is the following:

L (u, v) = ∂ 2 u ∂ x 2 ∂ 2 v ∂ y 2 + ∂ 2 u ∂ y 2 ∂ 2 v ∂x 2 -2 ∂ 2 u ∂ x∂ y ∂ 2 v ∂ x∂ y .
Circular cylindrical shell complete around the circumference is considered here. Let the shell be compressed along an element of the cylinder by uniformly distributed static axial forces after the shell has lost stability. One presents here the transverse loading as q = F cos ωt sin mπx l sin ny R .

It means that the external periodical force excites the primary vibration mode with respect to the space coordinates (x, y). Vibrations in the mode may be chaotic if we consider the post-buckling case. Two variants are possible. If the mode is orbitally stable then we obtain a localization of energy. But If the mode is unstable, the transfer of energy to other space mode is possible.

Companion modes in shells: analysis of stability

Let here the deflection w be approximated by the following truncated Fourier series [START_REF]Thin-Shell Structures. Theory, Experiment and Design[END_REF][START_REF] Volmir | Nonlinear Dynamics of Plates and Shells[END_REF][START_REF] Kubenko | Nonlinear interaction of the bending vibrations modes of cylindrical shells[END_REF]:

w = f 1 cos sy sin rx + f 2 sin sy sin rx + n 2 4R (f 2 1 + f 2 2 ) sin 2 rx, ( 13
)
where s = n R , r = mπ l ; n, m are integer positive numbers; f 1 , f 2 are unknown functions of time; the transverse loading is chosen of the form q = f cos ωτ cos sy sin rx.

To analyze here a principal problem of the vibration modes energy pumping, we suppose that an initial deflection w 0 of the cylindrical shell is equal to zero. Obtaining the stress function Φ in the form of truncated Fourier series with respect to the space coordinates and using the standard Bubnov-Galerkin procedure, we obtain the equations to determine f 1 (t) and f 2 (t) which were exposed, for example, in [START_REF]Thin-Shell Structures. Theory, Experiment and Design[END_REF][START_REF] Kubenko | Nonlinear interaction of the bending vibrations modes of cylindrical shells[END_REF].

One analyzes now a stability of the vibration mode f 1 = f 1 (t), f 2 = 0. Let us the variable f 2 emerges as a small variation. After a linearization by f 2 we can obtain the following equations:

d 2 f 1 dt 2 • 1 + 2χf 2 1 + ω 2 0 f 1 + 2χf 1 df 1 dt 2 + γ 1 f 3 1 + gf 5 1 = f cos ωt, d 2 f 2 dt 2 + f 2   ω 2 0 + 2χ   df 1 dt 2 + f 1 d 2 f 1 dt 2   + γ 1 f 2 1 + gf 4 1   = 0, ( 14
)
where

ω 2 0 = 1 ρ D h s 2 + r 2 2 + Er 4 R 2 (s 2 + r 2 ) 2 , 2χ = 2 3 n 2 2R 2 , γ 1 = 1 ρ E 16 r 4 + Dn 4 r 4 hR 2 - Er 4 s 4 (s 2 + r 2 ) 2 , g = 3E 16ρ n 2 r 4 s 6 1 (s 2 + r 2 ) 2 + 1 (s 2 + 9r 2 ) 2 .
One introduces the following change of variables in equations ( 14):

x = f 1 h , y = f 2 h , τ = ω 0 t, Ω = ω ω 0 , F = f ω 2 h , K = 2χh 2 , G = γ 1 h 2 ω 2 0 . ( 15
)
We neglect later the small order terms gf 5 1 and gf 4 1 f 2 . Besides, we introduce a coefficient α which is defined by the constant static axial force. In a case of α=0 the force is absent; in a case of α=1 (this is a critical loading) a trivial equilibrium position becomes unstable; in a case of α > 1 (this is a post-buckling case) new equilibrium positions appear. Under external periodic excitation a transfer from periodic motions to irregular, apparently chaotic motions is possible. Introducing a dissipation, we can obtain, as a result, the following equations instead of the system ( 14):

x = F cos Ωτ -(1 -α) x -Kxx 2 -Gx 3 -δx / 1 + Kx 2 , y = -y (1 -α) + K x 2 + xx + Gx 2 + δy . ( 16
)
The obtained equations ( 16) permit to investigate a stability both subcritical, and supercritical forced vibrations of cylindrical shells under the constant static axial force. 16) in the plane F, Ω: 0 ≤ F ≤ 1, 0 ≤ Ω ≤ 2, obtained by using the test (6) (x(0) = 2, α = 0, δ = 0.05, T = 1000).

x (0) = 0.5, 1, 2; initial perturbations y(0) are the following: y (0) = 0.01 • x(0). It is accepted that x (0) = y (0) = 0.

Calculations which were made by using the test (5) permit to determine stability/ instability regions of the vibration mode x = x(t), y = 0 in the plane F, Ω. One considers the following domain of the plane F, Ω: 0 ≤ F ≤ 1, 0 ≤ Ω ≤ 2. Mesh width in the domain was chosen as the following: ∆F =0.05, ∆Ω =0.09. In Figs. 2345a small part of the numerical calculation results is presented. Here points of the chosen mesh which correspond to the unstable regions, are marked by symbol. Time of the stability/ instability boundaries stabilization T is pointed out too. In the instability regions only the coupled (f 2 = 0) companion modes [START_REF] Ruelle | On the nature of turbulence[END_REF], which were investigated by Kubenko et al. [START_REF] Kubenko | Nonlinear interaction of the bending vibrations modes of cylindrical shells[END_REF], are realised.

An increase of dissipation implies a corresponding decrease of the instability regions. Besides, an increase of the axial static stress coefficient α implies a corresponding increase of the instability regions. Results are the same for x(0) = 0.5 and x(0) = 1.

In a case of α > 1 a trivial equilibrium position of the shell is unstable, and two new stable equilibrium positions appear, this is a post-buckling phase of the shell behavior. Under external periodical excitation chaotic vibrations can be appeared. The stability analysis by using the test ( 4) is not depends on the type (regular or chaotic) of the vibration mode under consideration. It is interesting that for α > 1 both regular or chaotic vibration mode f 1 = f 1 (t), f 2 = 0 is unstable, that is only the coupled companion modes [START_REF] Ruelle | On the nature of turbulence[END_REF] are realised in this case.

Note that a similar analysis can be realized in a problem of stability of the forced vibration mode of the sloping arch which dynamics is described within a framework of the Kirchhoff hypothesis. A boundary of the stability/ instability regions is obtained, and the corresponding results were presented in paper [START_REF] Yu | Stability of regular or chaotic post-buckling vibration modes of elastic shells, Nonlinear Dynamics of Shells and Plates[END_REF].

Mutual Instability of Phase Trajectories

Consider the well-known nonautonomous Dufffing equation, ÿ + δ ẏα y + β y 3 = f cos ω t [START_REF] Wolf | Quantifying chaos with Lyapunov exponents[END_REF] where the system parameters are: δ = 1, α = 10, β = 100, ω = 3.76 , and the value of the amplitude f is variable. Examine the following region of phase plane of the equation ( 17): 0 ≤ y ≤ 0.66; 0 ≤ ẏ ≤ 0.8. Introduce some mesh in the defined region using the increments: ∆y = 0.02, ∆ ẏ = 0.016. Points of the mesh P ij (y i0 , ẏj0 ) will be chosen as initial points for the selected phase trajectories y

(1) ij (t) of the equation [START_REF] Wolf | Quantifying chaos with Lyapunov exponents[END_REF]. Let us take other initial points near the chosen initial points P ij , namely, Q ij (y i0 + ∆y 0 , ẏj0 ) where the value ∆y 0 is sufficiently small, and consider the other phase trajectory, y

(2) ij (t). One compares trajectories outgoing from the near initial points, using the limited Lyapunov test [START_REF] Andronov | Theory of Oscillations[END_REF].

Instability of the outgoing trajectory is established if

y (1) ij (t) -y (2) ij (t) ≥ ρ | ∆y 0 | (0 ≤ t ≤ T ) . ( 18
)
We now select ρ = 10. Let ∆y 0 = 0.002. Results of the stability analysis (the time of stabilization here, T < 100) are presented in Fig. 6 (a,b,c) for different values of the external amplitude: f = 0.8; 0.856; 0.9. In Fig. 6 the initial points of the chosen mesh, which correspond to unstable trajectories, are marked by dark squares. The calculations (in the chosen mesh of the equation phase place) show that for f < 0.85 the mutual instability of phase trajectories can be observed near the separatrix branches. Instability regions begin to extend if values of the external amplitude f are slightly larger than 0.85, and the enlargement is very fast. The results generically correspond to known behavior of the nonautonomous Duffing system [START_REF] Holmes | A nonlinear oscillator with a strange attractor[END_REF][START_REF] Moon | Chaotic Vibrations[END_REF]. In Fig. 7 the results of stability analysis for the equation of the pendulum with excited point of suspension are presented. The equation of motion is the following:

ÿ + δ ẏ + (1 + f cos ωt) sin y = 0.
The calculations were made with the friction δ = 0.1, the force frequency w = 2, the time of calculations T = 50, the widths of the mesh ∆y = 0.1, ∆ ẏ = 0.1 and ∆y 0 = 0.01. ρ equals 10. Obtained results correspond to the behavior of this system described in book of Neymark, Landa [START_REF] Yu | Stochastic and Chaotic vibrations[END_REF].

Discussion

It is appropriate to compare the presented approach with the well known methods of the Lyapunov characteristic exponents calculations. We can note that it is not necessary to analyze here variational equations and calculate corresponding eigenvalues because only some direct consequence from the classical Lyapunov stability definition is used in the proposed approach. The Lyapunov exponents are introduced, for example, in [22, p.67 and further] as a limit of the eigenvalues when the independent variable tends to infinity. Of course, we can obtain a smaller information about the system behavior than from the Lyapunov exponents calculation approach, but we avoid some known computing difficulties [22, p.73 and further]. Note that contemporary procedures to calculate the Lyapunov exponents both in discrete and continuous nonlinear systems are described in [START_REF] Rugonyi | An evaluation of the Lyapunov characteristic exponent of chaotic continuous system[END_REF].

Next, we think that in numerous works on numerical approaches a problem of choice of the calculation time is not explained until now, although some discussion about a selection of the Lyapunov exponents calculation time is presented, for example, in [22, p.81]. Unfortunately, the principle "to calculate as long as possible" is dominanted in many works. It seems to us that we proposed a simple test for determination of the calculation time starting from the real calculations at points on some chosen mesh of the system under consideration parameter space.

Conclusions

We analyzed in this paper the stability of regular or chaotic vibration modes in two-DOF multistable nonlinear systems. The stability problem for chaotic modes has no analytical solution because the chaotic motions cannot be described without a computer simulation. A new stability test is proposed here. The test permits to implement on computer the classical Lyapunov stability definition for a case when initial variations are not arbitrary small, that is they are limited below . Note that a realization of this approach is not difficult because only the standard Runge-Kutta procedure with an additional stability condition is used. Calculations by computer which correspond to the test are conducted for points of some chosen mesh of the system parameter space and continued until boundaries of stability and instability regions (in the mesh scale) are stabilized. This is the principal criterion for a choice of the calculation time. Stability of regular or chaotic vibration modes of some nonlinear elastic system is considered.

The methodology presented in this work is sufficiently general to be applicable to stability problems for vibration modes in different types of nonlinear systems.

The approach of determination of the chaos onset, which is based on the Lyapunov stability definition, is also proposed. Mutual instability of phase trajectories is used as a criterion of the chaotic behavior of dynamical system. An effectiveness of this approach is confirmed by calculations of the chaos onset of the nonautonomous Duffing equation and pendulum system. The methodology can be considered as an alternative to the well-known Poincare' map approach, and can be applicable to a determining of the chaos onset in other dynamical systems.

Figure 1 :

 1 Figure 1: Stability (S) and instability (U) regions of the mode y = 0 of equations (11), obtained using the stability test[START_REF] Minorsky | Nonlinear Oscillations[END_REF]. Here q=0.9. Boundaries of stability and instability regions for the chosen region of parameters plane are stabilized on the interval of calculations, 0 ≤ t ≤ T (T ≈ 250). Results are the same for ρ = 10; 5; 2.5 .

Figure 2 :

 2 Figure 2: Instability regions of the vibration mode x = x(t), y = 0 of the equations (16) in the plane F, Ω: 0 ≤ F ≤ 1, 0 ≤ Ω ≤ 2, obtained by using the test (6) (x(0) = 0.5, α = 0.5, δ = 0.1, T = 400).

Figure 3 :

 3 Figure 3: Instability regions of the vibration mode x = x(t), y = 0 of the equations (16) in the plane F, Ω: 0 ≤ F ≤ 1, 0 ≤ Ω ≤ 2, obtained by using the test (6) (x(0) = 0.5, α = 0, δ = 0.05, T = 1000).

Figure 4 :

 4 Figure 4: Instability regions of the vibration mode x = x(t), y = 0 of the equations (16) in the plane F, Ω: 0 ≤ F ≤ 1, 0 ≤ Ω ≤ 2, obtained by using the test (6) (x(0) = 2, α = 0, δ = 0.05, T = 1000).

Figure 5 :

 5 Figure 5: Instability regions of the vibration mode x = x(t), y = 0 of the equations (16) in the plane F, Ω: 0 ≤ F ≤ 1, 0 ≤ Ω ≤ 2, obtained by using the test (6) (x(0) = 2, α = 0.5, δ = 0.05, T = 300).

Figure 6 :

 6 Figure 6: Mutual stability/instability of phase trajectories of the nonautonomous Duffing equation (17). Corresponding initial values of the trajectories are presented on the following region of phase place: 0 ≤ y ≤ 0.66; 0 ≤ ẏ ≤ 0.8. The test (18) are used to determine the mutual stability /instability of trajectories outgoing from the near initial points (the time of stabilization here, T < 100). Results of the stability analysis are obtained for f = 0.8 (figure a), f = 0.856 (figure b) and f = 0.9 (figure c). The initial points of the chosen mesh which correspond to unstable trajectories are marked by dark squares.

Figure 7 :

 7 Figure 7: Mutual stability/instability of phase trajectories of the pendulum equation. There -4 ≤ y ≤ 4; -2.4 ≤ ẏ ≤ 2.4. The test (18) are used to determine the mutual stability /instability of trajectories outgoing from the near initial points (the time of stabilization here, T = 50). Results of the stability analysis are obtained for a)f = 0.05; b)f = 0.2; c)f = 0.4; d)f = 0.6; e)f = 0.8; f )f = 1. The initial points of the chosen mesh which correspond to unstable trajectories are marked by dark dots.

Acknowledgments

The authors are very thankful to Dr. Konstantin Avramov for a very useful discussion on the paper.


Yu V Mikhlin

Normal Vibrations of a General Class of Conservative Oscillators

Keywords: Nonlinear systems, nonlinear normal modes of vibration, trajectories in configuration space, power series, internal resonances

This paper considers normal vibrations with curvilinear trajectories in a configuration space of systems which are close to systems permitting rectilinear normal modes of vibration. Analysis of trajectories of normal vibrations in the configuration space is used.

Introduction

Lyapunov [START_REF] Lyapunov | The General Problem of the Stability of Motion[END_REF] examined nonlinear finite dimensional systems of the form _ y = Ay + Ny y 2 R n N : R n ! R n N 2 C r r 1 [START_REF] Lyapunov | The General Problem of the Stability of Motion[END_REF] with an analytical first integral (where A is an n n constant matrix); he assumed that the linearized system, _ y = AY , possesses periodic solutions with natural frequencies which are not integrally related. Lyapunov proved that (1) possesses a one-parameter family of periodic solutions. To arrive at a solution, he used two approaches. One involved a power series of the amplitude of a variable, the series having coefficients periodic in time; the other was based on phase trajectories of periodic solutions. What is shown here is the possibility of generating the same tier of periodic solutions of Lyapunov systems by examining their trajectories in the configuration space. These solutions possess all the properties of normal vibrations.

Nonlinear normal vibration modes (NNMs) are a generalization of the normal (principal) vibrations of linear systems. In a normal mode a finite-dimensional system behaves like a conservative one having a single degree of freedom. In this case all position coordinates can be parametrized by any one of them as follows:

x i = p i x x x 1 i = 2 3 : : : n (2) p i x being analytical functions.

Kauderer [START_REF] Kauderer | Nichtlinear Mechanik[END_REF] became a forerunner in developing quantitative methods for analyzing NNMs. The first formulation and development of the theory of NNMs can be attributed to Rosenberg and his co-workers [START_REF] Rosenberg | On the natural modes and their stability in nonlinear two-degrees-offreedom systems[END_REF][START_REF] Rosenberg | The normal modes of nonlinear n-degrees-of-freedom systems[END_REF][START_REF] Rosenberg | Nonlinear vibrations of systems with many degrees of freedom[END_REF]. Rosenberg considered n degrees of freedom conservative oscillators and defined NNMs as 'vibrations in unison', i.e., synchronous periodic motions during which all coordinates of the system vibrate equiperiodically, reaching their maximum and minimum values at the same instant of time.

Rosenberg was the first to introduce a broad class of essentially nonlinear conservative systems allowing normal vibrations with rectilinear trajectories in a configuration space of the form:

x i = k i x i i = 2 3 : : : n : [START_REF] Rosenberg | On the natural modes and their stability in nonlinear two-degrees-offreedom systems[END_REF] For example, 'homogeneous systems' whose potential is an even homogeneous function of the variables belong to such a class. It is interesting to note that the number of modes of normal vibrations in the nonlinear case can exceed the number of degrees of freedom of the system. This remarkable property has no analogy in the linear (non-degenerate) case. In general, one expects the trajectories of normal vibrations of nonlinear systems to be curved instead of straight lines.

For some particular cases curvilinear trajectories were defined by Rosenberg and Kuo [START_REF] Rosenberg | Nonsimilar normal mode vibrations of nonlinear systems having two degrees of freedom[END_REF] and by Rand [START_REF] Manevich | Periodic solutions close to rectilinear normal vibration modes[END_REF]. In the paper by Manevich and Mikhlin [START_REF] Manevich | Periodic solutions close to rectilinear normal vibration modes[END_REF] the power series method was proposed for the construction of above mentioned trajectories.

Different results and new ideas concerning NNMs of conservative systems were performed by Vito [START_REF] Vito | An approximate method for treating the nonlinear vibration of certain two degrees of freedom systems[END_REF], Mikhlin [START_REF] Mikhlin | The joining of local expansions in the theory of nonlinear oscillations[END_REF], Vedenova et al. [START_REF] Vedenova | Normal oscillations of a string with concentrated masses on non-linearly elastic support[END_REF], Vakakis [START_REF] Vakakis | Non-similar normal oscillations in a strongly non-linear discrete system[END_REF], Vakakis and Rand [START_REF] Vakakis | Normal modes and global dynamics of a two-degree-of-freedom non-linear system. I. Low energies[END_REF], Shaw and Pierre [START_REF] Shaw | Nonlinear normal modes and invariant manifolds[END_REF][START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF], Nayfeh and Nayfeh [START_REF] Nayfeh | On nonlinear modes of continuous systems[END_REF], etc.

In [START_REF] Vito | An approximate method for treating the nonlinear vibration of certain two degrees of freedom systems[END_REF] NNMs are approximated by harmonic functions. In [START_REF] Mikhlin | The joining of local expansions in the theory of nonlinear oscillations[END_REF] Padé approximations are used for an analysis of the NNMs with large amplitudes. In [START_REF] Vedenova | Normal oscillations of a string with concentrated masses on non-linearly elastic support[END_REF] nonlinear mode localization is studied in discrete nonlinear systems with impact nonlinearities. In [START_REF] Vakakis | Non-similar normal oscillations in a strongly non-linear discrete system[END_REF] asymptotic methodology is used; for an oscillator with weak coupling stiffness in both localized and non-localized modes are detected. In [START_REF] Vakakis | Normal modes and global dynamics of a two-degree-of-freedom non-linear system. I. Low energies[END_REF] the global dynamics of strongly nonlinear systems are analyzed by means of Poincaré maps. In [START_REF] Shaw | Nonlinear normal modes and invariant manifolds[END_REF][START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF] the authors reformulated the concept of NNMs for a general class of nonlinear discrete oscillators without assuming the existence of an analytic first integral of motion. Their analysis is based on the computation of invariant manifolds of motion on which the NNM oscillations take place. In [START_REF] Nayfeh | On nonlinear modes of continuous systems[END_REF] a computationally efficient extension of the invariant manifold methodology (complex invariant manifold formulation) is proposed.

Here, a rigorous perturbation methodology for analyzing finite amplitude NNMs of broad classes discrete nonlinear systems is presented and the effect of internal resonances is considered.

Normal Vibrations in Lyapunov Systems

Let us now consider an n degrees of freedom conservative system of the form

m i x i + x i = 0 _ x i = dx i dt z = @ @z i = 1 2 : : : n (4)
= (x) being potential energy assumed to be a positively definite function and x = (x 1 x 2 : : : x n T ; power series expansion for x begins with terms having a power of at least 2. Without reducing the degree of generalization we assume that m i 1, since this can be always ensured by dilatation of coordinates.

The energy integral of the system (4) is given by

1 2 n X k =1 x 2 k x 1 x 2 x n h ( 5
)
where h is the fixed level of the total energy. Assume that within the region of the configuration space, bounded by a closed maximum equipotential surface h , the only equilibrium position is

x i 0 (i 1 2 n ).
In order to determine the trajectories of normal vibrations (2), the following relationships can be used [START_REF] Rosenberg | The normal modes of nonlinear n-degrees-of-freedom systems[END_REF][START_REF] Manevich | Periodic solutions close to rectilinear normal vibration modes[END_REF]:

2x 00 i h 1 + P n k=2 (x 02 k ) + x 0 i ( x ) = x i (i = 2 3 : : : n ; X X 1 ) (6)
where primes represent differentiation with respect to x.

These are obtained either as Euler equations for the variational principle in the Jacobi form [START_REF] Synge | Classical Dynamics[END_REF] or by elimination of time from the equations of motion (4) taking into account the energy integral [START_REF] Rosenberg | Nonlinear vibrations of systems with many degrees of freedom[END_REF]. Equations ( 6) are nonlinear and non-autonomous, and have removable singular points: they are not simpler in the linear case. Nevertheless, these relationships are suitable for the determination of nearly rectilinear trajectories of normal vibrations.

An analytical extension of the trajectories up to a maximum isoenergy surface = h is possible if the following conditions of orthogonality of a trajectory at = h are satisfied [4]:

x i x (X x 2 (X) : : : x n (X)) = x i (X x 2 (X) : : : x n (X))

(X x 2 (X) : : : x n (X)) being the end points of the trajectory on = h, where all velocities are equal to zero. If a trajectory x i (x) is defined, the law of motion with respect to time can be found using the relation:

x + x (x x 2 (x) : : : x n (x)) = 0:

This can be rewritten in the form x + V 0 (x) = 0:

The function x(t) now is the inversion of the quadrature

t + = 1 p 2 x Z X d p h V () :
Without loss of generality, the phase here is so selected that the initial velocity equals zero.

The amplitude and the energy are related by h = V (x). This equation is solvable with respect to the amplitudes X at a given energy h, provided the conditions of the closeness of all equipotential surfaces (at various values of h) holds and no equilibrium positions exist other than x i = 0. This condition also assures that the function x(t) is periodic [START_REF] Wintner | The Analytical Foundations of Celestial Mechanics[END_REF].

From the above discussion it is concluded that normal vibrations constitute a one-parameter (in the energy h) family of periodic solutions with smooth trajectories in a configuration space. A second arbitrary parameter, i.e. the phase , may also be considered, but this requires that the condition _ x(0) = 0.

Let us select terms of the potential energy that are quadratic in x j and define a generating linear system = 2 + Ñ, where the expansion of Ñ in terms of x j begins with at least a power of 3. The characteristic equation of the generating linear system has but purely imaginary roots of the form i! j (j = 1 2 : : : n ), i 2 = 1. A normal form of vibrations corresponds to a pair of roots. Let us select one of these vibration forms by fixing the pair i! (! = ! 1 ).

A basic assumption made by Lyapunov in his construction of NNMs was that there are no linearized normal modes with natural frequencies which are integrally related, i.e., which satisfy relations of the form ! i = p! k , p = 1 2 3 : : :

Transforming the quadric 2 to a sum of squares (this transformation is equivalent to the introduction of normal coordinates in a linear system described previously), one obtains:

= 1 2 ! 2 x 2 + 1 2 n X j=2 ! 2 j x 2 j + Nx x 2 : : : x n x x 1 : (8)
It is well known that such a transformation does not change the form of the kinetic energy K = 1=2 P n j=1 x 2 j . A generating of normal forms of vibrations in a linear system is now determined by equations x j = 0 (j = 2 3 : : : n ), x = X cos!t+ '. By selecting the time origin so that _ x = 0, one obtains ' = 0. Following Lyapunov, let us assume that all variables x j (j = 2 3 : : : n ) are of the order of magnitude c, where c = x0, _ x0 = 0. On having introduced c into the potential energy (8) by substituting x j ! cx j it is seen that every homogeneous constituent of contains c to a power equal to the degree of homogeneity as a factor. Thus we can write:

= c 2 2 0 @ ! 2 x 2 + n X j=2 ! 2 j x 2 j 1 A + 1 X k=3 c k N k x x 2 : : : x n (9)
where the functions N k represent the Okx k k terms of the potential energy. In the new coordinates we have x0 = X = 1. Later the symbol k will be indicated the same.

Let us now use equations ( 6) to find the normal oscillation trajectory of a nonlinear system which, on its linear limit, becomes rectilinear, x j = 0 (j = 2 3 : : : n ). The solutions of the equations will be presented as

x j = 1 X `=1 c `xj : (10)
Since we are dealing with trajectories in a configurational space, the variables x j `should be regarded as functions of x, rather than of time t.

Before proceeding to construct the solution it is worth noting that the singular points of equations ( 6) are roots of the equation: h = x x 2 : : : x n : [START_REF] Vedenova | Normal oscillations of a string with concentrated masses on non-linearly elastic support[END_REF] As the order of approximation varies so does the estimate of the energy h, since the solution x i x, i = 2 : : : n , becomes more accurate through the inclusion of higher order terms.

Suppose that the energy h 0 of the generating linear system does not coincide with the total energy h of the nonlinear system. Let us assume

h = c 2 h 0 + c 3 h i (12)
(h 0 has an order of c 2 while nonlinear terms are of greater order than c). In particular, having estimated the functions x i x, i = 2 : : : n , up to O" p , for a given amplitude of vibration x = X, the relation [START_REF] Vedenova | Normal oscillations of a string with concentrated masses on non-linearly elastic support[END_REF]

leads to: h = x p X `=1 c `x1` : : : p X `=1 c `xn`! :
This relationship enables the finding of the energy h ` ; since the energy of the linear system h 0 is known, one obtains the estimate for the energy, h 1 Substituting the series ( 10) into ( 6) and collecting the lowest order terms in c, one obtains the following equations:

2x 00 i1 h 0 1 2 ! 2 x 2 x 0 i1 ! 2 x + ! 2 i x i1 + N 3 x i x 0 : : : 0 = 0 ( 13
)
where N 3

x i x 0 : : : 0 i x 2 i = 2 3 : : : n :

At this point the boundary conditions ( 7) are used. These determine the orthogonality of the trajectory to the maximum equipotential surface. The amplitude values of x = X related to the energy h 0 by

h 0 = 1 2 ! 2 X 2 0 @ X 1 2 = s 2h 0 ! 2 1 A :
Retaining the lowest order terms in c, one derives the following set of n iOc 2 boundary conditions:

x 0 i` X! 2 X + ! 2 i x i` X + N 3 x i X 0 : : : 0 = 0 i = 2 3 : : : n : (14)
The homogeneous equations corresponding to ( 13) are hypergeometric equations (with two regular singular points) whose solutions have been thoroughly studied in the literature. A general solution of ( 13) exists in closed form. Alternatively, the analytical solution of ( 13) can be represented in terms of the Taylor series about the origin of the configuration space, which can then be analytically continued up to the maximum equipotential surface by satisfying boundary conditions [START_REF] Shaw | Nonlinear normal modes and invariant manifolds[END_REF]. To perform this computation the approximations x i1 x are expressed as

x i1 = 1 X j =0 a 1 ij x j : (15)
Substituting the series into (13), the following binomial recurrent relationships for the coefficients a 1 ik result: 4ha

1 i2 + ! 2 i a 1 i0 = 0 12ha 1 i3 ! 2 i a 1 i1 + ! 2 i a 1 i1 = 0 24ha 1 i4 2! 2 a 1 i2 2a 1 i2 ! 2 + ! 2 i a 1 i2 + i = 0 40ha 1 i5 6! 2 a 1 i3 3a 1 i3 ! 2 + ! 2 i a 1 i3 = 0
: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :

2j + 2j + 1ha 1 i j+2 jj 1! 2 a 1 ij + j! 2 a 1 ij + ! 2 a 1 ij = 0
: : : : : : : : :

(h = h 0 in this approximation).

The radius of convergence of ( 15) is found by examining the infinite recurrent set of equations ( 16),

R = lim j!1 a 1 ij a 1 i j+2 = lim j!1 2j + 2j + 1h jj i! 2 = 2h ! 2 :
Note here that, as shown in [START_REF] Rosenberg | Nonlinear vibrations of systems with many degrees of freedom[END_REF], any solution of x j x converging in a domain bounded by a surface x x 2 : : : x n = h, may be analytically continued up to this surface, provided that the orthogonality conditions [START_REF] Rand | A higher order approximation for nonlinear normal modes in two degrees of freedom systems[END_REF] hold, which in the first approximation are determined by [START_REF] Shaw | Nonlinear normal modes and invariant manifolds[END_REF].

The recursive formulas ( 16) provide non-unique solutions for the coefficients of the series [START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF]. Uniqueness of the solution is obtained by imposing the n 1 boundary conditions [START_REF] Shaw | Nonlinear normal modes and invariant manifolds[END_REF], i.e., by continuing the Taylor series expansions up to the maximum equipotential surface.

Using [START_REF] Nayfeh | On nonlinear modes of continuous systems[END_REF], one expresses the arbitrary coefficient a 1 ij in terms of coefficients a 1 i0 and a 1 i1 .

These last coefficients are computed by employing relations [START_REF] Shaw | Nonlinear normal modes and invariant manifolds[END_REF]. On substituting ( 15) into (14) (at x = X, where X 2 = 2h=! 2 ), one obtains the additional algebraic relations (at x = X, where X 2 = 2h=! 2 ):

! 2 x 2 4 1 X j=1 ja 1 ij X j 1 3 5 + ! 2 i 2 4 1 X j=1 a ij X j 1 3 5 + i X 2 = 0: (17)
Introducing at this point the quantities

K i j = j 2 2 i j + 2j + 1 i = ! i !
the solutions to the recursive relations ( 16) are expressed as:

a i 2k = K i 2k 2 K i 2k 4 : : : K i 4 2 4 K i 2 K i 0 a i i0 ! 2 2h ! k ! 2 2h ! k 1 i 12! 2 3 5 a i 2k+1 = K i 2k 1 K i 2k 3 : : : K i 3 K i 0 a i i0 ! 2 2h ! k : (18)
Expressions ( 16) relate an arbitrary coefficient a i ij to the leading coefficients a i i0 and a i i1 .

These last coefficients are determined by substituting [START_REF] Wintner | The Analytical Foundations of Celestial Mechanics[END_REF] into the boundary conditions [START_REF] Synge | Classical Dynamics[END_REF], resulting in the set of n non-homogeneous algebraic equations of the following form:

R i 0 a i i0 R i 1 R 1 i1 = R i 2 i = 2 3 : : : n (19)
where the computation of coefficients R i 0 and R i 1 requires some algebraic manipulations:

R i 0 = 2K i 0 + 4K i 0 K i 2 + 2 i 1 + K i 0 + K i 0 K i 2 + = 2 i + 2 2 i K i 0 + 4 2 i K i 0 K i 2 + = 2 i + 1 X m=1 2m 2 i K i 0 K i 2 : : : K i 2m 2 = 2K i 0 + 1 X m=1 [ 2m(2m 1) + ( 2m) 2 2 i ]K i 0 K i 2 : : : K i 2m 2 = 2K i 0 + 1 X m=1 [ 2m(2m 1) + ( 2m + 2)(2m + 1)K i 2m ]K i 0 K i 2 : : : K i 2m 2 = K i 0 K i 2 : : : K i 2m : : : R i 1 = 1 + 3K i 1 + 5K i 1 K i 3 + 2 i (1 + K i 1 + K i 1 K i 3 + ) = (1 2 i ) + ( 3 2 i )K i 1 + ( 5 2 i )K i 1 K i 3 + = (1 2 i ) + 1 X j=1 (2j 2 i )K i 1 K i 3 : : : K i 2j 1 = 6K i i + 1 X j=1 [ 2j(2j + 1) + ( 2j + 1) 2 2 i ]K i 1 K i 3 : : : K i 2j 1 = 6K i i + 1 X j=1 [ 2j(2j + 1) + ( 2j + 3)(2j + 2)K i 2j+1 ]K i 1 K i 3 : : : K i 2j 1 = K i 1 K i 3 : : : K i 2j+1 : (20)
In order to obtain unique and non-trivial solutions for the coefficients a i i0 and a i i1 , it is necessary that the coefficients of the homogeneous parts of (19) satisfy the conditions R i 0 6 = 0, and R i 1 6 = 0, i = 2 : : : n . Examining the analytical expressions [START_REF] Bellman | Introduction to Matrix Analysis[END_REF], it is concluded that in the critical case when R i 0 = R i 1 = 0, a subset of coefficients

K i j = j 2 2 i (j + 2)(j + 1)
vanishes, or equivalently, that the linearized natural frequencies of the system satisfy a resonance relation of the form ! i = j!for some positive integers j = 1 2 : : : These were precisely the cases which were eliminated from consideration in Lyapunov's analysis [START_REF] Lyapunov | The General Problem of the Stability of Motion[END_REF].

Let us now consider the higher order approximations. Suppose that in the series (9) singlevalued solutions x ik (x) are computed, where k l 1, which are analytical functions over a closed domain bounded by the surface = h. Considering O(c `) terms in (6) one obtains "

2x 00 i` h 0 1 2 ! 2 x 2 + 2x 00 i ` 1 h ` 1 1 x 0 i`! 2 x + x i`! 2 i + ` 1 X k=1 x 00 ik R 1 ` k i + ` 1 X k=1 x 0 ik R 2 ` k i [R 3 ` 1 i ] = 0
x 00 i1 ! 2 x + x i1 ! 2 i + ` 1 X k=1 x 0 ik R 2 ` k i R 3 ` 1 i # x=X j = 0
where x = X j (j = 1 2) are vibration amplitudes (one of these values, X 1 = 1).

The equations ( 21) may be rearranged as follows:

2x 00 i` h 0 1 2 ! 2 x 2 x 0 i`! 2 x + x i`! 2 i + F 1 i x = 0 (22)
where the terms F ` i x consist of already-computed functions of x. Expressing the solution of ( 22) in the series form:

x

i`= 1 X j=0 a ` ij x j
one obtains, similarly to the first approximation case, a non-homogeneous recurrent set of linear equations governing the coefficients a ` ij . The indefinite coefficients a ` i0 , a ` i1 are obtained from the boundary conditions for this approximation (conditions of orthogonality). Next the Oc ` approximation for the total energy, h ` is determined from [START_REF] Vedenova | Normal oscillations of a string with concentrated masses on non-linearly elastic support[END_REF], while the equation

h ` = c 2 h 0 + c 3 h `
1 is used to find the value of h ` 1 which appears in the calculations of the next order of approximations.

Needless to say, an alternative series of calculations would be also acceptable, namely, determining the maximum amplitude X given a fixed level of total energy h.

The conclusions concerning the convergence and unambiguous definitions of the coefficients arrived at for the series in the first approximation hold for the series in the `-th approximation as well.

Let us consider the convergence of the series [START_REF] Mikhlin | The joining of local expansions in the theory of nonlinear oscillations[END_REF]. As shown, a series of the form (10) represents a single-valued (provided K i j 6 = 0 holds) formal solution of a boundary problem ( 6), [START_REF] Rand | A higher order approximation for nonlinear normal modes in two degrees of freedom systems[END_REF] with the coefficients analytical in x. In this case x j (0) = a j1 , x 0 j (0) = a j2 , where a j1 , a j2 can be made sufficiently small by choosing a small value of the parameter c.

Over a domain < h all functions involved in (5) are analytical in x. Therefore, it follows from Poincaré's theorem on the small parameter series expansion [START_REF] Poincaré | Les méthodes nouvelles de la mécanique celeste. 3 Volumes[END_REF] that there is a value of c 0 > 0 such that, for all jcj < c 0 , the series of the form (10) converges in the domain, represents a unique solution of (6) analytical in c and x, and satisfies the conditions

x j (0) = a j1 , x 0 j (0) = a j2 ; moreover, as c ! 0 this solution becomes the trivial generating solution x j0 = 0 (j = 2 3 : : : n ). Since series of the form (10) also satisfy the conditions [START_REF] Rand | A higher order approximation for nonlinear normal modes in two degrees of freedom systems[END_REF], the solution can be analytically continued up to the domain boundary = h.

On having obtained the smooth trajectory, the problem of finding a periodic solution reduces to the integration of a conservative system with one degree of freedom.

Concluding this treatment of normal vibrations in Lyapunov systems, we note the following.

Firstly, the trajectories x j (x) can be derived not only in terms of power series in x, but also by the method of successive approximations. This method is described below for certain cases of nonlinearity of the generating system. Secondly, the requirement of the existence of an energy integral ( 5) is not essential. In autonomous systems of the form

x i = f i (x 1 x 2 : : : x n )
the equations for obtaining the trajectories x i (x), (x x 1 ) take the form

x 00 i _x 2 + x 0 i f 1 (x x 2 : : : x n ) = f i (x x 2 : : : x n ): (23)
If an analytical first integral H(x _x x 2 _x 2 : : : x n _x n ) = 0 exists, all that is needed is to obtain _x as a single-valued analytical function of x, x 2 , x 0 2 , : : : , x n , x 0 n , and to substitute the expression obtained into [START_REF] Zhupiev | Conditions for finiteness of the number of instability zones in the problem of normal vibrations of non-linear systems[END_REF]. Equations [START_REF] Zhupiev | Conditions for finiteness of the number of instability zones in the problem of normal vibrations of non-linear systems[END_REF] can then be analyzed similarly as above.

Finally, it should be noted that an approach involving the examination of trajectories in a configurational space is by no means more complicated than approaches employing series with coefficients periodic in time or examining phase trajectories.

Normal Oscillations in Conservative Systems Close to Systems Admitting Rectilinear Modes

Consider now a system of a more general class, namely one with a nonlinear generating system. The smallness of disturbances will be evaluated using a small parameter . The equations of motion are expressed in the form

q i F 0 q i q 1 q 2 q n F 1 q i q 1 q 2 q n 0 i 1 2 n where F 0 q i F 1 q
i is the analytical potential energy, and the unperturbed systems ( 0) possesses rectilinear modes of vibrations

q i k i q (q q 1 , i 2 3 n
). Consider one of these solutions as a generating one. Rotating the coordinate axes, so that the new x -axis is directed along the rectilinear trajectory, and the remaining coordinate axes are orthogonal to it. A generating solution in the new coordinates is represented as x i = 0 (i = 2 3 : : : n ), x 1 x = x(t).

In particular, for a system with two degrees of freedom q i + F qi (q 1 q 2 ) = 0 (i = 1 2) allowing the solution q 2 = k q 1 , the formulae for the axes rotations are of the form

x 1 = q 1 cos # + q 2 sin # x 2 = q 1 sin # + q 2 cos # (# = arctan k):

(

) 24
The equations of motion in the new coordinates take the following form:

x i + x i (x 1 x 2 ) = 0 (i = 1 2)
where

x 1 = F q 1 cos # + F q 2 sin # x 2 = F q 1 sin # + F q 2 cos #:
Let us assume that in the general case (with n degrees of freedom) after the rotation of the axes, the system

x i + 0 x i (x 1 x 2 : : : x n ) + " 1 x i (x 1 x 2 : : : x n ) = 0 ( 25
)
also admits a generating solution at x i = 0 at " = 0; but this means that 0 x i (x 1 0 : : : 0) 0 (i = 2 3 : : : n ):

(26)
It is initially assumed that the unperturbed system is homogeneous, i.e., 0 is an even homogeneous function of the power of r + 1 in all the variables (r may take the following values r = 1, 3, 5, : : : ).

Note that, in similarity to the Lyapunov systems examined in the previous section, the small parameter could be chosen, to scale the amplitude of vibration c = x(0). Substituting for x i ! cx i , one then selects a generating homogeneous system containing the smallest powers of the positional variables. This generating system may be nonlinear as well if r > 1.

In order to determine the trajectories of normal oscillations in the configurational space, another transition is required, this time from [START_REF] Nayfeh | Nonlinear Oscillations[END_REF] to (6) used in combination with the boundary conditions [START_REF] Rand | A higher order approximation for nonlinear normal modes in two degrees of freedom systems[END_REF].

A solution is sought in the form of a small parameter series:

x i = 1 X k=1 " k x ik (x) (x x 1 ): (27)
Although in this case the generating system is essentially nonlinear, all computations are similar to those performed for the linearizable case. Let us assume that the total energy h of the entire system and the energy h 0 of the generating system are related as follows:

h = h 0 + "h 1 :
Substituting (27) into (6), one isolates the `-th approximation equations in ":

2x 00 i`h h 0 0 (x 0 : : : 0) i + 2x 00 i ` 1 h h ` 1 1 + x 0 i` 0 x (x 0 : : : 0) i + n X k=2 0 x i x k (x 0 : : : 0)x k`+ ` 1 X j=1 x 00 ij R 1 ` j i + ` 1 X k=1 x 0 ik R 2 ` k i + R 3 ` 1 i = 0 (28)
where R s mi = X @ P s ri x 0 : : : 0 @x 2m 2 @x 3m

: x n " 1 + n X k=2 x 0 k 2 ## r P 3 ri = " x i x x 2 : : : x n " 1 + n X k=2 x 0 k 2 ## r :
The boundary conditions (conditions of orthogonality) corresponding to this approximation assume the following form: "

x 0 i` 0 x x 0 : : : 0 n X k=2

0 x i x k x 0 : : : 0x k+ ` 1 X k=1 x 0 ik R 2 ` k i R 3 ` 1 i # x=X j = 0 j = 1 2 (29)
with the amplitudes X j (j = 1 2) and the energy h being computed by h = x x 2 : : : x n where x k (k = 2 3 : : : n ) should be expressed by the series (27) worked out to an order of

" ` 1 .
Since the unperturbed system is homogeneous, the matrix B = 0 x i x k x 0 : : : 0 may be written as B = b ik x r 1 . Note that, owing to the conservative nature of the system, b ik = b ki , and a symmetric matrix b ik is reduced to a diagonal form by an invertible linear transformation of coordinates [START_REF] Bellman | Introduction to Matrix Analysis[END_REF]. Therefore, without loss of generality, one can assume that in (28) and (29) the function 0 x i x k x 0 : : : 0 = 0 for i 6 = k. Hence, only the terms 0 x i x i x 0 : : : 0x ià re retained in the expression P n k=2 0 x i x k x 0 : : : 0x k`, and the set (28), together with the boundary conditions (29) is 'split' in the variables x i`.

Having substituted the series x i`= P 1 j=0 a ` ij x 0 into (28), one finds that the coefficients a ` ij are interrelated by an infinite set of recurrent relationships 2h 0 r + j + 2r + j + 1a ` i r+j+2 jj + 12 0 1 0 : : : 0a ` i j+1 j + 1 0 x 1 0 : : : 0a ` i j+1 + 0 x i x i 1 0 : : :

0a ` i j+1 = ` ij (30)
where ` ij denotes the terms that depend on the solutions of preceding approximations. The recurrent relationships can be used to express all coefficients a ` ij in terms of the leading coefficients a ` i0 and a ` i1 (i = 2 3 : : : n ). The convergence of the series is ascertained precisely as it is done with Lyapunov systems over a domain x x 2 : : : x n < h. An analytical continuation of the solution up to the 'boundary' x x 2 : : : x n = h may be effected, provided that the boundary conditions are satisfied. Substituting the series x i` x into these boundary conditions (29), and in view of the recurrent relationships (30), one obtains the equations governing a ` i0 , a ` i1 .

n X i=2 R 0i a ` i0 R ` 1i a ` i1 = R 1 2i i = 2 3 : : : n :
As in the case of Lyapunov systems, the determinants of coefficients may be represented as products of an infinite number of factors K p = jq ij j q ij = j i h pp 12 0 1 0 : : : 0 + p 0 x 1 0 : : : 0 0 x i x i 1 0 : : :

0 i (31)
where j i are the Kronecker's delta, p = 0, 1, 2, 3, : : :

When the generating system is linear, the solvability conditions (31) can be shown to degenerate to the conditions of absence of internal resonances. Hence, the conditions K p 6 = 0 can be viewed as generalizations of the conditions of the absence of internal resonances derived in the linearizable case, and ensure that the analytical, asymptotic solutions x j = x j x are unique and single-valued.

The series x ` j = P k=1 " k xx j k (j = 2 3 : : : n ) may now be substituted into the equation h = x x 2 : : : x n j x=X and a refined value of the total energy approximation h ` may be obtained at a given amplitude x = X ( _x = 0). Hence, h `

1 may be deduced from h = h 0 + "h 1 to be employed in the next approximation in ".

The considerations concerning the convergence of the series (27) given in the preceding section still hold.

Note that in the absence of the assumption of 'splitting' of the first approximation equations, the conditions (31) for the generating of homogeneous systems become different. In the general case they take the following form: K p = jq ij j 6 = 0 q ij = j i h pp 12 0 1 0 : : : 0 + p 0 x 1 0 : : : 0 0 x i x j 1 0 : : : 0 i p = 0 1 2 3 : : : :

It is worth noting that at r = 1 the conditions (32) become ! i 6 = p! which for Lyapunov systems meant that in this case inner resonances in the generating linear system are not considered.

Proceeding to non-homogeneous generating systems with the potential energy 0 , we shall confine the discussion with the first approximation in ":

2x 00 i`h h 0 x 0 : : :

0 i + x 00 i`h 0 x x 0 : : : x i + n X k=2 0 x i x k
x 0 : : : 0 x

The homogeneous part of the set (33) contains variational equations for the generating of normal vibrations. Let us assume that the variational equation set may be 'split' into n independent equations by an invertible linear transformation with constant coefficients. The possibility of such 'splitting' is discussed in more detail in [START_REF] Hsu | On a restricted class of coupled Hill's equations and some application[END_REF][START_REF] Zhupiev | Stability and branching of normal oscillation modes of nonlinear systems[END_REF].

Here it will be shown that, for the case of two degrees of freedom, it is always possible to effect the 'splitting' by rotating the axes [START_REF] Mikhlin | Matching of local expansions in the theory of nonlinear vibrations[END_REF]. Indeed, introduce the transformation [START_REF] Mikhlin | Matching of local expansions in the theory of nonlinear vibrations[END_REF], and expand the potential energy 0 x 1 x 2 into an x 2 (x 1 x) power series:

0 x x 2 0 x 0 0 x 2 x 0 x 2 1 2! 0 x 2 x 2 x 0 x 2 2 1 3! 0 x 2 x 2 x 2 x 0 x 3 2
The condition of existence of a normal vibration form x 2 0 is given by 0 x 2 x 0 0, and the equations of motion take the form

x 0 x x 0 d dx 0 x 2 x 2 x 0 x 2 2 2 O x 3 2 0 x 0 x 2 x 2 x 0 x 2 d dx 0 x 2 x 2 x 2 x 0 x 2 2 2 O x 3 2 0:
Here, the variational equations in u and v for the normal vibrations x x t , x 2 0 are as follows:

u d dx 0 x 2 x 2 x 0 u 0 v 0 x 2 x 2 x 0 v 0: (34)
The set of equations is thus 'split'.

Introducing a new independent variable x instead of t, one obtains, in the place of (34), 2v 00

h h 0 x x 0 i v 0 h 0 x x 0 i 0 x 2 x 2 x 0 v 0: (35)
This equation corresponds to (33) at n 2 when the latter retains only a homogeneous constituent of the potential energy.

Since a rectilinear normal mode of vibration has only two cusps, the kinetic energy k h 0 x 0 vanishes twice over one period. Therefore, equation (35) has two regular finite singular points on the real axis.

There exist classes of potentials for which equation (35) has been studied so thoroughly that it is possible to develop a system of solutions either in power series (with due regard to singularities), in trigonometric series, or in the form of series in certain special functions. If 0 x 1 x 2 is a homogeneous even function, (35) may even be reduced to a hypergeometric equation by the substitution x p z, where p is the degree of homogeneity [START_REF] Zhupiev | Stability and branching of normal oscillation modes of nonlinear systems[END_REF]. If the potential energy 0 x x 2 contains terms of second and fourth powers of x and x 2 , (35) is the Lame equation [START_REF] Zhupiev | Conditions for finiteness of the number of instability zones in the problem of normal vibrations of non-linear systems[END_REF]. Denote by fw 1 x w 2 x g a fundamental set of solutions for (35), and it becomes possible to construct a general solution of the first approximation equation in " (33): should be determined from the regularity conditions for the solution in any point of the interval 0 x 0 h. Suppose that the powers of the singular points X j (j = 1 2) (zero kinetic energies) be equal to 0 1 and 0 2 , respectively. In particular, for the Lame equation 1 = 2 = 1=2. Then, in the vicinity of these points the general solution is decomposed as follows:

x 21 = A 1 C 1 C 2 f 1 x + A 2 C 1 C 2 f 2 xx X 1 1 x 21 = B 1 C 1 C 2 g 1 x + B 2 C 1 C 2 g 2 xx X 2 2 :
Evidently, the regularity conditions for this solution are of the form

A 2 C 1 C 2 = B 2 C 1 C 2 = 0:
A solution analytical in x is similarly derived in the next approximation in " obtained either through the construction of vibration normal modes as " series or by the iteration method.

A generalization to the case of n degrees of freedom (n > 2) is evident, provided the set of variational equations can be 'split'. Some applications of the asymptotic methodology for approximating the NNMs are presented in [START_REF] Vakakis | Non-similar normal oscillations in a strongly non-linear discrete system[END_REF][START_REF] Mikhlin | Matching of local expansions in the theory of nonlinear vibrations[END_REF]. In particular, in [START_REF] Vakakis | Non-similar normal oscillations in a strongly non-linear discrete system[END_REF] the asymptotic methodology is implemented to analyze the free vibrations of two unit masses connected by means of three strongly nonlinear stiffnesses with cubic nonlinearity. The zeroth order and O" approximations to NNMs are calculated.

The systems considered above can be obtained in calculations of nonlinear vibrations of rods, strings, plates, shells and other elastic systems (using the Bubnov-Galerkin technique) (see, for instance, [START_REF] Kauderer | Nichtlinear Mechanik[END_REF][START_REF] Nayfeh | Nonlinear Oscillations[END_REF]). Some examples of the NNMs analysis in nonlinear elastic systems using Bubnov-Galerkin discretization are presented in [START_REF] Manevich | The Method of Normal Oscillation for Essentially Nonlinear Systems[END_REF].

Conclusions

Normal modes of vibrations with curvilinear trajectories in the configuration space are considered. The corresponding boundary problem is formulated. It is shown that such normal modes may be constructed for the Lyapunov systems and for a more general class of conservative systems, neighboring systems possessing NNMs with rectilinear trajectories. The convergence of the power series obtained was considered. Under some conditions (that generalize the conditions of internal resonances in the quasilinear case) single-valued solutions do not exist.

  ) pi (s = 1 2 3), and the second summation sign in the above expression is carried out over all positive integer solutions of the equation are complemented by the following set of Oc 1 boundary orthogonality conditions:

3 :

 3 above expression is carried over all positive integer solutions of the equation where

2 x 0 w 2 dx w 2 w 1 ,

 221 w 2 are either hypergeometric or Lame functions. The constants C 1 1 , C 1 2

to be used in the next approximation.


J Nuño

M Gilles

M Guasoni

C Finot

J Fatome

All-optical sampling and magnification based on XPM-induced focusing

We theoretically and experimentally investigate the design of an all-optical noiseless magnification and sampling function free from any active gain medium and associated high-power continuous wave pump source. The proposed technique is based on the co-propagation of an arbitrary shaped signal together with an orthogonally polarized intense fast sinusoidal beating within a normally dispersive optical fiber. Basically, the strong nonlinear phase shift induced by the sinusoidal pump beam on the orthogonal weak signal through cross-phase modulation turns the defocusing regime into localized temporal focusing effects.

This periodic focusing is then responsible for the generation of a high-repetition-rate temporal comb upon the incident signal whose amplitude is directly proportional to its initial shape. This internal redistribution of energy leads to a simultaneous sampling and magnification of the signal intensity profile. This process allows us to experimentally demonstrate a 40-GHz sampling operation as well as an 8-dB magnification of an arbitrary shaped nanosecond signal around 1550 nm in a 5-km long normally dispersive fiber. The experimental observations are in quantitative agreement with numerical and theoretical analysis.

In modern photonic systems, the sampling process has widespread applications in the fields of optical communications, metrology, clocking, sensing, spectral comb or arbitrary waveform generation.

In this context, nonlinear effects have been demonstrated as potential key technologies to develop alloptical sampling devices [START_REF] Duguay | [END_REF][2][3][4][5][6][7][8][9][10][11][12][13][14] . In most of these techniques, an ultrashort pulse train acts as an optical gate and the basic physical phenomena under use include four-wave mixing (FWM) 10 , cross-phase modulation (XPM) 11 , nonlinear polarization rotation 13 and Raman soliton self-frequency shifting 14 . On the other hand, signal amplification is also a critical function in many area of physics. Basically, optical amplification process refers to multiply an incident signal by increasing its global energy through an active gain medium pumped by an external pump beam. The large range of current and mature amplification techniques at telecommunication wavelengths involves Erbium doped fiber amplifiers, Raman-based amplifiers, semiconductor amplifiers or parametric processes. However, it is noteworthy that most of these techniques lead to an inherent degradation of the signal-to-noise ratio induced by the detrimental spontaneous emission of photons. More recently, Azaña and coworkers have suggested and experimentally validated a different scenario in which a linear redistribution of energy into a periodic pulse train through the selfimaging Talbot effect can lead to a noiseless amplification process 15 . In this manuscript, we propose an alternative nonlinear approach enabling to simultaneously sample and magnify an arbitrary shaped incident signal due to a focusing effect occurring in a normally dispersive medium through the cross-phase modulation induced by an orthogonally polarized high-repetition-rate sinusoidal pump beam.

Fundamentally, the strong sinusoidal beating is responsible for a localized and periodic focusing of the signal energy into high repetition pulses generated at the pump modulation frequency 16 . This internal redistribution of energy within the signal wave leads to the generation of a temporal intensity comb whose amplitude is then directly proportional to the signal profile, thus inducing a simultaneous sampling and magnifying process. Based on this technique, we have experimentally demonstrated a 40-GHz sampling operation as well as an 8-dB magnification factor of an arbitrary shaped nanosecond signal around 1550 nm in a 5-km long normally dispersive fiber. We have also derived a theoretical estimation of the total magnifying factor of the system in good agreement with our experimental observations and numerical simulations based on a simplified Manakov modelling.

Basically, it is well-known that the anomalous dispersive regime of propagation in optical fibers can lead to a focusing effect owing to the interplay between chromatic dispersion and Kerr nonlinearity. This process has been exploited in numerous studies in order for instance to compress optical pulses or design high-repetition pulse sources [17][18][19] . More recently, we have reported the observation of a focusing effect occurring in the normal dispersion regime through a XPM interaction between a weak periodic signal and its orthogonally polarized interleaved high-power replica 16 . In this previous configuration, a 3.3-ps 40-GHz pulse train was then generated 16 . In this new contribution, we exploit a similar physical mechanism 20 to induce a temporal focusing of an arbitrary shaped signal so as to sample and magnify its intensity profile. As a rule of thumb, the signal to be processed is co-propagating in a normally dispersive optical fiber together with a sinusoidal pump beam. In the present configuration, the pump wave consists in a much more intense orthogonally polarized high-repetition-rate carrier-suppressed beat-signal. Due to the nonlinear defocusing dynamics of the fiber, the pump is reshaped into parabolic then broad and sharp square pulses. Subsequently, optical shocks associated to wave-breaking phenomenon appear in its intensity profile inducing deep singularities at its null points 21 . As these singularities propagate, they are characterized by steeper and steeper edges, and then progressively collapse. In parallel, the XPM nonlinear phase shift induced by these nonlinear dark structures combined with the group velocity dispersion of the fiber then turns out the normal dispersion regime into a periodic focusing dynamics for the co-propagating signal. As a consequence, a fast and periodic pulse train exactly proportional to the incident signal is then generated, thus inducing a discrete sampling and an overall magnification effect of its intensity profile.

Figures 1 (a-c) show a numerical illustration of this process. For these simulations, we consider two orthogonally polarized optical waves propagating in a normally dispersive randomly birefringent fiber.

The evolution of the complex slowly varying amplitudes of the signal u and pump v are described by the following set of two coupled nonlinear Schrödinger equations, corresponding to a simplified Manakov model in which nonlinear terms involved with the weaker signal have been neglected:

0 2 9 8 2 2 2 2 2     u i u v t u z u i        , 0 2 9 8 2 2 2 2 2     v i v v t v z v i        .
Here, z and t denote the propagation distance and time coordinates, respectively, whereas α indicate propagation losses. The factor 8/9 takes into account for random fluctuations of the intrinsic birefringence along the fiber length. In addition, we assume that third-order dispersion, losses or Raman effects play a (1) minor role in the nonlinear dynamics and can therefore be neglected. We also assume that the wavelengths of the pump and signal under investigation are relatively close so that the walk-off between the two waves can be neglected. An arbitrary shaped weak signal is first injected into the fiber with a moderate average power (red solid line in Fig. 1b). The fiber parameters correspond to the fiber segment used in our experiment and is characterized by a chromatic dispersion β2 = 3.2 ps 2 /km and a nonlinear Kerr coefficient γ = 1.7 W -1 km -1 . A high power sinusoidal pump wave (26.5 dBm average power) is then simultaneously injected with an orthogonal polarization state (red line in Fig. 1a). Note that to ensure an efficient signal processing, the temporal period of the pump wave has to be much shorter than the signal duration. Owing to the strong nonlinear defocusing regime, the carrier-suppressed pump wave is drastically reshaped into the fiber and sharp singularities gradually appear at its null points (see blue line in Fig. 1a after 5 km of propagation). Simultaneously, we clearly observe in Fig. 1b that the focusing effect imposed by the copropagating pump wave progressively leads to a redistribution of the signal energy. Indeed, at the output of the 5-km long fiber, an intensity comb at the pump repetition-rate is then generated upon the signal wave whose amplitude is directly proportional to its initial envelope, thus leading to both a discrete sampling and magnification effect (blue line in Fig. 1b). Finally, normalized profiles shown in Fig. 1c confirm that the profile of the magnified wave is a perfectly scaled version of the signal under test. In order to predict the magnification factor experienced upon the sampling process, it proves useful to develop a simplified model based on the analysis of the pump-signal interaction carried out in the optical spectral domain. The cosines-like input v(0,t)≡v(t)in pump reads as v(t)in=(P0) 1/2 cos(ωpt), to which correspond two spectral components centered at ±ωp. The weak input signal is assumed to be a slowlyvarying waveform u(0,t)≡uin(t), instead, whose corresponding optical spectrum is centered around ω=0 with a spectral extension well below ωp. The FWM interaction among the pump and signal leads to the generation of new spectral components centered around odd (even) multiples of ωp in the pump (signal) spectrum. On the other hand, the different power-level of pump and signal makes their longitudinal evolution deeply different. As long as the fiber length does not exceed few nonlinear lengths (defined as

1/γP0
), the level of the new spectral components generated by FWM in the pump spectrum remain 10 dB below the original components centered at ± ωp 16 . We can thus neglect as a first rough approximation the impact of these additional spectral components and therefore consider that the sinusoidal pump intensity profile is preserved upon propagation except for propagation losses, that is |v(z,t)| 2 =P0 cos(ωpt) 2 exp(-αz).

Note that this approximation remains still valid as long as the shock point undergone by the pump wave only occurs in the last stage of propagation, which is the case for most of the power values considered in the present experiment. In contrast, the generation of new spectral components cannot be neglected in the case of the weak signal, instead, as their amplitude rapidly grow up to the peak-level of the input signal spectrum U(ω)in. Replica of U(ω)in arise in the signal spectrum, centered around even-multiples of ωp, to which the following ansatz u(z,t) corresponds:

𝑢(𝑧, 𝑡) = 𝑢 𝑖𝑛 (𝑡) [∑ 𝑆 2𝑛 (𝑧) 𝑒𝑥𝑝(𝑖 2𝑛 𝜔 𝑝 𝑡)],
Insertion of the aforementioned ansatz in Eq. denotes the average pump power along the fiber length. In doing this,

we have exploited the slowly-varying hypothesis for uin(t), which allows neglecting all terms in ∂uin(z)/∂t and ∂²uin(z)/∂²t.

Solution of the LDE above requires truncation over a finite number of terms, that is -N ≤ n ≤ N, and the assignment of the input values S2n(0). In practice, N=4 (9 components) is sufficient to get a good estimation of S2n(z). Moreover the input spectrum U(ω)in is centered around ω=0, therefore we set S0(0)=1

and S2n(0)=0 for n≠0. Note that in the case of a symmetric spectrum, as in the present study, the modeling can be reduced to half of the truncated terms; that is to say only the n  0 spectral components can be considered. It is interesting to notice that according to Eq. 2 the input beam uin(t) is modulated by a function 𝑔(𝑧, 𝑡) = ∑ 𝑆 2𝑛 (𝑧)𝑒𝑥𝑝(𝑖2𝑛𝜔 𝑝 𝑡). Solution of the LDE show that g(z,t) evolves from a constantin-time function g(0,t)=1 towards a periodic train of pulses that are progressively compressed along the fiber length. Consequently, the signal u(z,t) experiences a simultaneous sampling and magnification, for which samples are more and more temporally compressed and intense upon propagation. In Fig. 1b, the maxima of the output intensity profile calculated by means of this analytical model have been reported

with black crosses and compared to numerical simulations. We can notice a good agreement between both results, thus validating our theoretical approach. Furthermore, in the inset of Fig. 1b, a close view of both intensity profiles has been depicted and also reveals an excellent prediction of the generated pulse shape.

Besides, as we will show in next section, this simplified analytic model well captures the features observed in experiments, and may allow optimizing and engineering the sampling and magnification processes here discussed. The experimental setup depicted in Fig. 2 was implemented in order to provide a proof-of-principle of this technique. For generating both the pump and signal waves, the scheme is initially divided into two arms.

In the upper side (pump generation), an external cavity laser emitting at 1550 nm is first intensity modulated thanks to an intensity modulator (IM). A 20-GHz frequency clock is used to drive the IM at its null-transmission point so as to generate a pure carrier-suppressed 40-GHz sinusoidal wave. Subsequently, the resulting beating is phase modulated to 100 MHz thanks to a phase modulator (PM) so as to significantly increase the Brillouin threshold above the power levels involved in our experiment. On the other side (signal generation), another external cavity laser emitting at the same wavelength is intensity modulated by means of a second IM, which is driven by an arbitrary waveform generator (AWG).

Therefore, the profile of this signal consists of a concatenation of two unbalanced Gaussian pulses with a duration of a few nanoseconds. Both the beat-signal (pump) and the arbitrary signal are independently amplified by two Erbium doped fiber amplifier (EDFA) and polarization multiplexed thanks to a polarization-beam splitter (PBS). The resulting signals are injected into a 5-km long optical fiber characterized by a chromatic dispersion D = -2.5 ps/nm/km at 1550 nm, an attenuation of 0.2 dB/km and a nonlinear Kerr coefficient γ = 1.7 W -1 km -1 . Finally, at the output of the system, the resulting pump and signal waves are polarization demultiplexed thanks to a second PBS and characterized in the time domain by means of two different oscilloscopes (an electrical sampling oscilloscope, ESO with a detection bandwidth of 70 GHz, and an optical sampling oscilloscope OSO from Alnair Labs having a temporal resolution of 1 ps). The output signal is also characterized in the spectral domain thanks to an optical spectrum analyzer (OSA).

The power of the incident signal is kept constant to 10 dBm while the pump power varied from 16 to 28 dBm. Figure 3a displays the normalized envelope of the arbitrary shaped nanosecond signal monitored at the output of the fiber when the sinusoidal pump wave is off (in red) and after the sampling process for an injected pump power of 26.5 dBm (in blue). We clearly observe a high-repetition-rate pulse train distributed along the signal wave and generated at the modulation frequency of the pump whose amplitude perfectly matches the initial signal envelope. We can thus notice that the sampling process acts in full strength whilst the shape of the signal is well preserved, in good agreement with our numerical predictions of Fig. 1c. It is noteworthy that for the results displayed in Fig. 3a, both signal and pump waves were synthetized thanks to the same initial RF clock and consequently appear temporally synchronized, which allows us to efficiently monitor the generation of the pulse train along the signal profile. Simultaneously, the sinusoidal pump wave has been broadened due to the combined effects of SPM and chromatic dispersion and undergoes an optical shock as illustrated in Fig. 3b for an injected power of 26.5 dBm, in full agreement with the trend reported in Fig. 1a. As it was previously mentioned, apart from the sampling process, the intensity profile of the signal wave is simultaneously magnified. The magnifying factor of the signal shape recorded at the output of the fiber as a function of the injected pump power has been reported in Fig. 4. Here, pump and signal waves are not temporally synchronized anymore. Therefore, the signal undergone an asynchronous sampling and only the envelope of the resulting magnified signal is then plotted. Fig. 4a compares the intensity profiles of a 1.5-ns arbitrary signal recorded at the output of the fiber with and without the presence of the 40-GHz sinusoidal pump wave, for different injected powers. We can notice that the output signal can reach a magnification factor close to 7 (8.5 dB) when the pump is activated. In particular, it is noteworthy that the second Gaussian pulse superimposed on the trailing edge of the arbitrary signal was hardly observable when the pump is off, whereas it can be easily detected when the signal is intensity magnified. In Figure 4b, we have reported the experimental magnifying factor (blue circles) of the output nanosecond signal as a function of the injected pump power. First of all, we can notice that the factor of magnification is not perfectly proportional to the pump power and that a maximum gain of 8.5 dB can be achieved owing to the present configuration. Experimentally, the pump power is here limited to 28 dBm because of the Brillouin backscattering effect. These experimental results have been compared both to numerical simulations (red solid line) and to our theoretical predictions (black crosses) assuming the simplified Finally, the output signal is analyzed in the spectral domain. The phenomenon can be characterized as the generation of new frequency components at the pump frequency (40 GHz) related to the time-dependence of the nonlinear phase shift induced by XPM as illustrated in Fig. 4c. Moreover, we can observe once again a good agreement between the output recording spectrum and our theoretical predictions. Note that the disagreement among theoretical and numerical/experimental spectra grows larger as the pump power increases. This discrepancy at large pump power is related to the generation of new spectral components in the pump spectrum through FWM which, combined to chromatic dispersion, leads to severe distortion of its intensity profile; effect neglected in our present theoretical approach. Nevertheless, these components mainly affects the phase profile of the generated pulse signal, while not the magnitude. For this reason a good agreement between theoretical and numerical/experimental intensity profiles is still achieved for large pump powers (see Fig. 1b). In summary, we have proposed and experimentally demonstrated a novel technique enabling to simultaneously sample and magnify an arbitrary shaped optical signal. The present system is based on a periodic focusing effect upon an incident signal which takes place in a normally dispersive medium by means of a cross-phase modulation interaction induced by an orthogonally polarized high-repetition-rate sinusoidal pump wave. This localized periodic focusing induces a redistribution of energy along the signal profile which leads to the generation of a temporal comb at the repetition frequency of the pump whose amplitude is directly proportional to the incident signal, thus enabling its sampling and magnification.

Using this technique, we have successfully implemented a 40-GHz sampling operation as well as an 8-dB magnification of an arbitrary shaped nanosecond signal around 1550 nm within a 5-km long normally dispersive optical fiber. Numerical and theoretical results are in good agreement with our experimental recordings. Alternatives based on wavelength multiplexed waves instead of polarization multiplexed waves could be developed, but the walk-off between the pump and signal may ultimately in such configuration degrade the efficiency of the magnification due to the focusing effect.

FIG. 1 .

 1 FIG. 1. Numerical simulations (a) Input (red) and output (blue) intensity profiles of the sinusoidal pump wave normalized to unity for an input average power of 26.5 dBm and 5 km of propagation. (b) Output intensity profiles of the orthogonally polarized signal wave when the pump is deactivated (in red) or activated (blue). The black crosses indicate the maximum values of the theoretical prediction. In the inset, a close view of the numerical and the theoretical estimations. (c) Same as (b) but normalized to unity.

FIG. 2

 2 FIG. 2 Experimental setup. IM: intensity modulator, PM: phase modulator, AWG: arbitrary waveform generator, EDFA: Erbium doped fiber amplifier, PC: polarization controller, PBS: polarization beam splitter, ESO: electrical sampling oscilloscope, OSO: optical sampling oscilloscope and OSA: optical spectrum analyzer.

FIG. 3 .

 3 FIG. 3. (a) Normalized output signal in the temporal domain. For the red line, the pump is off and, for the blue curve, the signal is sampled (pump power equal to 26.5 dBm). This figure is obtained with the ESO. (b) Output intensity profile of the pump beam (power equal to 26.5 dBm). This snapshot is obtained with the OSO.

  FIG. 4. (a) Output signal in the temporal domain. For the red line, the pump is off and the signal is magnified for different injected powers. This figure is obtained with the OSO and only the envelope is drawn. (b) Magnifying factor of the output signal as a function of the injected pump power. The experimental values (blue circles) are compared with the numerical simulations (red line) and theoretical predictions (black crosses). (c) Experimental spectra of the output signal (in blue) compared to theoretical (black crosses) and numerical estimations (red dashed line) for a pump power of 26.5 dBm.

Acknowledgements This research is funded by the European Research Council under Grant Agreement