Backbone curves, Neimark-Sacker boundaries and appearance of quasi-periodicity in nonlinear oscillators: application to 1:2 internal resonance and frequency combs in MEMS

Quasi-periodic solutions can arise in assemblies of nonlinear oscillators as a consequence of Neimark-Sacker bifurcations. In this work, the appearance of Neimark-Sacker bifurcations is investigated analytically and numerically in the specific case of a system of two coupled oscillators featuring a 1:2 internal resonance. More specifically, the locus of Neimark-Sacker points is analytically derived and its evolution with respect to the system parameters is highlighted. The backbone curves, solution of the conservative system, are first investigated, showing in particular the existence of two families of periodic orbits, denoted as parabolic modes. The behaviour of these modes, when the detuning between the eigenfrequencies of the system is varied, is underlined. The non-vanishing limit value, at the origin of one solution family, allows explaining the appearance of isolated solutions for the damped-forced system. The results are then applied to a Micro-Electro-Mechanical System-like shallow arch structure, to show how the analytical expression of the Neimark-Sacker boundary curve can be used for rapid prediction of the appearance of quasiperiodic regime, and thus frequency combs, in Micro-Electro-Mechanical System dynamics.

Nonlinear dynamical phenomena in mechanical systems are often connected with the description of complex features that have no counterpart in linear theory, like jump phenomena, hysteresis, quasiperiodicity and chaotic vibrations [18,39,54]. In recent years, a number of studies appeared in several fields of Physics and Mechanics, highlighting the occurrence of Frequency Comb (FC) in the observed dynamical responses, especially when the system is driven by a single harmonic component. In the frequency domain a FC is characterized by a collection of spectral lines of frequency mX AE nx NS , where m and n are integers, X is the driving frequency and x NS is a new incommensurate frequency. In the time domain a FC is associated with a signal showing both amplitude and frequency modulation.
FCs have been largely investigated in optics where they have important applications in metrology, see e.g. [10,29,60,63]. A strong link also exists with the socalled modulation instability studied in the physics of nonlinear waves [9,64].
In recent years several occurrences of FCs have been documented in Micro Electro-Mechanical Systems (MEMS). These devices display a wide range of complex nonlinear phenomena, mainly due to their large quality factors Q, typically ranging from 10 3 to 10 6 . Indeed MEMS are often monolithic structures encapsulated in near vacuum packages and their dissipation is drastically reduced with respect to macrostructures. Internal resonance (IR) plays an important role in triggering more complex motions and facilitate energy transfer between modes and also the onset of quasiperiodicity [1,2,6,22,23,44,58,59]. Indeed in most cases, FCs in MEMS are associated with IR phenomena. FCs have been experimentally observed in [7,8] addressing a MEMS resonator with strong 1:3 IR between bending and torsional modes. Another example is provided in [31], where a resonator is proposed in which the dissipation can be dynamically eliminated, leading to the appearance of a tunable FC. In [45] the formation, evolution, and tuning of FCs in a piezoelectric MEMS resonator, based on nondegenerate parametric pumping, is presented. The tuning mechanisms of FCs are discussed with specific attention to the dependence of the center frequency and the frequency spacing between the spectral lines on the external input. In [22] a 1:2 IR between the first two symmetric vibrational modes of a MEMS arch resonator is investigated experimentally and theoretically. The MEMS structure is electro-thermally tuned and electro-statically driven. Here FCs are observed amid a variety of other complex behaviours. Other noteworthy studies on IR in MEMS arch resonators, with experimental data supported by theoretical models, can be found in [46,47]. In the context of MEMS piezoresonators, similar phenomena have been demonstrated in [14][15][16]. In particular, if the system has two modes with eigenfrequencies f 1 and f 2 such that f 2 % f 1 =2, an excitation can activate parametric resonance and trigger a FC. However, it is worth mentioning that FCs in MEMS might be associated with different phenomena. An example is discussed in [19] where FCs arise as a consequence of intermittent contact between deformable beams in a MEMS structure.
The issue of predicting the onset of a FC is strongly linked to the stability analysis of the associated periodic response. Once it gets unstable, a Quasi-Periodic (QP) regime might arise thus creating a FC. Restraining ourselves to the case of nonlinear vibrations, this problem can be reformulated as finding the locus of bifurcation points in periodic responses originating quasi-periodicity. In most cases these can be classified as Neimark-Sacker (NS) bifurcations [41,48], which represent the main focus of this work. It is however worth stressing that the onset of QP solutions and FCs have been also observed in [7,8] as a consequence of a Saddle-Node-Invariant-Circle (SNIC) bifurcation [51]. Besides, the stability of QP solutions developing on a torus is a fundamental aspect for questioning the appearance of chaos related to the stability of the period-three solutions [3,4,21,28,42].
The appearance of NS bifurcations in systems with IR has already been reported for systems featuring 1:1 IR, see e.g. [56] for arches or [57] for circular cylindrical shells, as well as for systems with 1:2 IR [39,40,54] such as arches [55]. Other examples involve different physical phenomena like, for instance, the aero-elastic behaviour of a flexible elastic suspended cable with linear eigenfrequencies in an almost 1:2 ratio driven by the wind speed [30]. Finally, NS bifurcations can also originate in more complex 1:2:2 and 1:2:4 resonance scenarios as reported in [36].
One can also remark that the computation of the locus of specific bifurcation points, as parameters are varied, has recently stimulated intensive research also on numerical techniques. For instance, continuation methods can be suitably adapted with additional constraints to follow the evolution of a given bifurcation point in the parameter space [11,61,62].
This work focuses on the appearance of QP solutions in nonlinear oscillators featuring 1:2 internal resonance. In this context NS bifurcation points have been known for a long time, see e.g. [35,37,38,40] and references therein. Consequently the viewpoint of the present contribution is to revisit a classical example in nonlinear oscillations, bring new results and use them in the perspective of better understanding and predicting the occurrence of frequency combs in MEMS dynamics. The main outcomes of the present study are to derive an analytical formula for the locus of NS boundaries, exhibit the two different families of periodic solutions connected to the two backbone curves of the system, named them as parabolic modes, and clearly highlight the links between the backbones, the frequency-response curves, the existence of isolated solutions and the locus of NS points. Since applications to MEMS-like structure are targeted, an arch structure is used in order to show how the findings can be used for fast prediction of the onset of FC in MEMS dynamical solutions.
The formulas derived analytically give simplified models that can be applied in systems that admit order reduction. In this work we consider systems in which it is possible to apply the center manifold theorem and more generally invariant and inertial manifold theory (see [24,33,34,52]). In this framework, strong but reliable simplifications can be applied and even continuous systems may be reduced to a limited number of degrees of freedom.
Future investigations will address the cases of 1:1 and 1:3 IRs, framing all these systems within the same setting and offering a unified view on the appearance of quasiperiodicity and link to backbones in assemblies of nonlinear oscillators featuring IR.
The paper is organized as follows. Sect. 2 is devoted to the analysis of conservative systems with the multiple scales technique. Two families of parabolic modes are identified and associated with two different backbone curves. The behaviour of these families of solution, for system parameters variation, is investigated. Sect. 3 is concerned with the analysis of the forced-damped system. A connection with the limit values of the conservative solution is given, and the appearance of isolas is explained. The NS boundary curve is computed and its dependence upon parameters and shape variation is analysed in detail.
Finally, in Sect. 4, a MEMS-like arch structure is analysed and the analytical expression of NS boundary curves is utilized to predict the occurrence of QP solutions, which represents the primary goal of the paper. Reduced order models using implicit static condensation [13,25,50] are obtained starting from the Finite Element (FE) discretization. The dynamics of the system is studied numerically through harmonic balance and direct time integration techniques in order to validate and refine the analytical estimate.

