John T Melcher

Alan R Champneys
email: a.r.champneys@bristol.ac.uk

David J Wagg

The impacting cantilever: modal non-convergence and the importance of stiffness matching

The problem of an Euler-Bernoulli cantilever beam whose free end impacts with a point constraint is revisited from the point of view of modal analysis. It is shown that there is non-uniqueness of consistent impact laws for a given modal truncation. Moreover, taking an N -mode compliant, bilinear formulation and passing to the rigid limit leads to a sequence of impact models that does not converge as N → ∞. The dynamics of such truncated models are studied numerically and found to give rise to quite different dynamics depending on the number of degrees of freedom taken. The simulations are compared with results from simple experiments which show a propensity for multiple-tap dynamics, in which higher-order modes lead to rapidly cycling intermittent contact. The conclusion reached is that to derive an accurate model one needs to avoid the impact limit altogether, and take sufficiently many modes in the formulation to match the actual stiffness of the constraining stop.

Introduction

Mechanical systems with repeated impacts arise in a wide variety of applications. Examples include atomic force microscopy (AFM), radio-frequency switches, Braille printers, impact forging, percussive moling, backlash in gears, freeplay in control surfaces and numerous vibro-impact problems involving buzz, squeak and rattle. See for example [START_REF] Bernardo | Piecewise-smooth dynamical systems; theory and applications[END_REF][START_REF] Brogliato | Impacts in Mechanical Systems -Analysis and Modelling[END_REF] for general theory and numerous examples. A particular motivation for this work is an attempt to understand the complex dynamical properties of tapping mode AFM as observed, e.g. in [START_REF] Payton | Feedbackinduced instability in tapping mode atomic force microscopy: theory and experiment[END_REF][START_REF] Raman | Cantilever dynamics in atomic force microscopy[END_REF][START_REF] Zhao | Characterization of intermittent contact in tapping-mode atomic force microscopy[END_REF][START_REF] Stark | Bistability, higher harmonics, and chaos in afm[END_REF][START_REF] Melcher | Multiple impact regimes in liquid environment dynamic atomic force microscopy[END_REF]. Here, a flexible micro-scale cantilever is excited at its clamped end so that a narrow tip at the free end can repeatedly impact with a sample in order to capture surface properties. Fundamental questions arise as to what is the "correct" rational-mechanics formulation of the impact of a continuously flexible structure in order to faithfully capture its post-impact dynamics.

A commonly used formalism is that of rigid-body mechanics in which impact is instantaneous and lossy, with the energy dissipation being captured by a Newtonian coefficient of restitution. The study of such single-degreeof-freedom, impact oscillators goes back to the Eastern European literature beginning in the 1950s, see e.g. [START_REF] Feigin | Forced oscillations in systems with discontinuous nonlinearities[END_REF][START_REF] Babitskii | Theory of Vibroimpact Systems. Approximate methods[END_REF] and references therein, culminating in the work of Peterka in the 1970s [START_REF] Peterka | Theoretical analysis of n-multiple (1/n)-impact solutions[END_REF]. The theory was rediscovered in the West a decade or so later through the work of Shaw & Holmes [START_REF] Shaw | Periodically forced linear oscillator with impacts: Chaos and long-period motions[END_REF] and by Michael Thompson, to whom this volume is dedicated, together with his student Ghaffari [START_REF] Thompson | Chaotic dynamics of an impact oscillator[END_REF]. Much information is now understood about the influence of grazing and chattering on the global dynamics, see e.g. [START_REF] Chillingworth | Discontinuity geometry for an impact oscillator[END_REF] and [START_REF] Bernardo | Piecewise-smooth dynamical systems; theory and applications[END_REF]Ch.1] for reviews.

Experimental investigation of impact oscillators at the macro scale, has shown qualitative agreement with the theory, see e.g. [START_REF] Piiroinen | Chaos and periodadding; experimental and numerical verification of the grazing bifurcation[END_REF][START_REF] Wagg | An experimental study of the impulse response of a vibro-impacting cantilever beam[END_REF][START_REF] Emans | Cumulative effect of structural nonlinearities: chaotic dynamics of cantilever beam system with impacts[END_REF][START_REF] Oestreich | Analytical and experimental investigation of an impact oscillator[END_REF][START_REF] Thompson | An experimental study of low velocity impacts[END_REF] (the latter being a different Michael Thompson, albeit a former PhD student from the same UCL group!) However, precise quantitative agreement is hampered by the fact that Newtonian restitution is a gross approximation to the true process of energy dissipation. In reality, dissipation occurs through excitation of internal degrees of freedom of the objects undergoing impact. In the case of a clamped cantilever undergoing intermittent contact at its free end, this problem is particularly acute because of the propensity of each impact to excite higher-order modes of the structure.

Turning to tapping mode AFM, an often adopted modality is that of dynamic force spectroscopy in which the cantilever's deflection is used as a proxy for the interaction force between the tip and the sample. Then, using a model-based approach, the signal is inverted to reconstruct the tip-sample force as a function of the tip displacement, see e.g. [START_REF] Stark | Bistability, higher harmonics, and chaos in afm[END_REF][START_REF] Van Eysden | Frequency response of cantilever beams immersed in compressible fluids with applications to the atomic force microscope[END_REF] and references therein. Material properties of the sample, e.g. its stiffness, viscosity, etc., can then in theory be inferred from the data by fitting the reconstructed, tipsample forces to a model. Although such a model-based inversion technique can be used to capture single coefficients like that of Newtonian restitution, such an approach typically does not take into account the response of higher modes of the cantilever beam. Such higher modes can cause multiple-tap dynamics, where a single impact is followed by several others [START_REF] Melcher | Multiple impact regimes in liquid environment dynamic atomic force microscopy[END_REF] like a drummer's stick during a drum roll. Such dynamic behaviour and the question of how to model it, is the chief motivation for this paper.

To be concrete, we shall consider the problem of a cantilever that is subject to spatially distributed forcing and either viscous or material damping, whose free end is allowed to make point contact with a stiff constraint. The central question is how to faithfully model such a situation using finitedimensional impact mechanics. In particular, in such problems one often appeals to two separate limiting processes to derive a tractable reducedorder mathematical model for such a configuration. The first process is that of taking the infinite stiffness limit, to avoid the numerical problems associated with different time scales of contact and free motion. Indeed, in the single-degree-of-freedom case Nordmark [START_REF] Nordmark | Grazing conditions and chaos in impacting systems[END_REF] has shown theoretically and via simulation that the dynamics of a bilinear oscillator converges to that of the corresponding impact oscillator in the limit that one stiffness becomes infinite. The second limiting process is modal truncation. For a linear beam, a finite collection of modes can capture the response exactly, provided the forcing is entirely within the subspace spanned by those modes. But impact introduces strong nonlinearity, of a distributional kind, and one cannot appeal to such standard modal convergence theorems.

The rest of this paper is organised as follows. §2 introduces the model configuration to be studied. It is shown how, given a compliant constraint, modal decomposition and Galerkin projection can be used to formulate a reduced-order model up to any order. §3 goes on to study the same model in the rigid impact limit, with a given coefficient r. Two proposed methods from the literature for posing an impact law for a given order modal truncation are presented. §4 then revisits such impact laws by taking a succession of truncated models from §3 and passing to the large stiffness limit. It is shown that this process does not converge. In §5 it is shown via numerical simulation that the observed dynamics of these impacting N -mode truncations depend strongly on the choice of N . §6 then presents experimental results for an impacting cantilever and shows how to optimally capture the observations using reduced order models. Finally, §7 draws conclusions and considers practical implications for AFM operation.

the dimensionless equation

w tt (s, t)+2αw t (s, t) + w ssss (s, t) = f (s, t) (1) s = 0 : w = w s = 0 s = 1 : w ss = 0, w sss = 0, w d 0 k 4 (d 0 -w), w > d 0 .
Here w is the scaled transverse deflection of the beam, s ∈ [0, 1] is an axial co-ordinate, with s = 0 representing the clamped end and s = 1 the free end. Subscripts s or t on w denote partial differentiation with respect to that variable, 2α is a coefficient of mass-proportional viscous damping and f (s, t) represents an axially-distributed, time-varying loading. Note that the beam's bending stiffness and mass density have both been scaled to unity in this formulation, without loss of generality. The presence of a compliant, one-sided amplitude constraint at the free-end of the cantilever is captured via the piecewise-linear shear boundary condition, where d 0 is the equilibrium distance between the free-end of the cantilever and the onesided constraint, which (for convenience of later calculations) is assumed to have a scaled stiffness k/4. The equation of motion (1) represents a piecewise-smooth boundaryvalue problem (BVP) that can be considered as a nonlinear dynamical system on an appropriate space of functions defined on the interval [0, 1]. The only nonlinearity is the Heaviside discontinuity in the boundary conditions that occurs when w(1, t) = d 0 . But if we consider motion that does not cross this constraint, then (1) represents a linear, self-adjoint equation. In what follows we shall let the subscript + denote the region w(1, t) ≥ d 0 where the beam is constrained by the stop anddenote the region w(1, t) < d 0 where the cantilever is unconstrained.

Natural frequencies for the unforced, undamped problem can be obtained in the + andregions separately by restricting motion to the relevant region, setting α = f = 0 in (1) and looking for solutions of the form w(s, t) = ψ ± i (s)e iωt . This produces a countably infinite set of mutually orthogonal eigenfunctions,

ψ ± i (s) = cos(β ± i s) -cosh(β ± i s) -σ ± i sin(β ± i s) -sinh(β ± i s), (2)
where

σ ± i = cos β ± i + cosh β ± i sin β ± i + sinh β ± i , with β ± i being the ith positive solution to cos β cosh β + 1 + k ± 4β 3 (sin β cosh β -cos β sinh β) = 0. (3)
Here k ± is a bimodal stiffness, with k -= 0 and k + = k, and (β ± i ) 2 represents the ith eigenfrequency of the linear problem.

In the expression (3) the terms inside and outside the parentheses correspond to the characteristic equations for the eigenfrequencies of a clampedpinned beam and a simple cantilever beam, respectively. The factor k/(4β 3 ) represents a weighting between the two characteristic equations. If k ≪ (β - i ) 3 , then expanding the solution of (3) in the + region for small values of k about β - i we find

β + i 4 = β - i 4 + k + O(k 2 ), ( 4
)
where k is the modal stiffness of the constraint for all of the free cantilever modes and O denotes the order as k → 0. Equation [START_REF] Bernardo | Piecewise-smooth dynamical systems; theory and applications[END_REF] shows that the firstorder effect of k is to stiffen the modes proportionally without affecting the unconstrained eigenfunctions. On the other hand, if k ≫ β - i , then β + i tends to the ith root of the characteristic equation of a clamped-pinned beam. Nevertheless, for all finite k, there exists an i that is sufficiently large so that k ≪ β - i ≈ β + i .

