Antonio Papangelo
email: antonio.papangelo@poliba.it

Aurelien Grolet

Loic Salles

23 Norbert Hoffmann

Michele Ciavarella

Snaking bifurcations in a self-excited oscillator chain with cyclic symmetry

Keywords: snaking bifurcation, subcritical Hopf bifurcation, localised vibration, self-excitation, bistability, nonlinear dynamics

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Spatially localised states of dynamical systems have been studied in a large number of different fields in the sciences and in engineering. While for linear systems Anderson localisation was the key to quite a satisfactory understanding, in nonlinear dynamical systems the quest to understand localisation seems far from settled. For a long time progress seems to have been largely confined to conservative nonlinear systems, where solitons and breathers made their appearance. Only later, dissipative systems have come into focus, with first work based on tracing solitons into the driven and dissipative regime, introducing dissipative solitons. In parallel to the study of solitary states in conservative and dissipative systems, another breakthrough to the understanding of spatial localisation in dissipative localisation was accomplished in the study of subcritical bifurcations in pattern-forming systems, where the concept of branching has emerged and is a well-established field of study today.

Branching is today well known in a number of disciplines, amongst others in optics [START_REF] Champneys | Homoclinic orbits in reversible systems and their applications in mechanics, fluids and optics[END_REF], granular matter [START_REF] Umbanhowar | Localized excitations in a vertically vibrated granular layer[END_REF], structural mechanics [START_REF] Champneys | A multiplicity of localized buckling modes for twisted rod equations[END_REF][4] [START_REF] Lord | Computation of localized post buckling in long axially compressed cylindrical shells[END_REF][6] [START_REF] Grolet | Vibration analysis of a nonlinear system with cyclic symmetry[END_REF] [START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF], and mostly in fluid dynamics [START_REF] Champneys | Homoclinic orbits in reversible systems and their applications in mechanics, fluids and optics[END_REF] [START_REF] Niemela | Localized traveling-wave states in binary-fluid convection[END_REF][10] [START_REF] Batiste | Spatially localized binary-fluid convection[END_REF][12] [START_REF] Thual | Localized structures generated by subcritical instabilities[END_REF], and magnetohydrodynamics [START_REF] Blanchflower | Magnetohydrodynamic convectons[END_REF] [START_REF] Dawes | Localized convection cells in the presence of a vertical magnetic field[END_REF]. The first studies into the topic have probably emerged in the field of binary-fluid convection, where spatially localised convection rolls have been observed in water-ethanol mixtures [START_REF] Niemela | Localized traveling-wave states in binary-fluid convection[END_REF] or helium [START_REF] Woodcraft | Visualisation of convective flow patterns in liquid helium[END_REF]. There localised convection domains of arbitrary length are found to be stable, being surrounded by the conductive state.

In terms of bifurcation diagrams, the localised states have shown to be arranged in a unique and fascinating way, giving birth to what is called a snaking structure [START_REF] Batiste | Spatially localized binary-fluid convection[END_REF] [START_REF] Beaume | Homoclinic snaking of localized states in doubly diffusive convection[END_REF][15] [START_REF] Hunt | Cellular buckling in long structures[END_REF]. The snaking structure has e.g. been studied a lot in the one dimensional [START_REF] Burke | Homoclinic snaking: structure and stability[END_REF] and the two dimensional [START_REF] Avitabile | To snake or not to snake in the planar Swift-Hohenberg equation[END_REF] Swift-Hohenberg equation, which is a convenient and generic model system to study fundamental properties of the arising dynamics.

A typical snaking bifurcation diagram involves two snaking solution branches, intertwined into each other. Figure 1 gives an example from convection [START_REF] Beaume | Homoclinic snaking of localized states in doubly diffusive convection[END_REF] to illustrate the phenomenon. For the bifurcation diagram (left panel), the average kinetic energy "E" of the flow is plotted versus the Rayleigh number "Ra" and two intertwined branches appear. In the right panel, nineteen solutions are shown, which correspond to the numbers positioned close to the snaking structure, in which spatially localised convective rolls can be identified. Notice that the higher the energy of the solution, the larger the number of convection cells. Often the two snaking branches are also interconnected through a number of unstable branches, and a ladder like pattern emerges [START_REF] Avitabile | To snake or not to snake in the planar Swift-Hohenberg equation[END_REF].