Conservative system: backbone curves
Let us consider the normal form of a system of two coupled nonlinear oscillators featuring 1:2 IR: in which q i , ði ¼ 1; 2Þ denote the displacement of each oscillator, x i the eigenfrequencies, a 12 ; a 11 the nonlinear positive quadratic coupling coefficients, l i , ði ¼ 1; 2Þ the linear damping coefficients and X the angular frequency of the forcing term with intensity F. Since we are interested in the 1:2 IR, the eigenfrequencies are such that x 2 % 2x 1 . In a conservative context, where internal forces derive from a potential, we have a 11 ¼ a 12 =2. However, since the results discussed in what follows can be used in more general contexts such as e.g. economy, chemistry or biology, in Appendix 1 we will relax this assumption and show examples of the behaviour for more general cases of quadratic coupling coefficients.
In this section, we address conservative dynamics and discard consequently the damping and forcing terms from the equations of motion. We apply the method of Multiple Scales (MS) to system (1) following [40]. Only the main results are recalled with emphasis on new findings. The detailed MS procedure is given in Appendix 2 for the sake of completeness.

Multiple scales solution
The method of MS expresses the solution as a composition of different time scales T j ¼ e j t for j ¼ 0; 1, and decomposes the solution as q i ðtÞ ¼ q i0 ðT 0 ; T 1 Þ þ eq i1 ðT 0 ; T 1 Þ þ Oðe 2 Þ with i ¼ 1; 2 and e a small bookkeeping parameter. We also assume that the nonlinear coefficients are small and can be expressed as a 12 ¼ e a 12 and a 11 ¼ e a 11 . Starting from Eq. (1) and introducing the MS expansions, we sort non-linearities according to e and we ignore the forcing F and the damping terms l i (see details in Appendix 2). The closeness of the fulfilment of the 1:2 ratio between the eigenfrequencies is quantified by introducing an internal detuning parameter r 1 such that We introduce a polar form for q 10 ðtÞ where c.c. stands for complex conjugate; a 1 ðT 1 Þ, a 2 ðT 1 Þ are unknown amplitudes; h 1 ðT 1 Þ, h 2 ðT 1 Þ are unknown phases, varying at the slow time scale T 1 and i is the imaginary unit. Developing the solvability condition (see details given in Appendix 2), we arrive at the modulation equations for amplitudes and phases: Inspecting Eq. (4), three different classes of possible solutions can be identified: 1) uncoupled solution with a 1 6 ¼ 0 and a 2 ¼ 0, denoted A1-mode in what follows; 2) uncoupled solution with a 2 6 ¼ 0 and a 1 ¼ 0, denoted A2-mode; 3) coupled solutions with both a 2 6 ¼ 0 and a 1 6 ¼ 0. The A1-mode is not admissible because of the invariant-breaking term a 11 q 2 1 in Eq. (1) considering the conservative condition. Indeed inserting q 2 ¼ 0 in Eq. (1b) implies that also q 1 ¼ 0. The A2-mode is admissible, since inserting q 1 ¼ 0 Eq. (1b) reduces to a linear oscillator, while Eq. (1a) is trivially satisfied.
Considering now the coupled solutions with both a 1 and a 2 different from zero, we simplify Eq. (4), dividing by the non-zero amplitudes a 1 and a 2 : To analyse the permanent solutions corresponding to fixed points of Eq. (5), the system needs to be made autonomous. We observe that the phases appearing in the sine and cosine functions involve the same quantity 2h 1 À h 2 À r 1 T 1 . This means that the system can be made autonomous with only two amplitudes and one phase. However, if one adopts such a choice an indeterminate quantity will appear when reconstructing the whole solutions for q 1 and q 2 , since both h 1 and h 2 will need to have an initial phase. To solve this issue, it has been preferred to make the system autonomous using two angular variables.
The resulting autonomous system reads Obviously, Eq. (7b) is redundant since only c p appears elsewhere in the system.