Modal truncation

Given the solution (2) to the eigenvalue problem for (1) in the + andregions, it is tempting to pursue development of a reduced-order, lumpedparameter model using a truncated eigenvalue expansion in each region separately. For a purely linear problem, we know by the eigenfunction expansion theorem [START_REF] Meirovitch | Principles and Techniques of Vibrations[END_REF] that such a procedure will converge. Here however, there is a subtlety because the overall problem is non-linear, and so one should not expect a purely linear eigenfunction expansion to correctly capture the solution.

To resolve this apparent paradox, suppose the exact solution within each region, written in the form

w(s, t) = ∞ i=1 ψ ± i (s)q ± i (t),
is decomposed into a truncated series approximation and a residual error:

w N ± (s, t) = N i=1 ψ ± i (s)q ± i (t), e N ± (s, t) = ∞ i=N +1 ψ ± i (s)q ± i (t).
Indeed, this was the approach adopted by [START_REF] Moon | Chaotic vibrations of a beam with nonlinear boundary conditions[END_REF] where the first mode of the cantilever and the clamped-pinned beam were retained in a lumpedparameter, single degree-of-freedom, reduced-order model for the dynamics of a cantilever in the presence of a rigid, one-sided constraint. The orthogonality relations of the modes ensure only that w N + (s, t) is orthogonal to e N + (s, t) and w N -(s, t) orthogonal to e N -(s, t). But this is not the case for w N ± (s, t) and e N ∓ (s, t). Consequently, at each switching event between + andregions, a portion of the exact solution (the residual) is lost when changing basis from the unconstrained to the constrained eigenfunctions or vice versa. Thus, when crossing from the + region to the -, some of the information in w N + is projected onto w N -when we change basis, and vice versa when crossing fromto +. Hence, if simulating for a long time it is not clear a priori that this form of modal truncation will converge as N → ∞, as errors can accumulate upon repeated crossing of the switching condition. In what follows, we shall refer to the eigenbasis {ψ + i } as the set of constrained cantilever modes and {ψ - i } as the unconstrained cantilever modes. To simplify notation we shall let φ i = (-1) (i+1) ψ - i where the sign of φ i is chosen conveniently so that φ i (1) = 2 for all i.

A more tractable approach for modal truncation, for which the study of convergence is more straightforward, is to use a single eigenbasis. That is, a reduced-order model is obtained through a Galerkin decomposition using the unconstrained cantilever eigenfunctions {φ i (s)} as the trial basis. To do this, it is natural to incorporate the interaction force between the beam and the constraint directly in the PDE (1) as an additional forcing term rather than through the boundary conditions. This can be achieved as follows:

w tt (s, t) + 2αw t (s, t)+w ssss (s, t) = δ(s -1)h (w(1, t)) + f (s, t) (5) s = 0 : w = w s = 0 s = 1 : w ss = w sss = 0, where h(d) = 0, d ≤ d 0 -k 4 (d -d 0 ) d > d 0 . (6)
Now, we make a global decomposition

w(s, t) = ∞ i=1 φ i (s)η i (t) (7)
where the φ i are unconstrained eigenfunctions which are scaled such that orthogonality relation

1 0 φ i (s)φ j (s)ds = δ ij
is enforced with δ ij being the Kronecker delta function. Substitution of ( 7) into ( 5) and taking the inner-product with φ j yields

ηj (t) + 2α ηj (t) + k j η j (t) + φ j (1) k 4 ∞ i=1 φ i (1)η i (t) = kd 0 4 φ j (1) + f j (t), (8)
where

f j (t) = 1 0 φ j (s)f (s, t)ds, and k j = (β - j ) 4 ,
and a dot is used to denote differentiation in time. Letting

η(t) = [η 1 (t), η 2 (t), η 3 (t), . . . , η N (t)] T and φ(s) = [φ 1 (s), φ 2 (s), φ 3 (s), . . . , φ N (s)] T
an N degree of freedom (NDOF) reduced-order model may be expressed as

η(t) + 2α η(t) + K ± η(t) = F ± (t), (9)
where

F -(t) = [f 1 (t), f 2 (t), f 3 (t), . . . , f N (t)] T , F + (t) = F -(t) + kd 0 4 φ(1), K - ij = δ ij k j , and K + ij = δ ij k j + k,
The switching surface between the + andregions is now given by φ T (1)η(t) = d 0 .

Modal convergence

Figure 2 shows how the eigenvalues λ 1 , λ 2 , . . . , λ N of the stiffness matrix [K + ] approximate the modal stiffnesses (β + 1 ) 4 , (β + 2 ) 4 , . . . , (β + N ) 4 of the original continuum model. Notice how, for each value of N , the highest-order mode diverges from the exact solution as k → ∞ Recall that each β + i of the exact solution tends to the ith mode of the clamped-pinned beam when k ≫ (β + i ) 3 . In the reduced-order model however, as k becomes large, the highest mode corresponds to in-phase contributions of the unconstrained cantilever modes with eigenvector

v N = N -1/2 , N -1/2 , . . . , N -1/2
and eigenvalue λ N = kN . Therefore, if k exceeds O(β + N ) 3 , the highest eigenfrequency of the truncated model and the exact eigenfrequency diverge.

This non-convergence between eigenvalues of the truncated model and the BVP in the limit k → ∞ is in fact an unavoidable consequence of performing Galerkin discretization using the unconstrained eigenmodes, which satisfy natural (force and moment free) boundary conditions at the free end. In the limit that the stiffness of the constraint becomes large compared to that of the highest-order included eigenfunction, the stiffness of the constraint at the free end makes the boundary condition much closer to a geometric (pinned) one. Therefore, the unconstrained cantilever modes do not provide a good approximation to the true dynamics. Nevertheless, for a given constraint stiffness k, choosing N such that (β + N ) 3 ≫ k ensures that the modal description of the reduced-order model will accurately capture that of the full continuum BVP. In particular, for fixed k, we have modal convergence as N → ∞.

The impact limit; modal restitution laws

A common approach in modelling the dynamics of a flexible beam with a stiff, one-sided constraint has been to consider the so-called impact limit where the constraint becomes rigid. This approach imposes the classical Newtonian coefficient of restitution (CoR) law at a localized point on the continuous beam given by w t (s I , t I+ ) = -rw t (s I , t I-), w(s

I , t I ) = d 0 , (10a) w t (s, t I+ ) = w t (s, t I-), s = s I , (10b)
where 0 < r < 1 is the coefficient of restitution, t I± is the time instantaneously after (+) or before (-) an impact occurs and s = s I is the location of the impact. In what follows, we shall assume that s I = 1, so that impact occurs at the free end of the beam. Incorporating ( 10) into a modal model has the general form η(t

I+ ) = -[R] η(t I-), (11)
where [R] is a modal CoR matrix which allows the modal energy to be redistributed post impact. An NDOF, reduced-order, impact model is established once the modal CoR matrix [R] is determined .

via such an approach is simply

[R] C =        φ T (s 1 ) φ T (s 2 ) . . . φ T (s N -1 ) φ T (s I )        -1        -1 0 . . . 0 0 0 -1 . . . 0 0 . . . . . . . . . . . . . . . 0 0 . . . -1 0 0 0 . . . r               φ T (s 1 ) φ T (s 2 ) . . . φ T (s N -1 ) φ T (s I )        . (12)
Note, as pointed out in [START_REF] Vyasarayani | Modeling impacts between a continuous system and a rigid obstacle using coefficient of restitution[END_REF], that [R] C is non-unique due to the arbitrariness in the choice of the s i . Nevertheless, if the s i are chosen to be equally spaced along the beam, then, in the limit as N becomes large, one might expect [R] C to converge, making this approach somewhat analogous to finite element analysis.

Impulse-response restitution law

An alternative method for determining the modal CoR, proposed in [START_REF] Vyasarayani | Modeling impacts between a continuous system and a rigid obstacle using coefficient of restitution[END_REF], is to consider the impulse response of the reduced-order model and determine the magnitude of the impulse required to enforce (10)a. To obtain such a matrix, consider a cantilever subjected to an isolated impulsive point force

w tt (s, t) + w ssss (s, t) = p 0 δ(s -s I )δ(t), (13)
where p 0 is the magnitude of the impulse. Pursuing a solution to (13) via eigenfunction expansion involves writing w(s, t) = ∞ i=1 ψ - i (s)q - i (t) and truncating the series after N terms. Choosing s I = 1, projecting the solution onto the jth mode and solving for q - j (t) yields q - j (t) = p 0 ψ - i (1)(β - j ) -2 sin(β - j ) 2 t, t > 0. Substitution of the truncated modal expansion for the impulse response into (10) solves for the unknown p 0 and establishes the relationship qj (t

I+ ) = q- j (t I-) - r + 1 N N i=1 q- i (t I-).
Following this approach, the impulse-response

CoR [R] = [R] IR for the 2DOF, reduced-order model is given by [R] IR = 1 N      r + 1 -N r + 1 . . . r + 1 r + 1 r + 1 -N . . . r + 1 . . . . . . . . . r + 1 r + 1 r + 1 r + 1 r + 1 -N      (14)
4 Restitution law via infinite stiffness limit

Rather than impose the CoR law in [START_REF] Moon | Chaotic vibrations of a beam with nonlinear boundary conditions[END_REF] a priori, we shall now seek to derive a general form of CoR by considering the modal compliant model in the limit that the constraint stiffness k → ∞.

To this end, suppose we let k = ε -2 κ, where 0 < ε ≪ 1 is a scaling parameter and κ = O(1) as ε → 0. In this case it is straightforward to see that the time taken for a trajectory of the dynamical system to pass through the + region is O(ε -1 ). Thus, if we take the limit ε → 0 + , this time tends to zero and we will have an equivalent impact oscillator with a modal impact law [START_REF] Nordmark | Grazing conditions and chaos in impacting systems[END_REF]. The question arises how the impact law [R] N derived by this asymptotic approach applied to the compliant N -modal truncation will correspond to either

[R] C or [R] IR . The convergence properties of [R] N as N → ∞ are also of interest.
In order to answer these questions, we shall consider separately the cases N = 1 and N = 2, before tackling the general case. Finally, to simplify the analysis, we consider the special case of d 0 = 0.

Single-degree-of-freedom impact oscillator

Consider the single degree-of-freedom (1DOF) compliant model with initial conditions on the switching surface and impending motion in the + region:

η1 (t) + 2α η1 (t) + (k 1 + k)η 1 (t) = f 1 (t), η 1 (t) ≥ 0, ( 15
) η 1 (0) = 0, η1 (0) = ε -1 v + > 0.
Letting k = ε -2 κ, t = ετ and ξ = η 1 , we can rewrite [START_REF] Peterka | Theoretical analysis of n-multiple (1/n)-impact solutions[END_REF] as