Although snaking bifurcations are now generally known and studied in many fields of dynamical systems, it seems that there is hardly any study into the phenomenon in the context of structural vibrations in engineering. In many respects this is quite surprising, since non-linear oscillators with subcritical Hopf bifurcations, often coupled to neighbouring oscillators of the same type into chains or arrays, are actually very common models for a number of systems from engineering vibrations. And also the appearance of bi-or multi-stability, which is obviously at the core of the phenomenon [START_REF] Thual | Localized structures generated by subcritical instabilities[END_REF] [START_REF] Burke | Homoclinic snaking: structure and stability[END_REF] [START_REF] Avitabile | To snake or not to snake in the planar Swift-Hohenberg equation[END_REF], is well established in many of these engineering systems. Moreover, the emergence of spatially localised vibration states in structural dynamics is also a well known observational fact: e.g. in turbo-machinery, there is the so-called effect of 'mis-tuning in rotors' [START_REF] Whitehead | Effect of mistuning on the vibration of turbo-machine blades induced by wakes[END_REF] [START_REF] Ewins | The effects of detuning upon the forced vibrations of bladed disks[END_REF]. Traditionally, the origin of the localisation is thought to have its root in slight system inhomogeneities, leading to linear localisation in the sense of Anderson. From testing, strong localisation is confirmed, but proper validation of the theory has up to now not been accomplished in the linear framework. In a sense it is tempting to hypothesise that one of the key reasons behind might be the non-linearity involved, which definitely becomes substantial for the large local vibration amplitudes observed. To the best of our knowledge, in model systems for turbo-machinery dynamics, snaking behaviour has never been investigated. Also systems from fluid-structure-interaction, may show weak non-linearity, Hopf bifurcation, and bi-stability, like models for aerofoil flap dynamics [START_REF] Liu | Bifurcation analysis of airfoils in incompressible flow[END_REF][21] [START_REF] Chen | Supercritical as well as subcritical Hopf bifurcation in nonlinear flutter systems[END_REF] [START_REF] Pereira | Insights on aeroelastic bifurcation phenomena in airfoils with structural nonlinearities[END_REF]. Similarly in friction induced vibrations the emergence of snaking could be well expected, with all the necessary ingredients like flutter instability and bi-stability already known to exist, cf. e.g. [START_REF] Weiss | Friction induced dynamics of ball joints: Instability and post bifurcation behavior[END_REF], [START_REF] Gräbner | Nonlinearities in Friction Brake NVH-Experimental and Numerical Studies[END_REF].

We will thus consider a model system as simple as we can think of, but derived from models actually in use in the turbo-machinery community and the field of fluid-structure interaction and friction-induced vibration. We choose a chain of (weakly non-linear) oscillators coupled into a linear oscillator chain. For simplicity we close the chain into a cyclically symmetric ring, which moreover has the advantage of bringing it even closer to models used widely in turbo-machinery for rotors with a small but finite number of blades attached. To obtain or model instability and spatially local bifurcation, i.e. the bi-stability of the individual oscillator, we introduce non-linear damping terms, i.e. non-linear terms depending on velocities. This approach is heuristic and rather for simplicity at the present stage of understanding, but can be thought of bringing into our purely structural model the corresponding non-linear forcing and dissipation terms from surrounding flow, or an involved friction interface.

With a velocity dependent force arising from a fifth order polynomial representation, our system results in individual uncoupled oscillators showing subcritical Hopf bifurcations and bistability, while the oscillator chain, i.e. the coupled oscillators, yield snaking bifurcations that we determine by solving algebraic equations or time-integration. Interestingly, the results do turn out different to many of the hitherto reported snaking bifurcation patterns. The snaking observed in our study seems to show more than a single snake-and-ladders pattern, and much of the pattern as a whole seems to have disintegrated into isolated branches, now usually called isolas. Reviewing where these differences might come from, shows that the key features where our system is different, can easily be identified. First of all, our system is discrete. Individual self-excited oscillators are coupled, as e.g. in the work of Yulin and Champneys in [START_REF] Yulin | Discrete snaking: multiple cavity solitons in saturable media[END_REF], where a one-dimensional periodic array of optical cavities pumped by coherent light were studied. As in our results, the effect of discreteness was studied, showing that the pinning region (the parameter interval where the snaking occurs) gets progressively narrower as the continuum limit is approached. The second characteristic of our system that is slightly unusual is its finite size due to the cyclic symmetry, which does not allow arbitrary wave-numbers or wavelengths to appear, and so does put a constraint onto the system. Formally similar, Taylor and Dawes in [START_REF] Taylor | Snaking and isolas of localised states in bistable discrete lattices[END_REF], studied snaking and localised states in spatially discrete problems for modified periodic or phase shifted boundary conditions. Also in their case, isola rather than continuous snakes have been observed. Figure 1. An example for snaking in a convection system. Adapted from [START_REF] Beaume | Homoclinic snaking of localized states in doubly diffusive convection[END_REF]. Left: typical snaking pattern with two intertwined solution branches in the bifurcation diagram with kinetic energy "E" of the fluid plotted versus the Rayleigh number "Ra" . Right: The stream-functions for solutions as marked in the snaking pattern. The graphs show the spatially localised convection patches.