Coupled solutions: parabolic modes and backbone curves
The fixed point solutions of Eq. (7) exist only when sinðc p Þ ¼ 0 and consequently cosðc p Þ ¼ AE1 ¼ p, where the notation p ¼ AE1 is introduced for the rest of the paper. Consequently, the only possible coupled solutions correspond to limit cycles with a 1 and a 2 constants independent of time T 1 . The strong constraints imposed on the angular variables define two classes of solutions: one related to p ¼ 1 and one to p ¼ À1. We will show in the following that each of these solutions is associated to a family of periodic orbits with a given amplitudefrequency relationship (or backbone curve). To better describe the typology of these two solutions, let us first reconstruct the first order amplitudes q 10 and q 20 . Combining Eq. (3) with Eq. (6), the first order solutions are obtained as q 10 ðtÞ ¼ a 1 cosðc 1 þ x 1 tÞ; ð8aÞ q 20 ðtÞ ¼ a 2 cosð2ðc 1 þ x 1 tÞ À c p Þ: ð8bÞ When p ¼ þ1, then c p ¼ 2mp with m integer and the following relationship between the amplitudes holds: Eq. (9) dictates that the amplitudes lie on a parabola with positive concavity in the configuration plane ðq 10 ; q 20 Þ. In the whole four-dimensional phase space (including both velocities), we hence obtain a family of periodic orbits developing on a manifold such that its projection on ðq 10 ; q 20 Þ is a parabola. This solution is called parabolic mode, following also the designation given in [6,17,32] for the case of 1:1 resonance, where normal modes and elliptic modes are the two families of coupled solutions. More precisely, in the case p ¼ þ1, this solution is denoted as p þmode. When p ¼ À1 a second family of periodic orbit is obtained, where c p ¼ ð2m þ 1Þp, with m integer. The amplitudes now define a parabola with negative concavity in the configuration plane ðq 10 ; q 20 Þ: denoted as p À -mode . These two families of solutions are represented in Fig. 1. It should be noted that for the p þ -mode q 1 and q 2 have the same phase, while for the p À -mode the phase difference is equal to p.
In order to express the frequency-amplitude relationship (backbone curves) for these two families of coupled solutions, we elaborate on the relationship between the amplitudes a 1 and a 2 given by the fixed point equation associated to Eq. (7d): where we have replaced the detuning variable r 1 with its expression in terms of eigenfrequencies x 1;2 to have a more straightforward insight into the system. Denoting with q i the maximum amplitude of q i ðtÞ ( q i ¼ a i ; i ¼ 1; 2), we now replace the non-linear coefficients with the original a 12 and a 11 , and we express q 1 in terms of q 2 : which holds for both parabolic modes, setting either (b) (a) Fig. 1 Sketch of the parabolic modes corresponding to the coupled solutions in the configuration plane ðq 10 ; q 20 Þ. a) p þmode, b) p À -mode p ¼ 1 or p ¼ À1. To obtain the amplitude-frequency relationship we first solve the first-order equations at the slow time scale T 1 , Eq. (7), and get the following expressions for the phase angles: with / 1 and / p integration constants. Using Eq. (6), the original nonlinear coefficients a 12 and a 11 and using T 1 ¼ et we have: Inserting these equations in Eq. (3) we get From Eq. (15), substituting q 1;2 , one obtains the nonlinear oscillation frequencies x NLi of each oscillator: ; ð16aÞ We remark in particular that the relationship x NL2 ¼ 2x NL1 is always fulfilled, as can be shown inserting Eq. (12) into Eq. (16b). For the sake of clarity, x NL1 will be simply denoted x NL in what follows. Note that Eq. (16a), together with Eq. (12), identifies the solution manifold in the space ðx NL ; q 1 ; q 2 Þ where the backbone solutions are lying. The two families of solutions are discriminated by replacing p with AE1 in all expressions. Since the amplitudes q 1 , q 2 are assumed to be positive, one arrives at the conclusion that: when x NL \x 1 , only the p À mode exists, while when x NL [ x 1 , one only has the p þ -mode. This gives a first condition for the existence of the backbone curves.
An additional condition can be derived by inserting Eq. (16a) in the rhs of Eq. (12) that must be positive to guarantee the existence of a solution for q 1 : In Eq. (17), the third term is always positive due to the assumptions on a 12 . The first two factors of Eq. (17) impose that q 1 exists only for values of x NL outside the limit points The existence regions of each parabolic modes are summarized in Fig. 2 in terms of the shift of the nonlinear oscillation frequency x NL À x 1 (on the xaxis) and of the detuning x 2 À 2x 1 (on the y-axis). For a given system, the detuning is fixed so that a single horizontal line gives the existence region. The blue region (respectively the red region) corresponds to the existence of p À (resp. p þ ) mode. From Eq. (17), the lines x NL ¼ x 2 =2 (highlighted with a green color) and the y-axis x NL ¼ x 1 , delimit a region in which no coupled solutions exist. The green hatched region delimited by It is worth looking at the behaviour of the backbones at the limit point of their existence regions. At the first boundary x NL ¼ x 1 we get: Fig. 2 Schematic representation of the existence region for the coupled parabolic modes, in the plane ðx NL À x 1 ; x 2 À 2x 1 Þ. The abscissa represents the shifted nonlinear oscillation frequency while the ordinate is the detuning between the eigenfrequencies of the system lim meaning that close to x 1 , the backbone curve amplitudes are going to zero as expected. The situation is different in the vicinity of x 2 =2. Indeed at this second boundary In particular the limit value of q 2 is non-zero at the limit point of the backbone when x NL À! x 2 =2, which could appear conflicting with the coupled solution assumptions. One should note that this is only true in the limit, and that the uncoupled A2-mode solution still exists for the system. Consequently the limit point of the backbone connects with this solution at this specific point. This particular situation induces a jump in the backbone solution (see the detuned condition in Fig. 3 commented in Sect. 2.4) and, as shown in Sect. 3.2, will also lead to the emergence of isolated solutions in the forced and damped system.

Stability of coupled solutions
The stability of coupled solutions is governed by the eigenvalues of the Jacobian matrix J of the conservative system (see Appendix 2 for details). The determinant det ðkI À JÞ of the characteristic polynomial, expressing the angular variables in terms of q i ; i ¼ 1; 2, writes: The nonzero eigenvalues k are given by: and have negative real part if: Backbone curves of the two families p þ (red curves) and p À -modes (blue curves), in the plane ðx NL ; a 1 Þ (first row) and ðx NL ; a 2 Þ (second row), and for three different possible cases. First column: x 2 =2\x 1 , second column: The unstable region delimited by conditions in Eq. (22) is shown in Fig. 2 as the green hatched region. We observe that this region completely falls in a non-existence area of coupled solutions. Consequently, the conservative backbone curves are always stable.

Backbone curves
In this section, we show how the backbone curves of the two different parabolic modes organize when the detuning parameter is varied. Note that some of the results presented here are close to those presented in [27], where a different point of view was adopted on the same system. Fig. 3 shows the three different possible cases: when the detuning is negative, when it is vanishing, and when it is positive. The figure plots the projection of the solutions in the planes ðx NL ; q 1 Þ and ðx NL ; q 2 Þ. However, an optimal representation would be in the 3D space ðx NL ; q 1 ; q 2 Þ, as will be shown in Sect. 3 when the solutions will be fully reconstructed.
One can note that, consistently with Eq. (12)-(16), the backbones in the ðx NL ; q 2 Þ plane are always straight lines, while their projection in the ðx NL ; q 1 Þ plane are line segments only when the detuning vanishes, and otherwise have a parabolic shape. While increasing the detuning from negative to positive, one can observe the following features. When x 2 =2\x 1 , the backbone starting at x 1 is a p þ -mode, while the backbone starting at x 2 =2 is a p À -mode. According to Eq. (19b), the p À -mode backbone starts from a nonzero value in ðx NL ; q 2 Þ plane. In the case of perfect 2:1 resonance, x 2 =2 ¼ x 1 , the backbones of the two families start from the same point and from zero amplitude. Symmetrically, when x 2 =2 [ x 1 , the backbone emanating from x 1 is now a p À -mode, meaning that in the process of the transition a switch of family occurs. Consequently, the family of periodic solutions emanating from x 2 =2 is now a p þ -mode, and it starts at a non-zero value in the ðx NL ; q 2 Þ plane.

Forced system and NS boundaries
This section investigates the Frequency Response Functions (FRF) of the oscillators described by Eq. (1), with a forcing term applied on the lowfrequency oscillator and damping. While the appearance of quasi-periodic solutions has been already documented for this case e.g. in [40,55], our aim here is to provide an exact expression of the NS boundary curve. Note also that quasiperiodic solutions are usually not addressed in the case where the forcing is applied on the high-frequency mode since the main concern is typically to investigate the loss of stability of the uncoupled solution [40,53]. As for the conservative case, first-order perturbative solutions are derived with the MS method. Since this derivation is classical, most of the details are reported in Appendix 3, and emphasis is rather set on the derivation of the NS boundary curve, its dependence on parameters and the link between the FRF and the backbones.

Multiple Scales solution
First, the equations of motion in Eq. (1) are rewritten considering the MS expansion for time and oscillators degrees of freedoom using a bookkeeping parameter e (see Appendix 3). We also introduce the assumptions that the non-linearities, the forcing amplitude and the damping are small, and can be expressed as a 12 ¼ e a 12 , a 11 ¼ e a 11 , l 1 ¼ e l 1 and F ¼ e F. The forcing angular frequency is in the vicinity of the first eigenfrequency x 1 , and an external detuning parameter r 2 is introduced to describe this nearness as: The first-order system governing the modulation of the amplitude is derived in a similar way as in the previous section. At the slow time scale T 1 , the resulting system reads: When the forcing is applied on the low-frequency oscillator, only coupled solutions exist (see Appendix 3).
The system is made autonomous by introducing c 1 and c 2 such that: and reads: It is worth highlighting that the definitions of Eqs. (25) and (6) and the structure of Eqs. (26) and (7) are similar and only differ for the additional terms associated to F, l i and r 2 . The fixed points, associated to forced oscillations of constant amplitudes, can be expressed as function of the maximum amplitudes of ð q 1 ; q 2 Þ only. Internal and external detunings r 1 and r 2 can also be replaced by the frequency differences to have a more straightforward insight into the system. The amplitude equations, considering the original system parameters, reads: ð27bÞ can be seen as the solution manifold in the space ðX; q 1 ; q 2 Þ. Eq. (27a) is a polynomial of the third-order in q 2 and can give, depending on the parameters, from one to three solutions. The phases of fixed points are given by: These equations are not strictly necessary to build FRFs, but are needed to reconstruct the amplitude variations in time. They can also be used to establish existence conditions for the solutions exploiting the bounds of the trigonometric functions. Fig. 4 shows the type of solutions obtained when varying the frequency detuning and the link with the backbone curve and the p þ or p À modes. Note that in Fig. 4 emphasis is put on the shape of the solution, while stability is addressed in the following section.
In the first line, the case of a positive detuning between the eigenfrequencies is considered, In particular one can observe that close to X ¼ x 2 =2, the FRF has a local minimum for a 1 and a maximum for a 2 . The resonant branches follow the backbones of the two families of periodic orbits p þ and p À . Fig. 4 c) depicts the solution in the configuration plane ðq 1 ; q 2 Þ for the three points marked with coloured circles on the FRF. As expected, the blue point corresponds to p À -mode and one recovers the shape of the parabola in ðq 1 ; q 2 Þ, with a slight dephasing corresponding to the addition of forcing and damping (exact parabola is retrieved for the backbone only). The same reasoning applies for the red point corresponding to p þ -mode. The black point is selected exactly at X ¼ x 1 and is a p À -mode also.
The second line of Fig. 4 considers the tuned case Fig. 4 f), while p þ and p À mode are still perfectly well retrieved, the black point at X ¼ x 1 shows a combination of the two solutions resulting in a symmetrical shape in the ðq 1 ; q 2 Þ plane. Finally, the third line shows the case of a negative detuning, and behaves symmetrically as compared to the case of positive detuning. Increasing the detuning from negative to positive values, we see how the solution emanating from the forced low-frequency x 1 switches from the p À to the p þ mode, with the intermediate step being the symmetric shape shown in Fig. 4 f). This result is in line with those presented in [27] and complement their analysis by underlining the parabolic nature of the two families of periodic orbits, with different curvature.