ξ τ τ + 2αεξ τ + κξ = O(ε 2 ), ξ > 0 (16) ξ(0) = 0, ξ τ (0) = v + .
Substituting ξ = ξ 0 + εξ 1 + . . . into [START_REF] Piiroinen | Chaos and periodadding; experimental and numerical verification of the grazing bifurcation[END_REF] and collecting terms by order of ε yields

O(ε 0 ) : ξ 0 τ τ + κξ 0 = 0, O(ε 1 ) : ξ 1 τ τ + κξ 1 = -2αξ 0 τ , which has solution ξ(τ ) = v + κ -1/2 (1 -εατ ) sin κ 1/2 τ + O(ε 2 ), 0 ≤ τ ≤ ∆ + , (17)
where

∆ + = πκ -1/2 (18)
is the time of passage in the + region. Differentiating [START_REF] Raman | Cantilever dynamics in atomic force microscopy[END_REF] with respect to time and substituting τ = ∆ + , we find

v -= -r e v + , with r e = 1 -απκ -1/2 ε (19)
where v -is the velocity when the oscillator re-enters theregion. We therefore call r e the effective coefficient of restitution. Figure 3(a) shows a comparison between the return velocity v -predicted by the impact law [START_REF] Shaw | Periodically forced linear oscillator with impacts: Chaos and long-period motions[END_REF] and that obtained by simulating the 1DOF model [START_REF] Peterka | Theoretical analysis of n-multiple (1/n)-impact solutions[END_REF] for different values of the constraint stiffness. Specifically, letting k = ε -2 , η 1 (0) = 0 and η1 (0) = ε -1 , the effective restitution coefficient r 2 e is compared to the computed (v -/v + ) 2 , as ε is varied. Note from the figure, the close agreement for small ε. Also, note r e → 1 in the limit as ε → 0 + . However, in experiments, Thompson et al [START_REF] Thompson | An experimental study of low velocity impacts[END_REF] observed low values of the coefficient of restitution that were attributed to energy transfer to higher modes during contact. Hence, we will next consider the impact limit of the 2DOF model.

Two degrees of freedom

Next, we consider the effect of an additional mode. For a 2DOF reducedorder model to the + region, the governing equation reduces to

η(t) + 2α η(t) + K + η(t) = F -(t), φ T (1)η(t) > 0. (20)
Substituting k = κε -2 /2 into [K + ] and scaling by ε 2 , we find the eigenvectors are

[V ] = 1 √ 2 1 1 -1 1 + O(ǫ 2 ). ( 21
)
Letting η(t) = [V ] q(t), q(t) = [q 1 (t), q 2 (t)] T and t = ετ , (20) can be expressed as

q τ τ (τ ) + 2αεq τ (τ ) + 0 0 0 κ q(τ ) = O(ǫ 2 ), q 2 (τ ) > 0. ( 22
)
Notice that q 1 is orthogonal to the switching surface. The lowest-order mode imitates a clamped-pinned beam mode with zero amplitude at the free end. Letting ξ = q 2 we recover [START_REF] Piiroinen | Chaos and periodadding; experimental and numerical verification of the grazing bifurcation[END_REF] for the 1DOF oscillator. It follows that the velocity of q 2 obeys the impact law given in [START_REF] Shaw | Periodically forced linear oscillator with impacts: Chaos and long-period motions[END_REF]. Additionally, integrating the equation for q 1 and expanding in powers of ε, we relate q 1τ (0) = v + 1 to the return velocity q 1τ (∆

+ ) = v - 1 by v - 1 = v + 1 1 -2απκ -1/2 ε + O(ε 2 ). ( 23
)
Thus, the modal coefficient of restitution is given by

[R] 2 = [V ] 1 -2r e 0 0 r e [V ] T . (24)
Substituting [V ] from ( 21) into (24) we find

[R] 2 = 1 2 r e -1 1 -3r e 1 -3r e r e -1 . (25)
Notice that [R] 2 deviates in the off-diagonal terms from the expression ( 14) for the impulse-response restitution law [R] IR . The source of the discrepancy is that the energy loss captured in r e is from the external viscous damping, whereas r in ( 14) captures the ad-hoc energy dissipation modelled by the Newtonian impact law. We note that both converge in the limit as r and r e tend to unity.

The dynamics of the 1DOF and 2DOF impact models are fundamentally different post impact. To see this, consider pre-impact, initial conditions η(0) = 0, and η(t

) = v + = [ε -1 v + 1 , 0] T .
The oscillator returns to theregion with a velocity

v -= ε -1 v - 1 , v - 2 T , where v - 1 = -R 11 v + 1 , v - 2 = -R 21 v + 1 ,
and R 11 and R 21 are the elements from [R] 2 . Letting k = ε -2 , the postimpact velocities predicted by the impact law are compared to numerical simulations of the reduced-order, compliant model as ε is varied in Fig. 3(b,c). In the case of the 2DOF model, the effective coefficient of restitution of the first mode tends to zero as ε → 0 + . In other words, the motion of the first mode of the 2DOF reduced-order model in the impact limit becomes completely uncorrelated with the motion of the equivalent 1DOF model after the first impact.

General N degree-of-freedom case

For the N degree-of-freedom case, let κ = kN ε 2 . Following the approach for the 2DOF case, we substitute q(t) = [V ] T η(t) to diagonalize [START_REF] Melcher | Multiple impact regimes in liquid environment dynamic atomic force microscopy[END_REF] in the + region. Without solving for either [V ] or [Λ], we know that λ j = O(ε 2 ) for j < N and λ N = κ + O(ε 2 ) and the eigenvectors are mutually orthogonal with the eigenvector of the high-frequency mode given by

v N = N -1/2 , N -1/2 , . . . , N -1/2 T . ( 26
)
It follows that the equation of motion in the + region can be expressed as

q τ τ (τ ) + 2αεq τ (τ ) +      0 0 . . . 0 0 0 . . . 0 . . . . . . . . . 0 0 0 0 κ      q(τ ) = O(ε 2 ), q N (τ ) > 0. ( 27
)
Substituting ξ = q N into the equation for the N th mode in [START_REF] Wagg | A note on coefficient of restitution models including the effects of impact induced vibration[END_REF], we again recover [START_REF] Piiroinen | Chaos and periodadding; experimental and numerical verification of the grazing bifurcation[END_REF]. Accordingly, the time of passage in the + region and the return velocity of the N th mode are again given by ( 18) and ( 19), respectively. The equations for q j , j < N are again over-damped oscillators for which the return velocities are predicted by [START_REF] Thompson | An experimental study of low velocity impacts[END_REF]. It follows that the NDOF modal coefficient of restitution is given by

[R] N = [V ]      1 -2r e 0 . . . 0 0 1 -2r e . . . 0 . . . . . . . . . 0 0 0 0 r e      [V ] T . ( 28
)
From ( 26) and given that [V ] is an orthogonal matrix, it is possible to show

[R] N =      1 -2r e + 3re-1 N 3re-1 N . . . 3re-1 N 3re-1 N 1 -2r e + 3re-1 N . . . 3re-1 N . . . . . . . . . 3re-1 N 3re-1 N 3re-1 N 3re-1 N 1 -2r e + 3re-1 N .      (29)
Substituting N = 1 and N = 2 recovers the results for the 1DOF and 2DOF oscillators, respectively.

Modal non-convergence as N → ∞

Modal non-convergence in the impact limit can be demonstrated readily with the NDOF modal CoR given in [START_REF] Wagg | An experimental study of the impulse response of a vibro-impacting cantilever beam[END_REF]. The post-impact velocity profile is given by

w N t (s, t I+ ) = -φ T (s) [R] N v -. ( 30
)
For the truncated model to converge, w N t (s, t + I ) must converge for a given v -as N → ∞. However, increasing N in the [R] N excites higher-order modes to the same degree as the lower-order modes.

Fig. 3(d) shows the truncated approximation for the post-impact velocity distribution w N t (s, 0 + ) for N = 1, 2, 3, 10, r e = 1 and v -= [1, 0, 0, ..., 0] T .

Clearly, the post-impact velocity does not converge as modes are added to the restitution law. This non-convergence comes from taking an N -modal model, passing to the impact limit ε → 0, and then considering what happens as N → ∞. It would appear that these two separate limiting processes are incommensurate.

To gain further insight, it is instructive to consider again the results in Fig. 2. This shows that modal convergence in the impact limit is achieved only if the limits ε ∼ k -1/2 → 0 and N → ∞ are taken together. That is, for each N there is a maximum value of k = O(β + N ) 3 for which the model gives answers that are commensurate with other modal truncations. However, modelling the interaction with the constraint with an impact law also requires k ≫ O(β + N ) 4 . At least for a uniform cantilever beam, the two conditions are incompatible. Thus to accurately model the interaction between the cantilever and the constraint in a truncated model requires a compliant formulation.

Numerical studies

In the previous section we found that for the eigenfrequencies of the reducedorder model to converge to the exact solution of the continuous BVP required k ≪ O(β + N ) 3 . However, since the system is non-linear, this requirement is only a necessary condition for convergence. As k becomes large the excursion of the oscillator into the + region is characterised by increasingly small length scales where contributions from higher-order modes may not be considered negligible. To study convergence of the non-linear system we turn to numerical simulations. Numerical simulations using the model studied so far led to qualitatively similar results to those presented in what follows (although there was a quantitative difference in the amount of damping seen in the higher modes). However, to get an acceptable match with experimental results presented in the next Section, we found it necessary to include material damping in addition to the viscous damping from the surrounding fluid.

Inclusion of material damping

The simplifying assumption in the previous sections of linear, mass-proportional damping enabled two sets of uncoupled modes corresponding to the + andregions, which made the analysis tractable. If we now add material damping to [START_REF] Emans | Cumulative effect of structural nonlinearities: chaotic dynamics of cantilever beam system with impacts[END_REF] we have

w tt (s, t) + 2αw t (s, t) + γw sssst (s, t) + w ssss (s, t) = δ(s -1)h (w(1, t)) + f (s, t) (31)
s = 0 : w = w s = 0 (32) s = 1 : w ss = w sss = 0 ( 33
)
where γ is the coefficient of stiffness-proportional material damping, and h(•) is given [START_REF] Feigin | Forced oscillations in systems with discontinuous nonlinearities[END_REF]. Again, using a trial basis φ i (s) corresponding to the normal modes of the simple cantilever, the NDOF Galerkin model becomes