The mechanical system

We consider a cyclic system of   non-linear oscillators, see Figure 2, which are coupled via a weak linear spring of stiffness  ∆ . Each oscillator has mass  and is linked to the ground via a linear spring  and a non-linear damper which introduces a velocity proportional force of the form

  = - 1 •  +  3 •  3 - 5 •  5  (1)
Here  denotes the displacement of the individual oscillator,

•  the velocity, and we introduced the coefficients  1   3   5 to parametrise the velocitydependent force. The evolution equations for the individual oscillators read

 ••   +  1 •   - 3 •  3  +  5 •  5  +   - ∆ ( +1 +  -1 -2  ) = 0 (2)
where the stiffness  ∆ couples the -th mass with the neighbouring ones. We introduce the quantities

 0 = p   ∆ =  ∆    =   2 √    =  0 , divide the equation (2) using the group  2 0  0 , with  0 a reference displacement, obtaining  1 •• e   +  2 • e   - 3 • e  3  +  4 • e  5  +  5 e   - 6 (e  +1 + e  -1 -2e   ) = 0 (3)
where

 1 = 1  2 = 2 1   3 = 2 3  2 0  2 0   4 = 2 5  4 0  4 0   5 = 1  6 =  ∆  (4)
and the e ¤ superposed indicates that the new displacements are dimensionless, e  ( ) = ( ) 0  In (3) we defined a dimensionless time  =  0 , which allows to replace   with  0   . Notice that we choose  3   5  0, thus the third degree term of the velocity-dependent force introduces a destabilizing force into the system, while the fifth degree term tends to stabilize it. Figure 3(a) lists the arbitrary but characteristic parameters that will be used in the next sections. We will choose  1 , i.e. the linear damping coefficient, as our primary control parameter in a range from -04 to +06. In this range the velocity dependent force changes its shape in a way such that for low  1 values a negative damping is introduced, which is often used in the literature to model self-excited vibrations, such as in fluid-or friction-induced flutter [START_REF] Billah | Resonance, Tacoma Narrows bridge failure, and undergraduate physics textbooks[END_REF] or squeal [START_REF] Rudd | Wheel/rail noise-part II: wheel squeal[END_REF].

Numerical algorithm

Harmonic Balance Method (HBM)

The problem to solve is composed of   second order differential equations. In this study periodic vibrations will be considered only. We therefore apply the Harmonic Balance Method (HBM) as an efficient numerical technique to obtain an approximation to the steady-state solution of the system. In the following we briefly recall the main steps of the HBM, further details can e.g. be found in [START_REF] Nayfeh | Nonlinear oscillations[END_REF]. Consider a time dependent signal  () and express it in a Fourier series

() =  0 2 +   X =1 (  cos () +   sin ())  ( 5
)
where  is the fundamental frequency,   is the number of harmonics considered and

  = 2  Z 2 -2 () cos ()  (6)
  = 2  Z 2 -2 () sin ()  (7)
A set of non-linear algebraic equations is obtained by substituting the dimensionless form of eq. ( 5) into (3) and projecting the equilibrium equations on the Fourier basis 1 sin ( )  cos ( ) for  = 1     The projection gives back a system of algebraic equations in which the unknowns are the Fourier coefficients      of the   degrees of freedom (dof) considered. In our particular case the dynamical system is autonomous (there is no forcing term) and therefore the angular frequency of the solution  is also treated as an unknown.

Numerical solution and continuation