Isolas
Interestingly, the first-order solution obtained by MS analysis contains isolated branches of solutions, as Typical frequency-response functions and orbit in configuration space ðq 1 ; q 2 Þ, for different cases of detuning. First line, positive detuning, Stability is not reported. Selected values of the parameters: already remarked in [27]. In addition to their analysis, it is here underlined that the manifestation of these isolas is the consequence of the non-vanishing starting point on the q 2 axis when the internal detuning is different from zero. This particular feature is illustrated in Fig. 5 for the case of a positive detuning, x 2 =2 [ x 1 . As shown in the previous section, the limit value of the p þ -mode corresponds to a nonvanishing value of q 2 . Since the forced-damped solutions emanate from the backbone curves [5], this eventually creates the condition for the appearance of isolas.

Stability analysis and Neimark-Sacker boundary
The stability analysis of the solutions is based on the eigenvalues k of the Jacobian matrix of Eq. (26), analysed in Appendix 3. The characteristic polynomial can be written as: with all the coefficients detailed in Appendix 4. Since all expressions are available in closed form, one can explicitly compute the stability of the solutions and follow bifurcation points. We focus on the definition of the Neimark-Sacker (NS) boundary curve. The NS bifurcation requires that a pair of eigenvalues are complex conjugate with zero real part. Using the compact definition of the characteristic polynomial of Eq. (29), this implies that the four eigenvalues are iA; ÀiA; B; C where A 2 R and B; C 2 C. Then Eq. (29) can be rewritten as [26]: Equating the same power terms one obtains a system of four equations: Isolating A 2 as the ratio between Eq. (31a) and Eq. (31c), expressing the product BC using Eq. (31d) and inserting both terms in Eq. (31d), the condition for a NS bifurcation can be expressed in terms of the c i 's only, reading:  Fig. 6 FRF, backbone curves and NS boundary when x 1 ¼ 1, Inserting the detailed expressions of each coefficient in Eq. (32), the definition of the NS boundary curve is obtained as a polynomial in q 2 , and reads:  Fig. 7 FRF, backbone curves and NS boundary when and a 11 ¼ 2:5 Á 10 À2 . The black continuous lines are the FRF (forced and damped solutions). The red and the blue lines are the backbones corresponding to the p þ and p À mode respectively. The black dashed line marks the X ¼ x 2 =2 frequency.
The green line is the NS boundary. In 3D figure in the space X À q 1 À q 2 the orange surface is the manifold on which all the possible forced and damped solution lies and is given by Eq. (27b). Selected values of the forcing amplitude are F ¼ 2:5 Á 10 À2 ; 5 Á 10 À2 ; 7:5 Á 10 À2

FRFs, backbone curves and NS boundary
We now discuss the complete solutions including the FRF, their connection with the backbone curves and the location of the NS boundaries.
We select a fixed set of parameters x 1 ¼ 1, l 1 ¼ l 2 ¼ 1 Á 10 À3 , a 12 ¼ 5 Á 10 À2 and a 11 ¼ 2:5 Á 10 À2 . First, the case of a vanishing internal detuning, thus Fig. 6 collecting the FRF, Eq. (27), the NS boundary, Eq. (33); and the backbones given by Eqs. (16) and (12). Fig. 6 a) is a view in the 3D space ðX; a 1 ; a 2 Þ. The solution manifold given by Eq. (27b) is represented as an orange surface and all solution branches as well as the NS boundary curve lie on this manifold. In this case the figure is fully symmetric with respect to X ¼   The green line is the NS boundary. In 3D figure in the space X À q 1 À q 2 the orange surface is the manifold on which all the possible forced and damped solution lies and is given by Eq. (27b). Selected values of the forcing amplitude are F ¼ 2:5 Á 10 À4 ; 5 Á 10 À4 ; 1 Á 10 À3 b), c) show the projection of the 3D plot respectively on planes ðX; q 1 Þ and ðX; q 2 Þ. The FRF develops according to the backbone curves of the p þ and p À modes. When X approaches x 1 the FRF have a local peak for q 2 , and correspondingly q 1 approaches a nearly zero value. In this case, the NS boundary is symmetric and displays a unique minimum centred at X ¼ x 1 , as shown in the enlarged view Fig. 6 d) and e). Fig. 7 shows the case of a positive detuning, with  Fig. 9 Behaviour of the NS boundary curve when the detuning r 1 ¼ x 2 À 2x 1 is changed. The system parameter here considered are a 12 ¼ 1 Á 10 À2 , a 11 ¼ 5 Á 10 À3 , l 1 ¼ l 2 ¼ 5 Á 10 À3 . The square markers correspond to the minima predicted by the analytical expression given in Appendix 5, namely X ¼ x 2 =2 and X ¼ ðx 1 þ x 2 Þ=3 ( 2 )/2 ( 2 + 1 )/3 ( 2 )/2 Fig. 10 Behaviour of the NS boundary curve, for x 1 ¼ 1, x 2 ¼ 2:05, a 12 ¼ 1 Á 10 À2 , a 11 ¼ 5 Á 10 À3 and different values of l 1 ¼ l 2 as detailed in the legend. The blue and red continuous lines are the backbone of the p À and p þ mode while the black dashed lines correspond to X ¼ x 2 =2 and unchanged. In this case the symmetry is broken, as can be appreciated by inspecting the shape of the solution manifold in the 3D space, showing a minimum with respect to q 1 at X ¼ x 2 =2. The solution branches arising in the vicinity of X ¼ x 1 belong to the p À family, in line with previous remarks. The shape of NS boundary curve (see Fig. 7d) and e)), also shows important modifications, having two minima and one local maximum, particularly visible in the ðX; q 1 Þ projection. Interestingly, the projection on the ðX; q 2 Þ plane shown in Fig. 7c) underlines the closeness of the NS boundary to the non-zero point where the backbone of the p þ -mode emerges.
Finally, the case of a negative detuning with x 2 ¼ 2x 1 À 5 Á 10 À2 is addressed in Fig. 8. The solutions around x 1 now belong to p þ -mode. The NS boundary shows again the two minima and the local maximum, and its closeness to the emerging point of the p À -mode is symmetrically observed.
These results underline that a further investigation of the minima and maxima of the NS boundary curve, together with its relationship with the backbones, would provide a better understanding of the dynamics. This is the aim of the following sections. ( 2 )/2 ( 2 + 1 )/3 p + p - Fig. 11 Behaviour of the NS boundary curve, for x 1 ¼ 1, x 2 ¼ 2:05, a 12 ¼ 5 Á 10 À2 , a 11 ¼ 2:5 Á 10 À2 and different ratio l 1 =l 2 keeping l 1 ¼ 5 Á 10 À3 fixed. The blue and red continuous lines are the backbone of the p À and p þ mode while the black dashed lines mark X ¼ x 2 =2 and X ¼ ðx 1 þ x 2 Þ=3. The other lines corresponds to different damping levels as detailed in the legend