η(t) + [C] η(t) + K ± η(t) = F ± (t), (34)
where C ij = (2α + γk j )δ ij is the modal damping attenuation and K ± and F ± have the same meaning as in [START_REF] Melcher | Multiple impact regimes in liquid environment dynamic atomic force microscopy[END_REF]. While in theregion, the damping operator can be written as a linear combination of the mass and stiffness operators. However, in the + region, we find full coupling, which renders the perturbation analysis in the previous section less tractable. Nevertheless we should expect similar conclusions.

On the other hand, ( 9) and (34) are equally amenable to numerical analysis. Representing (34) in state-space form

ẋ = A ± x + F± (t) (35)
where

x(t) = η(t) η(t) , A ± = [0] [I] -[K ± ] -[C] , F± (t) = 0 F ± (t)
The solution of (35) is

x ± (t) = exp [A ± ](t -t 0 ) x ± 0 + t t 0 exp [A ± ](t -t ′ ) F± (t ′ )dt ′ , (36)
where x ± 0 = x 0 (t 0 ) is an initial condition either on the switching surface or, or in the case that f (t) is periodic with perioid ω, on the stroboscopic Poincaré section {ωt = 0 mod 2π}. In the latter case, t 0 can be reset to 0 without loss of generality. Note that the convolution integral in (36) can be evaluated by parts.

Simulation via root finding

The difficulty in using the explicit solution (36) for simulation arises arises when solving for the times-of-passage ∆ ± in the + andregion, which are the roots of a transcendental equation. These can be found using a numerical root-finding algorithm. To ensure that near-grazing events with the switching surface are detected, a golden-section search method is used to identify local extrema in the trajectory of d(t) on the bounded interval from [t 0 , 2π/ω]. The trajectory is then segmented such that d(t) is monotonic on each segment with at most one root. The bisection method is used to find the root of the first segment to cross the switching surface.

For the numerical simulations we choose the parameters d 0 = 0, α = 0 and define ζ 1 = γ β - 1 -2 /2 = 0.01, which is the damping ratio of the fundamental mode. Additionally, we choose the specific form of harmonic forcing given by

f (x, t) = k 1 ψ 1 (s) 2ζ 1 ψ 1 (1) cos ωt, with ω = 2 β - 1 2 , (37)
which projects purely to the first mode, and has frequency set to the effective resonance of the fundamental mode in the impact limit [START_REF] Shaw | A periodically forced piecewise linear oscillator[END_REF]. The amplitude of f (x, t) is chosen, without loss of generality, so that the unconstrained amplitude of w(1, t) is unity when

ω = β - 1 2 .
To simulate periodic orbits, we choose initial conditions x + (0) = 0 and discard the first 1000 excitation cycles. We then say that the resulting trajectory is a period-m orbit if ||x(t 0 + 2πm/ω)x(t 0 )|| < 1 × 10 -9 . Such a periodic orbit with n impacts is denoted (m, n).

Results

Figure 4 compares numerically simulated periodic orbits for the truncated, bilinear model with N = 1, 2, 3, 10 and k = 100k 1 = O(k 2 ). In each case, the above numerical procedure resulted in a unique periodic orbit. For the 1DOF model, this is a simple (1, 1)-orbit. However, for the 2-,3-and 10DOF models there are additionally complexities due to the higher-order modes contributing significantly to the interaction between the cantilever and the constraint. Thus, the 2DOF model predicts a (2, 4) orbit, whereas the 3DOF-model has converged to a (1, 3)-orbit. Essentially the first impact excites high frequency modes which result in additional taps on the constraining surface. In contrast the 10DOF model converges to a (1, 1) orbit with a complex interaction with the constraint. While m and n are equal

Effective energetic coefficient of restitution via simulation

When conducting numerical simulations, one possible approach for deriving an equivalent impact model of low degree of freedom is to compute an effective CoR by estimating the energy lost to higher-order modes through impact. We shall consider just such an approach to define a so-called energetic coefficient of restitution r E which captures such energy lost from the fundamental mode in a 1DOF impact model -see [START_REF] Stronge | Impact Mechanics[END_REF][START_REF] Wagg | A note on coefficient of restitution models including the effects of impact induced vibration[END_REF] and references therein for related approaches in the literature. Equating the energy lost during impact of a 1DOF impact oscillator to the energy lost from the fundamental mode of the NDOF bilinear oscillator in a single impact leads to an implicit definition of r E for that impact:

1 2 1 -r2 E (v + 1 ) 2 + t i +∆ t i F c [d(t ′ )] η1 (t ′ )dt ′ = 0. ( 38
)
Here t = t i is the time instant at which passage occurs from theregion to the + region in the NDOF system, at which time the first mode has velocity v + 1 , ∆ is the time of passage in the + region, η j is jth modal coordinate of the NDOF bilinear oscillator and

F c (d) = 0 d ≤ d 0 k 4 (d 0 -d)φ 1 (1), d > d 0 ,
is the interaction force between the beam and the constraint. We can then define a single energetic coefficient r E for a specific set of parameters and initial conditions, assuming repeated impacts occur via:

1 2 1 -r 2 E ∞ i=1 v +(i) 1 2 + lim t→∞ t t 0 F c [d(t ′ )] η1 (t ′ )dt ′ = 0, (39)
where t i and v

+(i) 1
are the set of impact times and corresponding velocities. Note that there can be numerical difficulties in evaluation of (39). In the case that motion converges to a periodic orbit, the limit can be evaluated simply by taking t 0 to be a time after a transient behaviour has decayed, and evaluating the integral and sum over a single period. But, for an arbitrary set of parameters, there is no guarantee that the limit in (39) will converge, nor that the limit will be unique over different initial conditions. For instance, it is known that even 1DOF impact oscillators can undergo chattering sequences, in which there are an infinite number of impacts in a finite time, and can feature coexisting attractors for the same parameter values -see [START_REF] Bernardo | Piecewise-smooth dynamical systems; theory and applications[END_REF] and references therein.

Figure 6 shows a comparison between the dynamics of 20DOF bilinear oscillator and the 1DOF impact oscillator with energetic coefficient of restitution computed using the above approach. In both cases, the motion converges to a periodic orbit, and there is therefore a unique, computable, solution to (39).

In the case depicted in Fig. 6(a), d 0 = 0.8, there is relatively large damping and parameters are chosen so that all initial conditions converge to a (1, 1) orbit that is close to grazing. In this case, the computation (39) predicts r E = 0.036 and there is a good quantitative agreement between the 20DOF model and the effective 1DOF impact oscillator. In this case, only a small portion of the total energy stored in the oscillator is redistributed to higher modes during impact which is subsequently dissipated as the damping is sufficiently high. Nevertheless, r E = 0.036 is significantly lower than that which would be predicted from a naive single-degree-of-freedom model using a typical value of 0.5-0.8 for the CoR for macroscopic collisions of steel on steel.

The case depicted in Fig. 6(b) is more like that in Fig. 4, in which d 0 = 0, and the damping is relatively weak. Here, (39) predicts r E = 0.74, which is closer to what one would expect from macroscopic impact theory. However, the results show that there is neither good quantitative nor qualitative agreement between the 1DOF impact model and the 20DOF bilinear model. Two particular effects that are due to the higher-order modes are the 'multiple-tap' dynamics of the impact where there are two separate intervals in theregion, and the large additional oscillations in the dynamics which are not damped out during the free motion away from the stop.

Experiments

In this section we present the results obtained from the experimental rig depicted in Fig. 7. A mild steel cantilever beam with dimensions 300 mm × 25.5 mm × 1 mm was mounted on a shaker table. A pin constraint was also mounted on a second cantilever beam with dimensions 300 mm × 25.5 mm × 12 mm. From the dimensions of the cantilever and the constraint, we estimate k/k 1 = O(10 3 ). The natural frequency ω 1 = 2π × 5.92 rad/s and damping ratio ζ 1 = 7.7 × 10 -4 of the fundamental mode was determined from a ring-down test. Velocity measurements at the point on the beam corresponding to the location of the constraint was made with a Polytec OFV-5000 laser Doppler vibrometer (LDV) in the ground frame. However, due to the low damping, the motion of the base is much smaller than the free-end of the beam and can be neglected. The velocity signal from the LDV is sampled at 50 kHz and has a resolution of 4 µm/s.

In addition to the velocity measurement, a simple electronic circuit was used to accurately determine when the stop and the beam are in contact. The circuit consists of two resistors in series which are supplied with a constant voltage V 0 . The cantilever beam and constraint act as a switch that shorts the first resistor when the beam is in contact with the constraint. The voltage over the second resistor was sampled at 100 kHz.

To excite the cantilever the shaker table is used to prescribe a harmonic motion to its base. The motion of the base imparts a uniformly distributed inertial force on the cantilever. The inertial force projects onto each eigenmode, however, the effect on higher eigenmodes is minimal because (i) the excitation frequency is closest to the fundamental mode and (ii) the higher modes are substantially stiffer than the fundamental. Consequently, the forcing from the base excitation of the physical system can be approximated by (37).

Figure 8(a) shows the results of a particular experimental run over one period T = 2π/ω of the forcing, obtained after transients had been allowed to decay. The motion was found to be approximately periodic, and the presented results are typical of what happens on each period. Approximately, when viewed at a macro scale, this is a (1, 1) periodic orbit. However, the detailed time series (shown in two successive zooms in Fig. 8(b,c)) reveal the motion actually to be more highly structured and to contain five separate time intervals in which the beam is in contact with the stop. There is also evidence of high-frequency, undamped motion between each of the impacts. It would be tempting to wonder if the high-frequency motion is due to measurement noise. However, a comparison between the excited motion of the beam without the stop in Fig. 8(d) reveals that the high-frequency response is several orders of magnitude above the noise threshold.

Figure 9 shows simulation results for the 20DOF bilinear model at the same parameters values as those estimated from the experimental data. Note the strong qualitative similarity with the experimental data; that there is similar multiple-tap dynamics with five separate impacts, and significant energy in the higher modes.

Acknowledgements

The authors acknowledge useful conversations with Arvind Raman, Oliver Payton and Chris Budd. The work of JM has been sponsored by the UK Engineering and Physical Sciences Research Council.


Marcio Gouveia

Jaume Llibre

Douglas D Novaes

Claudio Pessoa

Márcio R A Gouveia

Piecewise smooth dynamical systems: Persistence of periodic solutions and normal forms

Keywords: Piecewise differential system, Crossing periodic solutions, Limit cycles, Normal forms, Lyapunov-Schmidt reduction

We aim to study the persistence of periodic orbits contained in A when (1) is perturbed inside the class of all piecewise smooth vector fields with two zones separated by the hyperplane . A similar problem was addressed in [START_REF] Du | Bifurcation of periodic orbits in a class of planar Filippov systems[END_REF] but for planar piecewise smooth systems. In this case the authors developed the Poincaré map in power series of some small parameter, and then a Melnikov-like function was obtained just by taking the first coefficient of this series. In dimension n ≥ 3 the arguments are more delicate, and some recent techniques (see [START_REF] Llibre | Improving the averaging theory for computing periodic solutions of the differential equations[END_REF]) on Lyapunov-Schmidt reduction have to be used. For instance, in [START_REF] Liu | Poincaré bifurcation of a three-dimensional system[END_REF] the authors have used similar techniques to study the same problem but for 3D smooth systems.