The number of algebraic equations to be solved is (2    + 1). If only one harmonic is considered (  = 1) the system of equations for the -th mass is

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 5 8  4  5 ( 2  +  2  ) 2   -3 4  3  3 ( 2  +  2  )   +  ( 2   - 1   ) + +  ( 5 + 2 6 ) - 6 ( -1 +  +1 ) = 0 5 8  4  5 ( 2  +  2  ) 2   -3 4  3  3 ( 2  +  2  )   +  ( 2   +  1   ) + -  ( 5 + 2 6 ) +  6 ( -1 +  +1 ) = 0 (8)
where   (  ) is the first harmonic cosine (sine) coefficient. In subsection 4.1 it will be shown that in most parameter ranges even with just one harmonic the approximate solution is very close to the exact one, therefore we will mostly use one harmonic only, which leads to a system of 2  polynomial equations of the fifth order. The other equation comes from the projection of the equilibrium equation on 1 which would allow for evaluation of the mean displacement over a period. As the non-linearities are odd and there is no constant term in the equations, the mean displacement always vanishes. The system of 2  polynomial equations has been solved using a Newton-Raphson scheme implemented in the MATLAB R °function fsolve. For a fixed value of  1 , a set of initial conditions has been provided using steady state solutions obtained from a time integration algorithm. Localised as well as non localised solutions were used as a starting point for the Newton-Raphson algorithm. To solve the system we set for one mass  1 = 0, which allows us to solve for the fundamental frequency of the solution too. A continuation algorithm, also implemented in MATLAB R °, continued the solution using a pseudo arc-length scheme [START_REF] Nayfeh | Applied nonlinear dynamics: analytical, computational and experimental methods[END_REF] which allowed to follow the solution trajectory even when turning points were encountered [START_REF] Nayfeh | Applied nonlinear dynamics: analytical, computational and experimental methods[END_REF].

Simulation Results

Single oscillator dynamics

In this subsection we concentrate our attention on the behaviour of a single oscillator when it is isolated from the rest of the chain. The following parameters will be used:

 0 = 2  0 = 1  3 = 03  5 = 01  ∆ = 0 (9)
In Figure 3(b) we plot the maximum potential energy of the mass in dimensionless form: e

 max = 1 2 e  2 max
, where e  max is the vibration amplitude at the steady-state condition. The results have been obtained using three different approaches: circles represent the result of the Time Integration (TI), the solid line shows the result obtained using the HBM where  () is approximated with only the first harmonic, while the dashed line is the result obtained using the HBM with two harmonics, the first and the third. The choice of the first and third harmonics takes into account the fact that the non-linearities introduced are odd and symmetry breaking bifurcations are not considered in this work. We start focusing on the TI results (circles). When  1 & 02 only one solution exists in which for every initial condition the vibration is damped and vanishes. Decreasing  1 , in the range 0   1 . 02 another solution appears with e  max  0 which corresponds to a stable limit cycle of finite amplitude. The single oscillator experiences a subcritical Hopf bifurcation, and in the interval 0   1 . 02 two different stable solutions exist. Notice that using TI, only the stable solutions can be found (Figure 3(b) circles), while HBM also allows to obtain the unstable ones and to follow them. If  1  0 the steady state is not stable and the vibration amplitude grows up to the upper branch with e  max  0. In this range the only stable solution is a limit cycle for which the amplitude grows when  1 is further reduced. Figure 3(b) also shows the results obtained using HBM with one harmonic (solid line) and with two harmonics (dashed line). First we note that in the interval 0   1 . 02 another solution appears which represents an unstable limit cycle. Differences in the two curves obtained with HBM appear for the upper branch when  1 decreases, as higher vibration amplitude implies higher contribution of the non-linearities. The solution which includes two odd harmonics approximates very well the exact solution obtained with TI scheme. Nevertheless if we just focus on the interval 0   1 . 02 we can see that even the single harmonic approximation is never too far from the exact solution 1 . In the following sections the oscillator chain dynamics is studied in the bistability zone (0   1 . 02), and all the results presented will be obtained using HBM with one harmonic.

 1 = [-04 -025 -01  06]  3 = 03  5 = 01  0 = 2  0 = 1 (b)
Maximum potential energy of one single oscillator plotted against the bifurcation parameter  1  Blue dots: results of the time integration at the steady state. Solid line: HB solutions with a single harmonic. Dashed line: HB solutions with two odd harmonics.