NS boundary behaviour and extrema
In this section the parametric dependence of the NS boundary curve with respect to detuning and damping is investigated with specific attention to special points corresponding to extrema of the boundary. Fig. 9 first displays the shape of the NS boundary curve when the detuning is varied from negative to positive values, in both planes ðX; q 1 Þ and ðX; q 2 Þ. The plots have been obtained for x 1 ¼ 1, l 1 ¼ l 2 ¼ 5 Á 10 À4 , a 12 ¼ 1 Á 10 À2 and a 11 ¼ 5 Á 10 À3 and different values of x 2 . As already underlined, the boundary curve has a single minimum in the case of no detuning, obtained from Eq. (33) assuming x 2 ¼ 2x 1 ¼ 2X: For non-vanishing detuning two minima occur in the plane ðX; q 1 Þ, and they correspond to points where there is a change of curvature in ðX; q 2 Þ projection. This result is of prime importance since the locations of the minima convey meaningful information about the minimum forcing and amplitude levels that are needed to reach a quasiperiodic solution.
In Fig. 9, two points for each curve are marked with a black square, corresponding to X ¼ x 2 =2 and X ¼ ðx 1 þ x 2 Þ=3. These are minima of the NS boundary curve predicted by asymptotic analysis as detailed in Appendix 5.
These specific points are obtained in the limit of small damping and are thus only approximate values. This is further illustrated in Fig. 10, obtained for x 1 ¼ 1, x 2 ¼ 2:05, a 12 ¼ 1 Á 10 À2 , a 11 ¼ 5 Á 10 À3 and different values of damping assuming l 1 ¼ l 2 . The detuned frequency condition is considered because the loci of the minima change only in this case. The figure shows that the NS boundary tends to have sharper minima and values closer to the proposed approximation as damping tends to zero. On the other hand, the behaviour of the NS boundary curve at large values of amplitudes is asymptotically identical and not very sensitive to variations of the damping.
A different case is analysed in Fig. 11, obtained for x 1 ¼ 1, x 2 ¼ 2:05, a 12 ¼ 5 Á 10 À2 , a 11 ¼ 2:5 Á 10 À2 and different values of the ratio l 1 =l 2 , while keeping l 1 ¼ 5 Á 10 À3 fixed. The figure shows that increasing the damping l 2 makes the NS boundary smoother. It moreover underlines that in this specific case the asymptotic behaviour of the curve at large amplitudes  (1) when the detuning parameter vanishes, i.e. x 2 ¼ 2x 1 , and in the general case x 2 6 ¼ 2x 1 , the latter condition being more realistic from a practical point of view since process imperfections typically prevent a perfect match between the two eigenfrequencies. Finally a comparison between a reduced-order model including all nonlinearities, the corresponding MS solution and a full Finite Element Simulation is proposed using a custom code [43].