In smooth systems the normal form theory provides simplifications on the systems reducing the number of parameters and, eventually, simplifying the study of their dynamics. On the other hand, for nonsmooth systems formed by two regions separated by a hyperplane, more simplifications on the parameters can be done. In fact, piecewise continuous change of variables (smooth on each region determined by the hyperplane) is used to obtain normal forms for each one of the systems constituting the nonsmooth system. Some normal forms have been obtained in [START_REF] Freire | Canonical discontinuous planar piecewise linear systems[END_REF] for planar piecewise smooth systems separated by a straight line.

Taking all the comments above into account our first main result, Theorem A, develops a Melnikov-like function to control the persistence of periodic solutions of A for small discontinuous perturbations of system [START_REF] Buzzi | Chaotic planar piecewise smooth vector fields with non-trivial minimal sets[END_REF]. Our second main result, Theorem B, provides normal forms for all 3D piecewise linear systems which admits an invariant plane transversal to containing a period annulus A fulfilled by periodic orbits. Finally, in section 4, we apply Theorem A, with n = 3, to study discontinuous perturbations of the normal forms established in section 3.

We point out the importance of nonsmooth systems, particularly piecewise smooth systems, as a source of new interesting phenomena in dynamics (see, for instance, [START_REF] Buzzi | Chaotic planar piecewise smooth vector fields with non-trivial minimal sets[END_REF][START_REF] Novaes | Shilnikov problem in Filippov dynamical systems[END_REF]), and their relevance in modelling real phenomena (see, for instance, [START_REF] Bernardo | Piecewise-Smooth Dynamical Systems: Theory and Applications[END_REF][START_REF] Coombes | Neuronal networks with gap junctions: a study of piecewise linear planar neuron models[END_REF]).

Persistence of crossing periodic orbits

Consider the following discontinuous perturbation of system [START_REF] Buzzi | Chaotic planar piecewise smooth vector fields with non-trivial minimal sets[END_REF] Z(x, y, z; ε) = Z + (x, y, z; ε) = X + 0 (x, y, z)

+ εX + 1 (x, y, z) if z > 0, Z -(x, y, z; ε) = X - 0 (x, y, z) + εX - 1 (x, y, z) if z < 0, (2)
where

(x, y, z) ∈ D ⊂ R × R n-2 × R, X j = (X j,1 , X j,2 , .
. . , X j,n ) for j = 0, 1 and the set of discontinuity is given by the hyperplane = {z = 0}. As usual the following open regions are distinguished in :

(i) Crossing region: c = {p ∈ : (Z + f )(p)(Z -f )(p) > 0}; (ii) Escaping region: e = {p ∈ : (Z + f )(p) > 0, (Z -f )(p) < 0}; (iii) Sliding region: s = {p ∈ : (Z + f )(p) < 0, (Z -f )(p) > 0}.
Here

(Z ± f )(p) = Z ± (p), ∇f (p) .
A crossing periodic orbit of system ( 9) is a closed curve γ composed by trajectories of Z ± having the same orientation such that γ ∩ ⊂ c . We recall our main hypothesis:

(H ) the hyperplane = {x = 0} has a period annulus A = {(y, z) ∈ : r 0 ≤ |(y, z)| ≤ r 1 }
surrounding the origin and fulfilled by crossing periodic solutions of the unperturbed system ( ẋ, ẏ, ż) T = Z(x, y, z; 0).

Let ϕ ± (t, x, y, z; ε) = ϕ ± 1 (t, x, y, z; ε), ϕ ± 2 (t, x, y, z; ε), . . . , ϕ ± n (t, x, y, z; ε) be the solutions of the systems ( ẋ, ẏ, ż) T = Z ± (x, y, z; ε) such that ϕ ± (0, x, y, z; ε) = (x, y, z). From hypothesis, the solutions of (2), for ε = 0, contained in A reach transversally the set of discontinuity . So for a small neighborhood U ⊂ R n of A and |ε| = 0 small enough there exists a time t + (x, y; ε) > 0 (resp. t -(x, y; ε) < 0) such that an orbit of (2) starting in (x, y, 0) ∈ U ∩ returns, forward in time (resp. backward in time), to , that is ϕ ± n (t ± (x, y; ε), x, y, 0; ε) = 0. The solutions ϕ ± (t, x, y, z; ε) can be expressed in power series of ε, that is

ϕ ± i (t, x, y, z; ε) = ψ ± 0,i (t, x, y, z) + εψ ± 1,i (t, x, y, z) + O(ε 2 ), i = 1, 2, . . . , n such that ϕ ± (t, x, y, z; ε) = ψ ± 0 (t, x, y, z) + εψ ± 1 (t,
x, y, z) where ψ ± 0 (t, x, y, z) are the solutions of the unperturbed systems ( ẋ, ẏ, ż) T = Z ± (x, y, z; 0). The times t ± , introduced above, can also be written as power series of ε, that is

t ± (x, y; ε) = τ ± 0 (x, y) + ετ ± 1 (x, y) + O(ε 2 ).
However here we only have to assume that the expressions σ ± j (y) = τ ± j (0, y), j = 0, 1, are explicitly known. In fact, the next result gives the expressions of σ ± 1 and ψ ± 1 in terms of the solutions ψ ± 0 of the unperturbed system.

Proposition 1. Let p = (x, y, z) ∈ R × R n-2 × R
, and denote by Y ± (t, p 0 ) = D p ψ ± 0 (t, p 0 ) the derivative of ψ ± 0 (t, p) with respect to the initial condition p evaluated at p 0 = (x 0 , y 0 , z 0 ). Then the following equalities hold:

ψ ± 1 (t, p) = Y ± (t, p) t 0 Y ± (s, p) -1 X ± 1 ψ ± 0 (s, p) ds, ( 3
)
σ ± 1 (y) = - ψ ± 1,n (σ ± 0 (y), 0, y, 0) X ± 0,n ψ ± 0 (σ ± 0 (y), 0, y, 0) . ( 4
)
Proof. Computing the derivative in the variable t of both sides of the equality ϕ ± (t, x, y, z; ε) = ψ ± 0 (t, x, y, z) + εψ ± 1 (t, x, y, z) + O(ε 2 ) we obtain

X ± 0 ϕ ± (t, p; ε) + εX ± 1 ϕ ± (t, p; ε) = X 0 ψ ± 0 (t, p) + ε ∂ψ ± 1 ∂t (t, p) + O(ε 2 ).
Developing the left hand side of the above equation in power series of ε, and studying the coefficient of ε we obtain

∂ψ ± 1 ∂t (t, p) = DX ± 0 ψ ± 0 (t, p) ψ ± 1 (t, p) + X ± 1 ψ ± 0 (t, p) .
Moreover ψ ± 0 (0, p) = 0. Hence the solution of the above differential equation is given by (3). To see equality (4) we first develop the equality ϕ ± n (t ± (0, y; ε), 0, y, 0; ε) = 0 in power series of ε, being t ± (0, y; ε) = σ ± 0 (y) + εσ ± 1 (y) + O(ε). After that σ ± 1 can be isolated as (4). This concludes the proof of the proposition. 2

Now let V = {ν ∈ R n-2 : (0, ν, 0) ∈ U }. We define a Melnikov-like function M : V → R n-2 as M(ν) = (ν) - λ 1 (ν) ω 1 (ν) (ν), (5)
where

λ i (ν) = X + 0,i (ψ + 0 (σ + 0 (ν), 0, ν, 0))σ + 1 (ν) -X - 0,i (ψ - 0 (σ - 0 (ν), 0, ν, 0))σ - 1 (ν) + ψ + 1,i (σ + 0 (ν), 0, ν, 0) -ψ - 1,i (σ - 0 (ν), 0, ν, 0), ω i (ν) = ∂ψ + 0,i ∂x (σ + 0 (ν), 0, ν, 0) - X + 0,i ψ + 0 (σ + 0 (ν), 0, ν, 0) X + 0,n (ψ + 0 (σ + 0 (ν), 0, ν, 0)) ∂ψ + 0,n ∂x (σ + 0 (ν), 0, ν, 0) - ∂ψ - 0,i ∂x (σ - 0 (ν), 0, ν, 0) + X - 0,i ψ - 0 (σ + 0 (ν), 0, ν, 0) X - 0,n (ψ - 0 (σ - 0 (ν), 0, ν, 0)) ∂ψ - 0,n ∂x (σ - 0 (ν), 0, ν, 0), (6)
for i = 1, 2, . . . , n -1, and

(ν) = λ 2 (ν), λ 3 (ν), . . . , λ n-1 (ν) , (ν) = ω 2 (ν), ω 3 (ν), . . . , ω n-1 (ν) .
Let J M (ν) denote the determinant of the Jacobian matrix of M evaluated at ν. The next theorem is the main result of this section.

Theorem A. In addition to hypothesis (H ) we assume that ω 1 (ν) = 0 for every ν ∈ V. Then for each ν * ∈ V such that M(ν * ) = 0 and J M (ν * ) = 0, there exists a unique crossing periodic solution φ(t, ε) of system (2) such that φ(0, ε) → (0, ν * , 0) ∈ A when ε → 0. Moreover if M(ν) = 0 for every ν ∈ V then there are no crossing periodic solutions bifurcating from A for |ε| = 0 sufficiently small.

To prove Theorem A we shall use a version (see [START_REF] Llibre | Improving the averaging theory for computing periodic solutions of the differential equations[END_REF]) of the so-called Lyapunov-Schmidt reduction for finite dimensional function (see [START_REF] Chicone | Lyapunov-Schmidt reduction and Melnikov integrals for bifurcation of periodic solutions in coupled oscillators[END_REF]). In what follows, for positive integers k < d, the functions ξ : [START_REF] Llibre | Improving the averaging theory for computing periodic solutions of the differential equations[END_REF].) Assume that k ≤ d are positive integers. Let D and V be open bounded subsets of R d and R k , respectively. Let g 0 , g 1 and β :

R k × R d-k → R k and ξ ⊥ : R k × R d-k → R d-k will denote the projections onto the first k coordinates and onto the last d -k coordinates, respectively. A point ζ ∈ R d is denoted by ζ = (a, b) ∈ R k × R d-k .