Oscillator chain dynamics: linear system

Before studying the dynamical behaviour of a cyclically symmetric chain of   = 12 non-linear oscillators we analyse the response of the underlying undamped linear system with a dimensionless coupling stiffness  ∆ = 001. Due to the cyclic symmetry this system presents mostly pairwise degenerate eigenfrequencies that can be computed as [START_REF] Samaranayake | Subharmonic oscillations in harmonically excited mechanical systems with cyclic symmetry[END_REF]

  = 1 + 2 ∆ µ 1 -cos µ 2   ¶ ¶  ( 10
)
where N and 0 ≤  ≤   2 for even   or 0 ≤  ≤ (  -1) 2 for odd    Due to the weak coupling among the oscillators the natural frequencies will lie on a narrow band. A possible set of orthogonal normal modes is [START_REF] Samaranayake | Subharmonic oscillations in harmonically excited mechanical systems with cyclic symmetry[END_REF]

 0 = [1 1  1]   (11)
   = [cos (  )  cos (2  )   cos (    )]   (12)
   = [sin (  )  sin (2  )   sin (    )]   (13)
The following table lists the natural frequencies [rads] and Figure 4 shows the normal mode shapes of the undamped linear system.

 0 = 1  1 = 1001  2 = 1005  3 = 1010  4 = 1015  5 = 1019  6 = 1020

Oscillator chain dynamics: non-linear system

We now study a cyclically symmetric chain of   = 12 non-linear oscillators in the bistability zone. The parameters used are the same we used for the single oscillators above, except for  ∆ = 001, which introduces a small coupling between the oscillators. The results obtained from time integration were used to derive initial conditions for the continuation algorithm. Figure 5 shows, in the left panel, the sum of the maximum potential energy of each mass e

 max = 1 2 P  =1 e  2 max = 1 2 P  =1 ³ e  2  + e  2
 ´plotted against the linear damping coefficient  1  Many solutions appear to be entangled, making it almost impossible to distinguish one from the other. Looking more closely at the overall structure created by the superposed solutions we can observe trajectories similar to snake and ladder branches [START_REF] Batiste | Spatially localized binary-fluid convection[END_REF][12][6] [START_REF] Dawes | Localized convection cells in the presence of a vertical magnetic field[END_REF] and twelve 'steps' (corresponding to the number of oscillators) can be identified. Each step is labelled with a red circle and the corresponding energy distribution is plotted in the twelve bar plots on the right-hand side of Figure 5. In each bar plot the mass number has been reported on the x-axis, while the bar height is computed as

e   = e  2 max max ¡ e  2 max ¢ =112  (14)
such that the energy of each oscillator is normalized with respect to the one which is vibrating with the largest amplitude. The bar plots 1 -12 in Figure 5 shed some light on the system dynamics: each step can be easily related to the number of masses on which the vibration is localised. From the bottom to the top, at the first step one mass is moving while the others are more or less motionless, at the second step two masses are moving and so on up to the twelfth. This resembles the usual snaking behaviour, where, for example in fluid dynamics, the steps can be related to the increasing number of convection rolls. Especially in the bottom part of our snaking (Figure 5, left panel) it is possible to see very close similarities with the classical snaking picture. On the other hand, when increasing the energy, and thus the influence of the non-linearities, the picture gets more and more distorted and only seems to be bounded by the very last branch in which all the masses are moving.

Consequently one might conjecture that the relative phases between the individual oscillators play a special role in our system. Figure 6 thus shows the vibration mode shape, indicating the relative phasing between the oscillators, in 12 bar plots numbered from 1 to 12, each-one related to one of the red circles in the bifurcation diagram (Figure 5, left panel). It seems that there is an additional sub-structure due to relative in-or out-of-phase behaviour within our localised vibration zones. Sometimes an oscillator is out of phase with both of its neighbours, sometimes two adjacent oscillators go in phase, but out of phase to their respective other neighbours. At the present state of our study, we have not yet succeeded to obtain a deeper understanding of this result, but think it is at the very heart of the multiplicity of solutions observed. Further studies will need to follow.