Reduced order model
The arch geometry and dimensions are reported in Fig. 12.
The device is made of polycrystalline silicon with density q ¼ 2330 kg/m 3 . A linear elastic constitutive model (Saint-Venant Kirchhoff) is considered, with Young modulus E ¼ 167000 MPa and Poisson coefficient m ¼ 0:22 [49]. The double-arch has a large radius of curvature and it can be considered as shallow. The structure is inspired by the design proposed in [12], and it exploits two clamped arches connected at their center points in order to shift the anti-symmetric modes to higher frequencies. The first six eigenfrequencies of the structure, obtained from a Finite Element (FE) analysis using the mesh in Fig. 12 b, are reported in Tab. 1, and the corresponding mode shapes are shown in Fig. 13. In particular, one can observe that it is possible to achieve a 1:2 resonance between the first and the sixth eigenmodes, Fig. 13 a and f.
A reduced-order model (ROM) is obtained from the full FE model using a static condensation approach [13,25]. The ROM retains as master modes only the two eigenmodes involved in the 1:2 resonance. The approach uses static loadings proportional to the inertia of the master modes to create a stress manifold ðkÞ ijs are the linear, quadratic and cubic coefficients, respectively, for eigenmode k. The load multiplier considered in our benchmark is F ¼ 0:0201 lNls 2 /ng. The first and second order derivatives with respect to the nondimensional time are denoted with 0 and 00 respectively. Even if X is now a nondimensional frequency, the same symbol of the previous sections has been used since we fixed x 1 ¼ 1 in all the examples presented therein.
The damping factor is l i ¼ x i Q i x 1 with Q 1 and Q 2 factors set to 500 and 1000 respectively, which is consistent with a mass proportional damping. The coefficients of the polynomial are reported in Tab. 2.
Arch displacements can be reconstructed with good accuracy from the linear modes and in particular the midspan deflection is uðtÞ ¼ 0:11q 1 ðtÞ þ 0:0707q 2 ðtÞ ½lm: ð36Þ We notice that the selected damping values are slightly lower than the ones assumed in the previous examples concerning the MS solutions. This condition is typical of MEMS structures that are often packaged in near vacuum. The impact of these choices will be discussed in the following sections.

Multiple scales model
In order to apply the MS formulas, the nonlinear terms of Eq. (35) that are not included in Eq. (1) are neglected. First we consider the perfectly tuned case by enforcing x 2 ¼ 2x 1 . We compare the FRF given by Eq. (27) and the solution computed via direct numerical continuation of periodic orbits applied to Eq. (1) using the continuation package MANLAB [20]. The numerical continuation solution is also cycles are simulated to reach the steady state regime and collect time history data. The time step adopted is 1/500 of the forcing period 2p=X.
Results are plotted in Fig. 14a). The MS and continuation (MANLAB) FRF are superimposed everywhere except close to the peaks. The numerical solution of Eq. (1) predicts a larger amplitude at the lower frequency peak and a smaller amplitude at the high frequency one. Indeed at the peaks the assumption of the MS solution that nonlinear terms are small starts being violated. On the contrary, the NS boundary intersections with the FRF predicted by the MS have a nearly perfect agreement with the bifurcation analysis performed by MANLAB (see the enlarged view in Fig. 14 b). This is reasonable since the NS boundary crosses the FRF at low amplitudes, where the MS approximation is fully respected.
Inspection of the time-marching solution reveals additional information. The match with the FRF obtained thanks to numerical continuation is almost exact in the periodic regions, the only difference being that the RK4 results present jumps in the solution (highlighted with arrows in the plot) since unstable solutions cannot be simulated. The time-marching solution allows appreciating amplitude modulation in the QP region. For a fixed frequency value a cloud of points is observed, as expected in a QP regime. The incommensurate frequency can be estimated from the Fourier Transform (FT) of the response and results are reported in Table 3 for both the upward and downward span.
The MS solution developed in this work cannot predict the behaviour of x NS inside the QP region, but it can be used to estimate the values at its boundary, since here the purely imaginary eigenvalues define the oscillating frequency of a i . The NS bifurcations occur at X ¼ 0:99838 and X ¼ 1:00161 respectively, corresponding to the eigenvalues k ¼ AEi5:49 Á 10 À3 , that have the same value at both frequencies. The computed analytical value of k can be compared with the numerical results of Tab. 3. It is worth stressing that the agreement is good at the left of the NS in the upward sweep and at the right of the NS in the downward sweep, i.e. at the onset of the instability. The agreement then changes proceeding in the NS region since the numerical solution follows the QP solution branch on which the additional frequency of the torus is varying, a feature not provided by the MS analysis. Examples of the FT and time histories encountered in the frequency span are illustrated in Fig. 15 a and b referring to X ¼ 0:99838 and X ¼ 0:98711 respectively.
In Fig. 15 a, the QP regime has just settled down and one can observe a clear modulation of the envelope in the time domain, and a well shaped spectrum with clean spacings between all the frequency peaks, following the rule nX AE mx NS with n and m integers. Fig. 15 b has been obtained for X ¼ 0:98711 in the upward frequency sweep, thus the QP regime is further along the branch of QP solution. One can appreciate how the QP is modified along the solution branch, with a more complex envelope modulation, and a more important number of frequency peaks in the spectrum. One can also observe the enlargement of the frequency peaks, indicating that the QP solution will probably lose its stability in favour of a chaotic behaviour further in the solution branch. One can also note that the direct time integration results underline the facts that the NS bifurcation is supercritical with emergence of a stable torus at both side of the branching points, symmetric, with probably a small portion of unstable torus close to the branching points.
Before moving to the detuned condition, phase space representation from time-marching solution are inspected.
The torus emerging from the NS bifurcation is represented in Fig. 16, obtained from direct RK4 integrations with X ¼ 1. A 3D representation in the space (q 1 ; _ q 1 ; q 2 ) and in the planes ðq 1 ; _ q 1 Þ and ðq 1 ; q 2 Þ are shown. For sake of clarity, only a short portion of the time history, corresponding to few quasi-periods, is considered. Fig. 16 a represents the phase space ðq 1 ; _ q 1 ; q 2 Þ while Fig. 16b) is a cross-section of the four dimensional torus considering the portion included in the green and the blue planes. The red dots mark the crossing region. Figure 16c) and d) are (1) in the tuned case x 2 ¼ 2x 1 . Only a short portion of the time history is plotted for sake of clarity the 2D version of the phase space considering its two projections in the planes q 1 ; _ q 1 and q 1 ; q 2 . The detuned case is now investigated with x 2 ¼ 1:989x 1 , and the same analysis is repeated. The FRF results are plotted in Fig. 17 and the x NS values are collected in Tab. 4. By comparing the MS and the numerical continuation (MANLAB) solutions, the FRFs and the NS bifurcation crossing points share the same degree of accuracy of the previous case. Focusing on the incommensurate frequencies, the MS solution predicts x NS ¼ 4:37 Á 10 À3 at X ¼ 0:99372 and x NS ¼ 6:19 Á 10 À3 at X ¼ 0:99679. Fig. 17b) shows how the QP time-marching solution is completely embedded in the region estimated by MS. An outstanding match is found close to X ¼ 0:99679 while the values on X ¼ 0:99372 are lower than expected. The underlying reason is that at X ¼ 0:99679 a supercritical NS bifurcation appears, which is also confirmed by the identical values of x NS found close to this point when sweeping the forcing frequency in the increasing or decreasing direction. On the other hand, a subcritical NS bifurcation takes place at X ¼ 0:99679, consequently the added frequency predicted by the MS method corresponds to unstable quasi-periodic solutions that are not found with direct time integration. Instead, the numerical solution jumps to the branch of stable QP solutions and then travels up to the other supecritical bifurcation point.