Lemma 2. (See

V → R d-k be C 2 functions, consider g : D × (-ε 0 , ε 0 ) → R d as g(ζ, ε) = g 0 (ζ ) + εg 1 (ζ ) + O(ε 2 ),
and take

Z = {ζ ν = (ν, β(ν)) : ν ∈ V } ⊂ D. We denote by ν the upper right corner k × (d -k) matrix of D g 0 (ζ ν ), and by ν the lower right corner (d -k) × (d -k) matrix of D g 0 (ζ ν ). Assume that for each ζ ν ∈ Z, det( ν ) = 0 and g 0 (ζ ν ) = 0. We define the bifurcation function f 1 : V → R k as f 1 (ν) = -ν -1 ν ξ ⊥ g 1 (ζ ν ) + ξg 1 (ζ ν ). ( 7
)
If there exists

ν * ∈ V such that f 1 (ν * ) = 0 and J f 1 (ν * ) = 0, then there exists ν ε such that g(ζ ν ε , ε) = 0 and ζ ν ε → ζ ν * when ε → 0.
Remark 3. The proof of Lemma 2, done in [START_REF] Llibre | Improving the averaging theory for computing periodic solutions of the differential equations[END_REF], is based on the implicit function theorem. So it is easy to check the uniqueness of its predicted branch of zeros

ζ ν ε such that ζ ν ε → ζ ν * when ε → 0. Moreover, it is also a direct consequence of the implicit function theorem that if f 1 (ν) = 0 for every ν ∈ V then g(•, ε)
does not admit branches of zeros converging to V when ε goes to 0.

Proof of Theorem A. For (x, y) ∈ R × R n-2 such that (x, y, 0) ∈ U we define the map δ(x, y; ε) = δ 1 (x, y; ε), δ 2 (x, y; ε), . . . , δ n-1 (x, y; ε) as

δ i (x, y; ε) = ϕ + i (t + (x, y; ε), x, y, 0; ε) -ϕ - i (t -(x, y; ε), x, y, 0; ε), i = 1, 2, . . . , n -1. So δ i (x, y; ε) = δ 0 i (x, y) + εδ 1 i (x, y) + O(ε 2 )
, where

δ 0 i (x, y) = ψ + 0,i (τ + 0 (x, y), x, y, 0) -ψ - 0,i (τ - 0 (x, y), x, y, 0), δ 1 i (x, y) = X + 0,i (ψ + 0 (τ + 0 (x, y), x, y, 0))τ + 1 (x, y) + ψ + 1,i (τ + 0 (x, y), x, y, 0) -X - 0,i (ψ - 0 (τ - 0 (x, y), x, y, 0))τ - 1 (x, y) -ψ - 1,i (τ - 0 (x, y), x, y, 0). ( 8
)
Here O is one of the Landau's symbols, that is g(ε) = O(ε ) for some positive integer if there exists constants

ε 1 > 0 and M > 0 such that |g(ε)| ≤ M|ε | for -ε 1 < ε < ε 1 .
Clearly δ(x ε , y ε , ε) = 0 for some (x ε , y ε , 0) ∈ U if and only if the solution of system (2) passing through (x ε , y ε , 0) is periodic.

In the sequel, for the purpose of proving Theorem A, we identify the elements of Lemma 2.

Take d = n -1, k = n -2, D = {(a, b) ∈ R n-2 × R : (b, a, 0) ∈ U }, V = {ν ∈ R n-2 : r 0 < |ν -ν * | < r 1 }, β = 0, ζ ν = (ν, 0), g 0 (a, b) = δ 0 (b, a), δ 0 1 (b, a) , g 1 (a, b) = δ 1 (b, a), δ 1 1 (b, a) , with δ j = (δ j 2 , δ j 3 , . . . , δ j n-1 ) for j = 0, 1, and Z = {ζ ν = (ν, β(ν)) = (ν, 0) : ν ∈ V } ⊂ D. From (8), g 1 (ζ ν ) = (ν), λ 1 (ν) . Now for each ν ∈ V the solution of system (2), for ε = 0, passing through (0, ν, 0) ∈ A is periodic. Therefore δ 0 1 (0, ν), δ 0 (0, ν) = δ(0, ν, 0) = (0, 0) and g 0 (ζ ν ) = g 0 (ν, 0) = δ 0 (0, ν), δ 0 1 (0, ν) = (0, 0). Moreover Dg 0 (ν, 0) = ⎛ ⎜ ⎜ ⎜ ⎝ ∂δ 0 ∂y (0, ν) ∂δ 0 ∂x (0, ν) ∂δ 0 1 ∂y (0, ν) ∂δ 0 1 ∂x (0, ν) ⎞ ⎟ ⎟ ⎟ ⎠ . So ν = ∂δ 0 1 ∂x (0, ν) = ω 1 (ν) and ν = ∂δ 0 ∂x (0, ν) = (ν). Indeed ∂δ 0 i ∂x (0, ν) = X + 0,i ψ + 0 (σ + 0 (ν), 0, ν, 0) τ + 0 ∂x (0, ν) -X - 0,i ψ - 0 (σ + 0 (ν), 0, ν, 0) τ - 0 ∂x (0, ν) + ∂ψ + 0,i ∂x (σ + 0 (ν), 0, ν, 0) - ∂ψ - 0,i ∂x (σ - 0 (ν), 0, ν, 0),
and computing implicitly the derivative in the variable ε of the equality ϕ ± n (t ± (x, y; ε), x, y, 0; ε) = 0 we get

∂τ ± 0 ∂x (0, ν) = - ∂ψ ± 0,n ∂x (τ ± 0 (0, ν), 0, ν, 0) X ± 0,n (ψ ± 0 (τ ± 0 (0, ν), 0, ν, 0))
.

Therefore

∂δ 0 i ∂x (0, ν) = ∂ψ + 0,i ∂x (σ + 0 (ν), 0, ν, 0) - X + 0,i ψ + 0 (σ + 0 (ν), 0, ν, 0) X + 0,n (ψ + 0 (σ + 0 (ν), 0, ν, 0)) ∂ψ + 0,n ∂x (σ + 0 (ν), 0, ν, 0) - ∂ψ - 0,i ∂x (σ - 0 (ν), 0, ν, 0) + X - 0,i ψ - 0 (σ + 0 (ν), 0, ν, 0, 0) X - 0,n (ψ - 0 (σ - 0 (ν), 0, ν, 0)) ∂ψ - 0,n ∂x (σ - 0 (ν), 0, ν, 0) = ω i (ν),
for i = 1, 2, . . . , n -1. Now we compute the function [START_REF] Colombo | Nondeterministic chaos and the two-fold singularity in piecewise smooth flows[END_REF] as

f 1 (ν) = -ν -1 ν ξ ⊥ g 1 (ζ ν ) + ξg 1 (ζ ν ) = - λ 1 (ν) ω 1 (ν) (ν) + (ν) = M(ν).
Since ω 1 (ν) = 0, the proof of this theorem follows by applying Lemma 2 and Remark 3. 2

Normal forms of piecewise linear systems

Consider the following piecewise linear differential system

Z(x, y, z) = X + (x, y, z), if z > 0 X -(x, y, z), if z < 0, (9)
(x, y, z) ∈ R 3 . The discontinuity set is given by = {z = 0}. Here

X ± = (X ± 1 , X ± 2 , X ± 3 )
, where

X ± j (x, y, z) = a ± j x + b ± j y + c ± j z + d ± j ,
for j ∈ {1, 2, 3}. Note that = f -1 (0), where f (x, y, z) = z. Moreover, we denote

+ = {p ∈ R 3 ; f (p) > 0} and -= {p ∈ R 3 ; f (p) < 0}.
Firstly we assume the existence of a plane transversal to which is invariant through the flow of system [START_REF] Christopher | Limit Cycles of Differential Equations[END_REF]. Without loss of generality, after a linear change of variable, we may assume that = {x = 0}, equivalently b

± 1 = c ± 1 = d ± 1 = 0.
In what follows we assume, without loss of generality, that system (9) has a center at the origin turning the orbits around it in counterclockwise sense, that is b + 3 > 0 and b - 3 > 0. Therefore (0, y, 0) ∈ c , for y = 0 small enough, and d ± 3 = 0. Now (0, y, 0) ∈ c if and only if

(X + f )(0, y, 0)(X -f )(0, y, 0) = b + 3 b - 3 y 2 > 0.
In order to simplify the computations it is convenient to make a continuous piecewise linear change of variables which transforms system (9) into a simple form. This change of variables is a homeomorphism h which keeps invariant the discontinuity set , the plane = {x = 0} and the sets + and -. Furthermore, h restricted to + ∪ -will be a topological equivalence between the vector fields. More precisely, following closely the ideas of [START_REF] Freire | Canonical discontinuous planar piecewise linear systems[END_REF], we have the next result.

Proposition 4. Consider system [START_REF] Christopher | Limit Cycles of Differential Equations[END_REF]

with b ± 1 = c ± 1 = d ± 1 = d ± 3 = 0, b + 3 > 0 and b - 3 > 0. Then the homeomorphism (u, v, w) = h(x, y, z) given by (u, v, w) T = M ± (x, y, z) T if (x, y, z) ∈ ± ∪ ,
where

M ± = ⎛ ⎜ ⎜ ⎜ ⎝ 1 0 0 0 1 c ± 3 b ± 3 0 0 1 b ± 3 ⎞ ⎟ ⎟ ⎟ ⎠ ,
transforms system (9) into the normal form

( u, v, ẇ) T = Y ± (u, v, w) = A ± (u, v, w) T + (0, d ± 2 , 0) T , ( 10
)
where

A ± = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ a ± 1 0 0 a ± 2 b ± 3 +a ± 3 c ± 3 b ± 3 b ± 3 + c ± 3 -(b ± 2 c ± 3 -c ± 2 b ± 3 ) a ± 3 b ± 3 1 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠
.

The discontinuity set and the plane {x = 0} are invariant by h. Moreover the following statements hold.

(i) h transforms the crossing and sliding sets, tangency points and boundary equilibria from the original system (9) into sets and points of same type for system (10);

(ii) h is a topological equivalence between systems (9) and (10) for all their orbits which do not have points in common with e ∪ s ; (iii) h preserves the sets e and s , that is h( e ) and h( s ) are the escaping and sliding regions of system [START_REF] Du | Bifurcation of periodic orbits in a class of planar Filippov systems[END_REF], respectively.

Proof. With straightforward computations system (9) becomes system [START_REF] Du | Bifurcation of periodic orbits in a class of planar Filippov systems[END_REF]. Obviously the set and the plane {x = 0} are invariant by h. The statements (i), (ii) and (iii) follow immediately from the equalities

Y ± f (h(x, y, 0)) = Y ± f (u, v, 0) = 1 b ± 3 X ± f (u, v, 0), because Y + f • Y -f | differs from X + f • X -f | by a positive constant.
Clearly, orbits totally contained in ± are transformed in a homeomorphic way. Furthermore, the topological equivalence is not lost at the crossing set because orbits arriving at the crossing set are continued by the natural concatenation. 2

A point p ∈ where the flow of X + | (resp. X -| ) has a parabolic contact with of the form

X ± f (p) = 0 and (X ± ) 2 f (p) = 0 is called fold point of X + (resp. X -| ). Moreover, when the orbit of X + | (resp. X -| ) through p is locally contained in -(resp. + ), that is (X + ) 2 f (p) < 0 (resp. (X + ) 2 f (p) > 0), the fold point p is called invisible, otherwise, it is called visible.
A point p ∈ ∩ is a center of system ( 9) if there exists a neighborhood V of p such that V ∩ \ {p} is fulfilled by crossing periodic orbit of system [START_REF] Christopher | Limit Cycles of Differential Equations[END_REF]. This means that (X + f )(q)(X -f )(q) > 0 for all q ∈ (V ∩ {x = 0}) ∩ \ {p}, and (X + f )(p) = (X -f )(p) = 0. For system (9) there are three different types of centers in a point p ∈ : focus-focus (FF), focus-parabolic (FP) and parabolic-parabolic (PP). Now we give a briefly description of each one below.

In the focus-focus type both systems X ± | have a singular point at p which has eigenvalues with non-zero imaginary part. In the focus-parabolic (resp. parabolic-focus) type the system X + | (resp. X -| ) has a singular point having eigenvalues with non-zero imaginary part, while the system X -| (resp. X + | ) has an invisible fold. Finally, in the parabolic-parabolic type both systems X ± | have an invisible fold at p.

In order to state the main result of this section, which classifies the centers of system ( 9), we need some preliminary notations:

S ± = (b ± 1 ) 2 + (c ± 1 ) 2 + (d ± 1 ) 2 + (d ± 3 ) 2 , η ± = b ± 2 + c ± 3 , χ ± = -[(b ± 2 + c ± 3 ) 2 -4(b ± 2 c ± 3 -c ± 2 b ± 3 )], and γ ± = η ± /χ ± provided χ ± = 0. We note that S ± = 0 if and only if b ± 1 = c ± 1 = d ± 1 = d ± 3 = 0.
Theorem B. System (9) has a center in the plane = {x = 0} of type:

(i) focus-focus if and only if

S ± = 0, b + 3 b - 3 > 0, d ± 2 = 0, χ ± > 0 and γ + + γ -= 0; (ii) focus-parabolic (resp. parabolic-focus) if and only if S ± = 0, b + 3 b - 3 > 0, d + 2 = 0, d - 2 > 0 (resp. d - 2 = 0, d + 2 < 0), χ + > 0 (resp. χ -> 0) and η ± = 0; (iii) parabolic-parabolic if and only if S ± = 0, b + 3 b - 3 > 0, d + 2 < 0, d - 2 > 0, either η ± = 0 or η + = 0, η -d + 2 = η + d - 2 and (b - 2 c - 3 -c - 2 b - 3 )(d + 2 ) 2 = (b + 2 c + 3 -c + 2 b + 3 )(d - 2 ) 2 .
Remark 5. A point p ∈ is a T-singularity of the vector field associated to system (9), if p is an invisible fold of both vector fields X ± and the fold-curves X ± f = 0 are transversal at p. Therefore in addition to the hypotheses of statement (iii), the parabolic-parabolic center is a T-singularity if and only if a -

3 b + 3 -a + 3 b - 3 = 0.
The proof of Theorem B follows the same ideas of [START_REF] Han | On Hopf bifurcation in non-smooth planar systems[END_REF] for proving its Theorem 3.1. Indeed our proof could be obtained by applying this theorem to the restriction of system ( 9) to . So a major part of it could be omitted here. However, since the conditions on the original parameters that characterize the types of centers are not given explicitly in [START_REF] Han | On Hopf bifurcation in non-smooth planar systems[END_REF], and also for sake of completeness, we prefer to write the whole proof here.

In order to prove Theorem B we shall apply the normal form [START_REF] Du | Bifurcation of periodic orbits in a class of planar Filippov systems[END_REF]. To do that the next lemma is needed. Lemma 6. (See [START_REF] Filippov | Differential Equations with Discontinuous Righthand Side[END_REF].) Consider the equation

dy dx = ax + by + G(x, y), a = 0, ( 11
)
where G is of class C 4 at the origin satisfying

G(x, y) = cx 2 + dxy + ey 2 + f x 3 + gx 2 y + hx 4 + O(x 4 + y 2 ),
for (x, y) near the origin. Let y = Y (x) be the solution of [START_REF] Filippov | Differential Equations with Discontinuous Righthand Side[END_REF] satisfying

Y (-ρ) = Y (σ ) = 0 for -ρ < 0 < σ, aY (x) < 0 for -ρ < x < σ.
Then for ρ > 0 small

σ = ρ + μρ 2 + μ 2 ρ 3 + kρ 4 + O(ρ 5 ),
where

μ = 2 3 b - c a , k = 10 11 μ 3 + μ 5 d + 2 15 L, L = bc 2 a 2 - 2c 3 a 3 -2ae - 2bf a + 5cf a 2 + g - 3h a .
Proof of Theorem B. If system (9) has a center in the plane {x = 0} it satisfies the conditions b

± 1 = c ± 1 = d ± 1 = d ± 3 = 0 and b + 3 b - 3 > 0.
Conveniently we assume the counterclockwise direction for the flow, i.e. b - 3 > 0 and b + 3 > 0. Hence, by Proposition 4 we write system (9) in the form [START_REF] Du | Bifurcation of periodic orbits in a class of planar Filippov systems[END_REF]. Denoting the dependent variables of system (10) by x, y, z, system (10), restricted to the plane {x = 0}, reads

ẏ = (b ± 3 + c ± 3 )y -(b ± 2 c ± 3 -c ± 2 b ± 3 )z + d ± 2 , ż = y, (12)
for z ≷ 0, where (0, y, z) ∈ ± ∪ . Now we start the proof of (i). In this case the upper and lower systems have a singular point at the origin, i.e. d ± 2 = 0, with eigenvalues given by

η ± 2 + i χ ± 2 and η ± 2 -i χ ± 2
, where χ ± > 0.

Let ϕ ± (t, (0, r, 0)) = (ϕ ± 1 (t, (0, r, 0)), ϕ ± 2 (t, (0, r, 0))) be the flow of system [START_REF] Freire | Canonical discontinuous planar piecewise linear systems[END_REF] where

ϕ ± 1 (t, (0, r, 0)) = re η ± 2 t cos χ ± 2 t + γ ± sin χ ± 2 t , ϕ ± 2 (t, (0, r, 0)) = 2r χ ± e η ± 2 t sin χ ± 2 t .
Therefore the half Poincaré maps of upper and lower systems are given by

π ± (r) = ϕ ± 1 2π χ ± , (0, r, 0) = -re γ ± π , where 2π χ ± is such that ϕ ± 2 2π
χ ± , (0, r, 0) = 0. Hence the Poincaré map is

π(r) = π -• π + (r) = re (γ -+γ + )π .
Thus we have a center if γ -+ γ + = 0, and statement (i) is proved. Under the assumptions of statement (ii) the upper system (12) has a singular point at the origin, i.e. d + 2 = 0 and the eigenvalues of Jacobian matrix are

η + 2 ± i χ + 2 with χ + > 0.
Analogous to the previous case the half Poincaré map of the upper system is π + (r) = -re γ + π . On the other hand, the lower system (12) has an invisible fold at the origin, i.e. d - 2 > 0. Consider the differential equation

dz dy = y (b - 3 + c - 3 )y -(b - 2 c - 3 -c - 2 b - 3 )z + d - 2 , ( 13
)
associated to the lower system [START_REF] Freire | Canonical discontinuous planar piecewise linear systems[END_REF]. Since d - 3 > 0 we can expand the right hand side of (13) around y = 0 and obtain the following expression

dz dy = 1 d - 2 y - (b - 3 + c - 3 ) (d - 2 ) 2 y 2 + (b - 2 c - 3 -c - 2 b - 3 ) (d - 2 ) 2 yz + (b - 3 + c - 3 ) (d - 2 ) 2 y 3 - 2(b - 3 + c - 3 )(b - 2 c - 3 -c - 2 b - 3 ) (d - 2 ) 3 y 2 z - (b - 3 + c - 3 ) 3 (d - 2 ) 4 y 4 + O(y 4 + z 2 ). ( 14
)
Let z = z(y) be the solution of equation ( 14) with initial condition z(-r) = 0. Then z(π -(-r)) = 0 and by Lemma 6

π -(-r) = r + μ -r 2 + (μ -) 2 r 3 + k -r 4 + O(r 5 ),
where

μ -= 2 3 b - 2 + c - 3 d - 2 , L -= -2(b - 2 c - 3 -c - 2 b - 3 )(b - 2 + c - 3 ) (d - 2 ) 3 , k -= 2 3 (b - 2 + c - 3 ) d - 2 ⎡ ⎣ 40 99 b - 2 + c - 3 d - 2 2 - 1 5
(b - 2 c - 3 -c - 2 b - 3 ) (d - 2 ) 2 ⎤ ⎦ = μ - ⎡ ⎣ 40 99 b - 2 + c - 3 d - 2 2 - 1 5
(b - 2 c - 3 -c - 2 b - 3 ) (d - 2 ) 2 ⎤ ⎦ . ( 15
)
Thus the Poincaré map is

π(r) = π -• π + (r) = π -(-re γ + π ) = e γ + π r + μ -e 2γ + π r 2 + O(r 3 ).
Hence a necessary condition for having a center at the origin is e γ + π = 1 and μ -e 2γ + π = 0, i.e. γ + = μ -= 0, and so η ± = 0. In this case we get π + (r) = r, and the lower system [START_REF] Freire | Canonical discontinuous planar piecewise linear systems[END_REF] becomes

ẏ = b - 3 (b - 2 + c - 2 )z + d - 2 , ż = y,
which is invariant under the change of variables (y, z, t) → (-y, z, -t). Therefore π -(-r) = r and so π(r) = r, i.e. the origin is a center. Now the case parabolic-focus can be reduced to the case focus-parabolic trough the change of variables (x, y, t) → (-x, -y, t). Hence statement (ii) is proved.

Finally we prove statement (iii). Then the lower and the upper systems ( 12) have an invisible fold at the origin, i.e. d + 2 < 0 and d - 2 > 0. As in the previous statement applying Lemma 6, we obtain

(π + ) -1 (-r) = r + μ + r 2 + (μ + ) 2 r 3 + k + r 4 + O(r 5 ),
where μ + , L + and k + are the same as in [START_REF] Jacquemard | Piecewise smooth reversible dynamical systems at a two-fold singularity[END_REF] just changing the sign "-" by "+". Hence for r > 0 we obtain

π -(-r) -(π + ) -1 (-r) = (μ --μ + )r 2 + [(μ -) 2 -(μ + ) 2 ]r 3 + (k --k + )r 4 + O(r 5 ).
Then, necessary conditions for having a center at the origin are μ -μ + = 0 and k -k + = 0. Observe that the first condition implies that d + 2 η -= d - 2 η + . Now if η + = 0 we have η -= 0 and k + = k -= 0. For η + = 0 the condition k

+ -k -= 0 implies that (d + 2 ) 2 (b - 2 c - 3 -c - 2 b - 3 ) = (d - 2 ) 2 (b + 2 c + 3 -c + 2 b + 3
). In the first case, i.e. η + = 0, system [START_REF] Freire | Canonical discontinuous planar piecewise linear systems[END_REF] becomes

ẏ = -(b ± 2 c ± 3 -c ± 2 b ± 3 )z + d ± 2 , ż = y,
which is invariant under the change of variables (y, z, t) → (-y, z, -t). Therefore (π + ) -1 (-r) = π -(-r) = r and so π --(π + ) -1 ≡ 0, i.e. the origin is a center.

In the second case the upper and lower systems [START_REF] Freire | Canonical discontinuous planar piecewise linear systems[END_REF] become the systems

ẏ = (b + 3 + c + 3 )y -(b + 2 c + 3 -c + 2 b + 3 )z + d + 2 , ż = y, (16)
when (0, y, z) ∈ + ∪ , and

ẏ = d - 2 d + 2 (b + 3 + c + 3 )y - d - 2 d + 2 2 (b + 2 c + 3 -c + 2 b + 3 )z + d - 2 , ż = y, (17)
when (0, y, z) ∈ -∪ , respectively. Note that by the change of variables

(y, z, t) → y, d - 2 d + 2 z, d - 2 d + 2 t , (18)
system (16) becomes system [START_REF] Liu | Poincaré bifurcation of a three-dimensional system[END_REF]. Hence if ϕ -(t, (0, -r, 0)) = ϕ - 1 (t, (0, -r, 0)), ϕ - 2 (t, (0, -r, 0))

is the solution of system (17) through the point (0, -r, 0), then r), (0, -r, 0)) = 0. Thus, using the inverse of the change of variables [START_REF] Llibre | Averaging theory for discontinuous piecewise differential systems[END_REF], it follows that