At this point it also seems useful to have another look at the numerical accuracy of our computational approach. Figure 7 shows in each column results obtained from TI for four solutions using  1 = 01 and different initial conditions. In particular from the left to the right the vibration is localised on one mass, five masses, eight masses and on all the twelve masses. In the first row the displacement time history is shown at the steady-state, in the second row the solution is shown in the phase space and in the third row the dimensionless maximum potential energy of each mass is given (TI results, dark blue). To assess the effect of the first harmonic truncation on the energy distribution among the different oscillators, HBM results are computed for the same parameter sets and using the TI solutions as starting point. The maximum potential energy, shown in Figure 7 (HBM results, pale yellow bars), proves that the energy is distributed exactly in the same manner and that the low-order approximation causes only a very small reduction in the potential energy of each mass (due to the loss of solution details) which seems acceptable for our purposes. Notice also that the shape of the solution with all the masses vibrating corresponds very well to the last mode of the linearised system, in which the motion is out-of-phase.  vs e ), in the third row the dimensionless maximum potential energy of each mass in steady state condition is shown (dark blue). In the third row, the energy distribution obtained with single harmonic HBM is plotted (pale yellow).

Isolas: closed solution branches

Figure 5 has suggested that our snaking picture in the bistability zone is more complex than the usual one. In Figure 8 (left panel) the snaking picture observed before is thus plotted again, but this time the ensemble of all the solutions is drawn in the background in pale gray, while four solutions are highlighted using thick lines. It turns out that they have a well defined shape resembling a figure eight. Hereafter we will refer to these branches as isolas, since they are isolated. Despite the whole picture seems really intricate, it is made up of a number of such isolas. Closer study shows that our isolas are actually not strictly speaking isolated in the sense that there are no other solutions branching off, but that there are further bifurcations and interlinking solutions and connecting branches involved. We will come back to this point below.

In obtaining the snaking picture we found and plotted 53 isolas and connecting branches. Following the branches of the isolas, localised solutions with a different number of vibrating masses are linked together, which is conceptually linked to the snaking phenomenon. In Figure 8 four isolas are selected with the peculiarity to have all the masses moving in out-of-phase manner. On the left panel eight points, lying on the straight line  1 = 015, are marked with a red circle. The dimensionless vibration shape of each point is graphed in the correspondent bar plot on the right side. Focusing on the pair of points that belongs to the same isola, we can see that increasing the energy of the vibration leads to a higher number of masses involved: for example the first isola links solutions localised on one and three masses vibrating, the second isola links solutions localised on three and five masses vibrating, and so on. Compare the points 2 -3 or 6 -7: they belong to different isolas, nevertheless they are close in the bifurcation diagram. It can be seen that in fact those points have a slightly different vibration shape. Notice that the points 4 and 5 deserve some further attention: both of them involve five vibrating masses, but the solution has a different vibration shape (compare the bar plots of Figure 8). These considerations indicate how complicated the snaking phenomenon for our vibration system is, as not only the number of vibrating masses matters, but also the relative phases, i.e. the vibration shape plays a fundamental role.

Again, and just to double-check, in the bifurcation diagram (Figure 8, left panel) points obtained via time integration are denoted as red stars to show which branch of the isola is stable. Discrepancies between TI and HBM results are apparent and due to the 1 harmonic approximation, but small.

Isolas and connecting branches superposition

In this subsection we study in some more detail the bifurcation structures of the isolas found. In Figure 9 the left panel shows on the background the snaking picture (pale gray), while six solutions are emphasised to show on which part of the pattern they lie. We chose isolas in the bottom part of the pattern because they are less distorted than those in the upper part, and the picture is thus easier to explain. Nevertheless any other choice of the isolas would show the same main features. The six solutions are drawn in pairs, using a solid and a dashed line, into the subplots (a), (b), (c) to show to the reader how they are linked together. Looking at the subplot (a) one can see that the left hairpin bends of the isolas are connected through an independent branch which bifurcates from the isola. The same behaviour can be observed in Figure 9(b) where the branch drawn with the dashed line bifurcates from the isola drawn with the solid line. Figure 9(c) shows two solutions which are very close to each other along one branch. All these six solutions are then superposed to each other in Figure 9(d) to show how the complex overall snaking pattern emerges. Looking at the six solutions as a whole, it can be clearly seen that the snaking branches appear and even the peculiar structure of the 'ladder' (in the sense that they connect two different points of the same isola) is there, indicated by arrows in Figure 9(d).