Reduced order model and full order simulation
Finally we consider a comparison between the MS solution and the ROM retaining all nonlinear terms. Results obtained with the ROM are further validated with analyses performed on the full FE model using two different techniques. First, an Harmonic Balance approach (HB) is used, from which a reference FRF is computed. However, this technique in our implementation cannot detect QP solutions. As a consequence, a costly direct time integration of the full FE model is also performed with a Newmark-b scheme for specific values of X. The time-marching solution adopts a time step equal to 1/300 of the forced eigenmode period and 3000 cycles are simulated for each frequency. The full FEM model uses as external forcing a body force having the same shape of the first eigenmode, which is consistent with the ROM loading in Eq. (35). The comparison between the MS solution and the ROM is plotted in Fig. 18. As expected and mostly because of the presence of the cubic terms in the ROM, some quantitative differences appear between the two solutions. However no qualitative difference is reported since the simplified system studied in the MS development contains the most important resonant monomial terms that convey the important bifurcation information. Moreover in the low amplitude region and close to X ¼ 1 the MS solution is nearly exact. Also the prediction of the NS boundary intersections with the FRF shows a very good match with numerical results on the ROM. This is expected since the MS approach is reliable in a low amplitude region where the smallness assumptions still hold, and thus underlines that the analytical formula of the NS boundary curve can be used for rapid prediction of the occurrence of FC in such a case.
Finally we compare the ROM solution with a Full FEM simulation considering the maximum absolute midspan displacement max juj. The results are plotted in Fig. 19 where the NS boundary curve from the MS solution is also included. We notice the excellent agreement between the FRF obtained using numerical continuation and HB for the full FEM model and for the ROM, respectively, validating the reduction technique adopted. As expected, the NS boundary from the MS predicts correctly the onset of the QP regime also in this case. Fig. 20 shows three representative quasi-periodic solutions obtained from direct time integration on the full FE model. Note that these solutions are obtained from a decreasing sweep so that the solutions follow the QP branch of solution by order of decreasing frequencies. The first line in Fig. 20 corresponds to X=0.9965, which is the onset of the QP regime, just after the NS bifurcation point. The amplitude modulation displays a simple pattern in Fig. 20 a) and the frequency spectrum shows the appearance of the extra peaks indicating the birth of the FC. At this point one can compare the x NS value predicted by the MS at the boundary (see Table 4) and that from from direct numerical integration. The comb spacing in Fig. 20b reveals the value x NS ¼ 5:93 Á 10 À3 , very close to the values reported in Table 4 and thus underlining again the good predictive capacity of the simple analytical model. The second line in Fig. 20 corresponds to X=0.9956, a point further along the branch of QP solutions. One can observe that the amplitude  Fig. 18 Comparison between the analytical solution (black line) and the static condensation ROM (red line). The analytical NS boundary is the green dash-dotted line. The red star markers define the bifurcation points identified by MANLAB modulation is more complex but still clearly periodic, resulting in a well defined frequency comb. In this case the new frequency x NS cannot be determined from the MS analysis. Finally, the third line shows the results obtained for X=0.9939. In this case one can clearly observe that the envelope modulation has no clear periodicity at the reported time scales, underlining that the torus is on the way to lose its stability by creating longer and longer periods which will result in a chaotic solution. The frequency spectrum is more densely filled. A complete study of the chaotic nature of this solution could be performed by computing the Lyapunov spectrum or a more involved stability analysis of the torus solution. This point is however not the scope of the present study which is focused on the appearance of quasiperiodic solutions and their link to Frequency combs in MEMs dynamics.

Conclusion and further developments
In this paper we have developed analytical results for 1:2 internally resonant coupled oscillators. Even though this analysis is classical in nonlinear vibration theory, new insights have been reported, in particular with regard to the identification of the two backbones and their corresponding parabolic mode shapes and the analytical derivation of the NS boundary curves. All these results have been also put under the frame of explaining the appearance of FC in MEMs dynamics, a subject that give rise to numerous investigations in the recent years, and all the results have been applied to a MEMS arch structure in order to underline the predictive capacity of the analytical NS boundary curve. This result can be helpful in a context of MEMs design for accurate and fast predictions since the QP solutions are sometimes intentionally searched for.
Starting from the normal form of the dynamical system we inspected the conservative solution providing a clear identification of the mode shapes and solution branches. By extending our study to the forced and damped condition we were able to retrieve a closed form solution for the NS boundary curve. We inspected in depth the behaviour of the bifurcation boundary by varying the system parameters and comparing it with the conservative and the nonconservative solutions. We proved the validity of our findings in the challenging simulation of a MEMS structure example. The analytical solution has been validated against several numerical methods including numerical continuation procedures (MANLAB), reduced order models solved with direct time integration or the Harmonic Balance Method and eventually also with full FEM approaches.     Fig. 22 Behaviour of the backbones when a 11 =a 12 is changed keeping fixed a 11 ¼ 1 Á 10 À2 . In these plots the line color marks with blue the p À and with red the p þ mode, the black dashed line marks X ¼ x 2 =2. Figure a) and b) represent the backbone q 1 and q 2 respectively when x 2 À 2x 1 ¼ 0. Figure c) and d) represent the backbone q 2 and q 2 respectively when x 2 À 2x 1 ¼ 5 Á 10 À3  Fig. 21 Behaviour of the backbones when a 12 =a 11 is changed keeping fixed a 12 ¼ 1 Á 10 À2 . In these plots the line color marks with blue the p À and with red the p þ mode, the black dashed line marks X ¼ x 2 =2. Figure a) and b) represent the backbone q 1 and q 2 respectively when x 2 À 2x 1 ¼ 0. Figure c) and d) represent the backbone q 1 and q 2 respectively when x 2 À 2x 1 ¼ 5 Á 10 À3 A good agreement between the analytical results and the numerical approach is always observed. In particular, the NS boundary provides a nearly perfect prediction of the quasi-periodic regime arising and a very good estimate of the incommensurate frequency. The formulas proposed prove useful for design applications on real devices and structures. The closed form solution allows exploring in depth the resonance phenomena with low computational effort. Further development are currently undertaken in order to reframe the cases of 1:1 and 1:3 internal resonance within the same analyses, underlining the existence of different families of periodic orbits , deriving accurate predictive solutions for the NS boundary curve, in order to shed new light and unify analyses of the emergence of QP solutions in nonlinear oscillators with applications to MEMs dynamics.
Acknowledgements The authors thanks Yichang Shen for his help in checking the MS developments and Andrea Opreni for the development of the latest version of the continuation HB FEM code used in the full order simulation.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