π -(-r) = ϕ - 1 (t -(-r), (0, -r, 0)), where t -(-r) is such that ϕ - 2 (t -(-
ϕ + (t, (0, -r, 0)) = (ϕ + 1 (t, (0, -r, 0)), ϕ + 2 (t, (0, -r, 0))) = ϕ - 1 d + 2 d - 2 t, (0, -r, 0) , d + 2 d - 2 ϕ - 2 d + 2 d - 2
t, (0, -r, 0) , is the solution of system (16) through the point (0, -r, 0). Therefore r) and (π + ) -1 (-r) = π -(-r), i.e. the origin is a center. 2

(π + ) -1 (-r) = ϕ + 1 (-t + (-r), (0, -r, 0)), where t + (-r) is such that ϕ + 2 (-t + (-r), (0, -r, 0)) = 0. Then d + 2 t + (-r) = -d - 2 t -(-

Perturbation of the normal forms

In this section we use the theory developed in section 2 to study the persistence of crossing periodic solutions of the normal forms given by Theorem B, when it is perturbed inside the class of piecewise linear differential systems with two zones separated by the plane = {z = 0}.

Parabolic-parabolic center

We consider the unperturbed systems X a 0 (x, y, z) = x, -zsign(z), x + y and X b 0 (x, y, z) = x, zsign(z), x + y .

From statement (iii) of Theorem B it is easy to see that the origin is a parabolic-parabolic center of the unperturbed systems ( ẋ, ẏ, ż) T = X a 0 (x, y, z) and ( ẋ, ẏ, ż) T = X b 0 (x, y, z). Indeed for both systems

S ± = 0, b + 3 b - 3 = 1 > 0, d + 2 = -1 < 0, d - 2 =
1 > 0, and η ± = 0. The singular points of the upper and lower vector fields associated to X a 0 are invisible centers and so the period annulus A is all the plane {x = 0}, see Fig. 2. Now, the singular points of the upper and lower vector fields associated to X b 0 are visible saddles and so the period annulus A is bounded by the separatrices of the saddles, see Fig. 3.

The theory developed in the previous sections allows us to study the limit cycles of the following perturbed systems X a (x, y, z; ε) :

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ẋ = x + ε α + 0 + α + 2 y + α + 3 z ẏ = -z -1 + εβ + 2 y ż = x + y + ε κ + 0 + κ + 3 z if z > 0, ẋ = x + ε α - 0 + α - 2 y + α - 3 z ẏ = -z + 1 + ε β - 0 + β - 2 y + β - 3 z ż = x + y + ε κ - 0 + κ - 2 y + κ - 3 z if z < 0. ( 19
)
So we have the following result.

+ α + 0 -2α + 0 sech π √ 3 -3(β - 0 + β + 0 + 2κ - 0 ) and B = 2(α - 2 + α - 3 + 3(α + 2 + α + 3 )) √ 3 e π √ 3 -1 + tanh π √ 3 -1 (α - 2 + α - 3 + 3(α + 2 + α + 3 )) √ 3 - 2e - √ 3π (α - 2 + α - 3 ) √ 3 + 2πα - 2 -πα - 3 + 2 √ 3α + 2 -2πα + 2 + 2 √ 3α + 3 -πα + 3 + 3πβ - 3 -2πβ + 2 -πβ + 3 -3πκ - 2 + 3πκ - 3 + πκ + 2 -πκ + 3 .
Moreover if A = 0 and B = 0 then there are no crossing periodic solutions bifurcating from periodic orbits contained in {x = 0}.

It is worth to say that even assuming B = 0 and |ε| = 0 small enough this proposition does not provide the global uniqueness of crossing periodic solutions of system [START_REF] Llibre | Improving the averaging theory for computing periodic solutions of the differential equations[END_REF] close to the plane {x = 0}. Indeed some crossing periodic solutions may bifurcate purely from infinite (see [START_REF] Gouveia | On limit cycles bifurcating from the infinity in discontinuous piecewise linear differential systems[END_REF]).

Proof of Proposition 9. The proof of the first part of this proposition is analogous to the proof of Proposition 7, noting that the solutions of ( ẋ, ẏ, ż) T = X c 0 (x, y, z) are ψ + 0 (t, 0, ν, 0) = 0, νe Clearly if A = 0 and B = 0, then M(ν) = C 0 = 0 for all ν ∈ V. It concludes this proof. 2

Fig. 2 .

 2 Fig. 2. Phase Portrait of unperturbed vector field X a 0 .

Proposition 9 .

 9 If B = 0 then, for |ε| = 0 small enough, there is a unique crossing periodic solution φ(t, ε) of system[START_REF] Llibre | Improving the averaging theory for computing periodic solutions of the differential equations[END_REF] bifurcating from periodic orbits contained in {x = 0}, such that φ(0, ε)

  in this case, the Melnikov-like function (5) is linear, i.e.M(ν)= C 0 + C 1 ν, where

Acknowledgments

We thank to the referee for the helpful comments and suggestions. The first and fourth author are partially supported by a FAPESP grant 2013/34541-0. The first and fourth authors are supported by a CAPES PROCAD grant 88881.068462/2014-01. The second author is partially supported by a MINECO grant MTM2013-40998-P, an AGAUR grant number 2014SGR-568, the grants FP7-PEOPLE-2012-IRSES 318999 and 316338, and by a CAPES CSF-PVE grant 88881.030454/2013-01. The third author is supported by a FAPESP grant 2015/02517-6.

Proposition 7. There exist parameters α ± i , β - i , κ - i , i = 1, 2, 3, β + 2 , κ + 0 , and κ + 3 such that system [START_REF] Llibre | On the birth of limit cycles for non-smooth dynamical systems[END_REF] admits, for |ε| = 0 small enough, at least 8 limit cycles converging, when ε goes to 0, to some of the periodic orbits contained in {x = 0}. Proposition 8. There exist parameters α ± i , β - i , κ - i , i = 1, 2, 3, β + 2 , κ + 0 , and κ + 3 such that system (20) admits, for |ε| = 0 small enough, at least 8 limit cycles converging, when ε goes to 0, to some of the periodic orbits contained in {x = 0}.

Proof of Proposition 7. In order to prove this proposition we have to identify the elements of Theorem A. Computing the solutions of ( ẋ, ẏ, ż) T = X a 0 (x, y, z) we obtain

Expanding the solutions of system [START_REF] Llibre | On the birth of limit cycles for non-smooth dynamical systems[END_REF] in power series of ε we get

and

Analogously we compute σ ± 0 (t, 0, ν, 0) = ±2 arctan ν,

and

Using the formulae (6) we conclude

e 2 arctan ν 2

,

where

Now expanding the function M(ν) (see [START_REF] Chicone | Lyapunov-Schmidt reduction and Melnikov integrals for bifurcation of periodic solutions in coupled oscillators[END_REF]) in power series of ν we get

Here the parameters C i , i = 0, . . . 8, are linear combinations of the parameters

2 , κ + 0 , and κ + 3 . Moreover they are linearly independent. So we can choose them in order to obtain 8 positive simple zeros of the equation M(ν) = 0. The proof of this theorem follows by applying Theorem A. 2

Proof of Proposition 8. The proof of this proposition is analogous to the proof of Proposition 7, noting that the solutions of ( ẋ, ẏ, ż) T = X b 0 (x, y, z) are ψ + 0 (t, 0, ν, 0) = 0, ν cosh tsinh t, ν sinh tcosh t + 1 , ψ - 0 (t, 0, ν, 0) = 0, ν cosh t + sinh t, ν sinh t + cosh t -1 , and

Here again, expanding the function M(ν) (see [START_REF] Chicone | Lyapunov-Schmidt reduction and Melnikov integrals for bifurcation of periodic solutions in coupled oscillators[END_REF]) in power series of ν we get

where the parameters C i , i = 0, . . . 8, are independent linear combinations of the parameters α

Focus-focus center

We consider the unperturbed systems X c 0 (x, y, z) = x, -z + sign(z)y, x + y . From statement (ii) of Theorem B it is easy to see that the origin is a focus-focus center of the unperturbed system ( ẋ, ẏ, ż) T = X c 0 (x, y, z) in such way that their periodic orbits fulfill all the plane {x = 0} (see Fig. 4). Indeed for both systems

3, η ± = ±1 and then γ 1 + γ 2 = 0.

We study the limit cycles of the following perturbed system X c (x, y, z; ε) :