In sum it almost seems that the usual snaking picture to be found in more strongly non-linear, and perhaps more strongly dissipative systems, is kind of broken into smaller elements, i.e. isolas and connecting branches. Further work on clarifying these aspects is necessary. In the subplot (d) all the six solutions are plotted together to show how they arrange in the overall snaking structure.

Vibration shapes for different solution branches

Here an analysis of the vibration shapes is carried out to show how the solution branches shown in Figure 9 relate to different vibration shapes. In Figure 10 we plot in the left panel the solution branches again (from Figure 9(a)) in the plane

³ e  max   1
´. On the isola of solutions (Figure 10) eight points labeled from 1 to 8 are marked with a red circle. We note that the points 8-1-2-3 belong to a part of the isola on which the vibration remains mostly localised on one single mass. After the solution number 3 a hairpin curve brings us to another part of the isola on which the solutions 4-5 show a vibration localised on three masses which vibrate in an out-of-phase manner. Then, for the solutions 6 and 7 the vibration of two of the masses reduces again and the vibration comes back to be localised one on just one mass.

We look now at the thick branch of solutions (Figure 10 left panel) which connects two points of the same isola (thin solid line). The solution 11 shows a vibration localised on one mass in a similar manner as for the solutions 7 and 8 on the previous isola. Moving towards the solutions 12-13-...-18 the neighbouring mass starts vibrating in phase with the previous one and the vibration remains localised on two masses (Figure 10) up to the intersection between the connecting branch and the isola. To sum up: while the isola comprises solutions localised on one or three masses, the connecting branch corresponds to solutions localised on one or two masses. Figure 11 shows in the left panel the branches presented in Figure 9(c). In this case the patterns have both a figure eight shape and they approach each other in the middle zone. The subplots on the right-hand side of Figure 11 give the shape of the vibration, while in Figure 12 the dimensionless energy of each mass e   is plotted for the points marked with squares/circles in Figure 11 (left panel). The isola of solutions drawn with a thin solid line connects solutions localised on two masses with others localised on three masses (Fig. 11 -Fig. 12). The thick solid line isola, instead, links solutions with three vibrating masses with others with four vibrating masses. Moving to the top of the isola, a smooth transition happens from one kind of localisation to the other. It is interesting to compare the solutions 4-5-6-7 with the solutions 18-11-12-13. The two branches are very close together but are not coincident. The solutions (18-4) (11-5) (12-6) (13-7) have the same number of masses involved in the vibration, but with a different shape. Figure 12 shows that even if the overall energy of the solutions (in pairs) is broadly the same, it is just distributed differently among the masses due to the different shape. For example if in the subplot 18 of Figure 12 the mass 7 is exchanged with the mass 8 the same distribution that appears in the subplot 4 is obtained, and the same reasoning can be done for the other three pairs, and so on for the whole two branches of the two isolas. Note that this is a feature that enriches the general picture of the snaking phenomena in vibrating systems as the appearance of the many different branches could be thought of to arise from some symmetries in the system. This could explain even why in the upper part of our snaking structure (Figure 5 left panel) the solutions appear more entangled: the larger the number of the masses involved in the vibration, the larger the possibilities to arrange them in different ways. ´. On each trajectory eight points are marked and labelled with a number. For each number the corresponding bar plot shows the shape of the vibration.

Conclusions

In this work we have studied snaking bifurcations of a non-linear cyclically symmetric oscillator chain. Bistability has been introduced in a heuristic manner through non-linear velocity dependent forces. Solutions have been obtained by time integration and harmonic balance techniques. The bifurcation diagrams resulting in the bistability zone resemble typical snaking patterns, but also show marked differences. The solution branches are composed of isolas, which have a figure eight shape in the bifurcation diagram. When the isolas are put together, they picture a typical snaking pattern, which has been observed in many other fields of physics. Still, our findings suggest that the snaking behaviour in structural dynamics could be more complicated due to the superposition of different non-linear mode shapes: solutions which have very different shapes present almost the same energy content and thus the corresponding solution branches overlap or are very close to each other in the bifurcation pattern.

Physically the snaking phenomenon is due to the weak nearest neighbor coupling among the oscillators and to the positive linear damping coefficient ( 1  0) which damps the small oscillations. In fact when one or few masses are in the high amplitude limit cycle the neighbors can be still stable under small oscillations and don't necessary jump on the high energy solution (limit cycle). We obtain therefore a multitude of solutions (from which, the snaking behavior) because of this bistability between low and high amplitude limit cycles.