General nonlinear coefficients in Eq. 1
This section shows how the formulas proposed in the paper can be used into the general case of uncorrelated a 12 and a 11 . We will mainly focus on the effect of a different sets of nonlinear coefficients on backbones and NS boundary. Considering the backbone curves, keeping the frequency detuning x 2 À 2x 1 fixed and varying the non-linear coefficients a 12 and a 11 ratios, the results are shown in Figs. 21 and 22.  Fig. 23 The figures show the NS boundary behaviour corresponding to x 1 ¼ 1, x 2 ¼ 2:005, l 1 ¼ l 2 ¼ 1:5 Á 10 À3 for different ratio a 12 =a 11 keeping a 12 ¼ 5 Á 10 À2 Fig. 21 assumes a 12 ¼ 1 Á 10 À2 and it considers different values of the ratio a 12 =a 11 . Figure 21 a) and b) refer to a detuning value r 1 ¼ 0, while Fig. 21c) and d) refer to r 1 ¼ 5 Á 10 À3 . The backbones of the second oscillator response is invariant with respect to the ratio a 12 =a 11 . This is reasonable because only a 12 appears in Eq. (16a). The q 1 backbone tends to grow with the ratio a 12 =a 11 preserving the same shape . The r 1 value introduces a curvature in the q 1 mode response and left unchanged the starting point of the q 2 backbone. Fig. 22 assumes a 11 ¼ 1 Á 10 À2 and considers different values of the ratio a 11 =a 12 . Fig. 22 a) and b) refer to r 1 ¼ 0, while Fig. 22c) and d) to r 1 ¼ 5 Á 10 À3 . The backbone of both oscillators grow with the ratio a 11 =a 12 preserving their shape. The r 1 value introduces a curvature in the q 1 mode response and shifts to higher amplitudes the starting point of the q 2 backbone.
We now inspect the NS boundary behaviour. We consider the system parameters x 1 ¼ 1, x 2 ¼ 2:005, l 1 ¼ l 2 ¼ 1:5 Á 10 À3 and we plot in Fig. 23 the NS boundaries corresponding to different ratio a 12 =a 11 , keeping a 12 ¼ 5 Á 10 À2 . The plots show that in the q 1 amplitude the NS boundary tends to enlarge with a lower ratio (higher a 11 ) while it is identical in the q 2 response. This is reasonable because only a 12 appears in Eq. (33). Considering the same system, in Fig. 24 the NS boundaries for different ratio a 11 =a 12 , keeping a 11 ¼ 5 Á 10 À2 are plotted. The plots show that in both the amplitudes the NS boundary tends to enlarge with lower ratio (higher a 12 ).

Conservative system
This section presents the details of the MS solution in Sect. 2.1 for the conservative case, i.e. with no forcing and damping. We insert the MS approximation q i ðtÞ ¼ q i0 ðT 0 ; T 1 Þ þ eq i1 ðT 0 ; T 1 Þ into Eq. (1) and split the resulting terms using e as sorting parameter. Consequently, we get a first order system: D 2 0 q 10 þ x 2 1 q 10 ¼ 0; D 2 0 q 20 þ x 2 2 q 20 ¼ 0; ð37Þ and a second order system: D 2 0 q 11 þ x 2 1 q 11 ¼ À2D 0 D 1 q 10 þ a 12 q 10 q 20 ; in which D k i denotes the order k derivative with respect to the time scale i ¼ 0; 1. Eq. (37) describes two linear uncoupled oscillators with motion given by: Inserting the solution of Eq. (37) in Eq. (38) we get an expression with A i as unknowns. Nullifying secular terms we get the conditions:  Fig. 24 The figures show the NS boundary behaviour corresponding to x 1 ¼ 1, x 2 ¼ 2:005, l 1 ¼ l 2 ¼ 1:5 Á 10 À3 for different ratio a 11 =a 12 keeping a 11 ¼ 5 Á 10 À2 a 12 A 2 A 1 e iT 1 r 1 þ 2ix 1 A 0 1 ¼ 0; ð40aÞ A 2 1 a 11 þ 2ix 2 A 0 2 e iT 1 r 1 ¼ 0: To solve Eq. (40) we need to expand A i in polar form as A i ðT 1 Þ ¼ 1 2 a i ðT 1 Þ exp ðih i ðT 1 ÞÞ; i ¼ 1; 2. Insert A i in Eq. (40) and imposing that the real and imaginary part vanish independently, we get Eq. (4).
As required by the stability analysis in Sect. 2.3 we finally compute the Jacobian matrix of system. (7), i.e.:

Forced and damped system
Starting from Eq. (1) we sort the non-linearities using a book-keeping parameter e. We assume a first order approximation, thus we rewrite Eq. (1) in the form: € q 2 þ x 2 2 q 2 ¼ e½À2 l 2 _ q 2 À a 11 q 2 1 : ð42bÞ The external forcing angular frequency is related to the eigenfrequency through Eq. (23). Proceeding as for the conservative case (see from Eq. (38) to Eq. (40)), we get the solvability condition: Inserting the polar form of A i ; i ¼ 1; 2 (see Appendix 2) in Eq. (43) and imposing that the real and imaginary parts vanish, we get the system of four equations in Eq. (24). Considering the A1-mode condition a 1 6 ¼ 0 and a 2 ¼ 0, Eq. (1) admit a nonzero solution for a 1 only if a 11 is zero, but this is a degenerate condition. Thus, the A1-mode is not allowed.
Considering the A2-mode condition a 2 6 ¼ 0 and a 1 ¼ 0, Eq. (1) implies a zero forcing value, and consequently nonzero solutions exist only if the damping term also is zero. This is a degenerate condition since it is the one predicted by the conservative case, thus the A2-mode in damped and forced conditions is impossible.
The coupled solutions are obtained from Eq. (24), as detailed in Sect. 3, and are given by Eq. (26).
In order to address the stability of the coupled solutions derived in Sect. 3, one needs the Jacobian matrix for Eq. (26), i.e.: with: J 11 ¼ a 12 a 2 sin c 2 ð Þ À 4 l 1 x 1 4x 1 ; ; J 14 ¼ a 12 a 1 a 2 cos c 2 ð Þ 4x 1 ; J 24 ¼ À a 12 a 2 sin c 2 ð Þ 4x 1 ; and d 3 , are relevant. This provides a linearised version of the NS minima problem and has the solution X ¼ ðx 2 þ x 1 Þ=3 independent of its amplitude.