From the present results it has become clear that more work on snaking in engineering structures is necessary. Future work will need to focus both on conceptual aspects, as well as on evaluating engineering relevance and impact. As for conceptual understanding, more work is necessary on understanding analogies and differences between the systems under study here and elsewhere. The larger number of states found in the present work needs to be clarified, and the relevance of symmetries, symmetry breaking and imperfect bifurcations, as observed here, needs to be understood better. As for engineering relevance, probably first another set of simplified and idealised model systems, ideally extracted from fluid-structure-interaction, frictioninduced vibration, or similar fields, will need to be studied. Thereafter more realistic models, e.g. derived from larger scale computer assisted modelling systems (like finite element analysis or computational fluid mechanics) could be studied.

As for our personal impression about future impact, a deeper understanding of snaking phenomena in engineering systems could make engineers able to better and quantitatively more accurately predict localised non-linear vibration states, which are known to be a source of numerous issues in engineering and technology, like fatigue, strength, or noise. For these reasons we hope that this work might serve as a starting point to conduct further studies on localisation phenomena and snaking in vibration engineering.

Figure 2 .

 2 Figure 2. The model system under study.

Figure 3 .

 3 Figure 3. (a) Dimensionless velocity dependent force plotted versus the velocity for 1 = [-04 -025 -01  06]  3 = 03  5 = 01  0 = 2  0 = 1 (b)Maximum potential energy of one single oscillator plotted against the bifurcation parameter  1  Blue dots: results of the time integration at the steady state. Solid line: HB solutions with a single harmonic. Dashed line: HB solutions with two odd harmonics.

Figure 4 .

 4 Figure 4. Mode shapes of the underlying linear system.

12 Figure 5 .

 125 Figure 5. Left: bifurcation diagram for the non-linear oscillator chain. The subcritical Hopf bifurcation of the individual oscillator is indicated, and the complex snaking pattern linking the spatially homogeneous stationary static state with the state where all oscillators are vibrating fills the zone of bistability. Middle and Right: the average dimensionless energy of each mass for the 12 equilibrium solutions which are marked with a red circle in the snaking pattern.

12 Figure 6 .

 126 Figure 6. In each subplot from 1 to 12 the vibration shape, i.e. with the corresponding phase information, of the solutions in Figure5(left panel, red circles) is shown.

1 ẼFigure 7 .

 17 Figure 7. Each one of the four columns corresponds to one solution obtained with TI. In the first row the displacement time history, in the second row the space state trajectories (  e vs e ), in the third row the dimensionless maximum potential energy of each mass in steady state condition is shown (dark blue). In the third row, the energy distribution obtained with single harmonic HBM is plotted (pale yellow).

Figure 8 .

 8 Figure 8. Left: The overall snaking pattern is kept in the background, while four isolas are highlighted. On each isola two points are marked with a red circle and a number. Stars indicate time integration results. Right: eight bar plots visualise the shape of the vibration for the corresponding solution marked.

Figure 9 .

 9 Figure 9. Left: Bifurcation diagram in the plane ³ e  max   1 ´ The snaking structure is left in the background, while six families of solutions are highlighted (thick line). On the right-hand side the six solutions are plotted in pairs (using one dashed and one solid line) into the subplots (a),(b),(c).In the subplot (d) all the six solutions are plotted together to show how they arrange in the overall snaking structure.

18 Figure 10 .

 1810 Figure 10. Left: two solution branches, or isolas, are drawn (Figure 9(a)) in the plane ³ e  max   1´. On each branch eight points are marked and labelled. For each point the corresponding bar-plot shows the shape of the vibration.

18 Figure 11 .

 1811 Figure 11. Left: two solution branches (Figure 9(c)) in the plane

18 Figure 12 .

 1812 Figure 12. Dimensionless energy plotted for each solution marked with the corresponding number from Figure 11.

  . is grateful to the Italian Ministry of Education, Universities and Research, which funded his PhD research. Part of the work has been supported by Deutsche Forschungsgemeinschaft in project HO 3852/11-1.

Also in Figure8(left panel) we will plot in the bifurcation diagram TI solutions superposed to the HB solutions to demonstrate that the first order truncation is able to capture the basic features of the exact solution. We are aware that this approximation can smooth out some fine details of the true solution, but the aim here is to focus on the overall result.