# Normal form analysis of nonlinear oscillator equations with automated arbitrary order expansions

André de Figueiredo Stabile, Cyril Touzé, Alessandra Vizzaccaro

► To cite this version:

André de Figueiredo Stabile, Cyril Touzé, Alessandra Vizzaccaro. Normal form analysis of nonlinear oscillator equations with automated arbitrary order expansions. Journal of Theoretical, Computational and Applied Mechanics, 2025, 10.46298/jtcam.13234. hal-04505455v3

HAL Id: hal-04505455

https://hal.science/hal-04505455v3

Submitted on 19 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**Identifiers**

DOI 10.46298/jtcam.13234

HAL hal-04505455v3

**History**

Received Mar 15, 2024

Accepted Feb 11, 2025

Published May 19, 2025

**Associate Editor**

Stéphanie CHAILLAT

**Reviewers**

Ulrich J. RÖMER

Ghislain RAZE

Anonymous

**Open Review**

HAL hal-04996854

**Supplementary Material**

Dataset
10.5281/zenodo.15197276

Software

github.com/adfs-dev/MORFE_Symbolic

See addendum

**Licence**

CC BY 4.0

©The Authors

# Normal form analysis of nonlinear oscillator equations with automated arbitrary order expansions

André DE FIGUEIREDO STABILE<fnref n="1" />, Cyril TOUZÉ<fnref n="1" />, and Alessandra VIZZACCARO<fnref n="2" />

Arbitrary order expansions for the automatic reduction and solutions of nonlinear vibratory systems have been developed successfully within the realm of the direct parametrisation of invariant manifolds. Whereas the method has been used with high-order expansions and large dimensional systems, this article proposes to look at the same problem from the opposite point of view. By using low-dimensional systems, symbolic computations, analytical developments and numerical verifications, this contribution analyses the reduced dynamics appearing in cases where a single master mode is involved, reviewing typical scenarios in nonlinear vibrations: primary resonance, sub- and superharmonic resonances and parametric excitation. To achieve this task, the normal form style is preferentially used. A symbolic open-source package is also provided to generalise the presented results to other styles, higher orders, and different scenarios. It is shown how the low-order terms allow recovering the classical solutions given by perturbation methods, and how the automated expansions allow one to generalise the analysis to arbitrary orders. When analytical solutions are not tractable anymore, numerical solutions are employed to underline how converged solutions are at hand when the validity limit of the expansions is not reached. All the results presented in this paper can thus be used to better understand the nonlinear dynamical solutions occurring in nonlinear vibrations, as well as from a system identification perspective, since the normal form is the simplest dynamical system displaying a given resonance scenario.

**Keywords:** nonlinear oscillations, normal form, parametrisation method, geometric nonlinearity, nonlinear resonance, asymptotic expansion, symbolic calculation

## 1 Introduction

Since its introduction in an abstract framework in (Cabré et al. 2003a; Cabré et al. 2003b; Cabré et al. 2005; Haro et al. 2016), the parametrisation method for invariant manifolds has been extensively used in order to produce arbitrary order expansions for model order reduction of systems with smooth nonlinearities. It has been employed in the field of nonlinear vibrations in (Haller and Ponsioen 2016), allowing demonstration of the existence and uniqueness of spectral submanifolds (SSMs). Since then, it has then been successfully extended for reduced order modeling of large dimensional problems discretized by the finite element (FE) procedure, see e.g. (Jain and Haller 2022; Li et al. 2022; Vizzaccaro et al. 2022; Opreni et al. 2023; Martin et al. 2023; Vizzaccaro et al. 2024).

Another powerful result given by the method is also to unify the two main approaches that have been used in the past to compute nonlinear normal modes (NNMs). The centre manifold technique has been exploited in the works by Shaw and Pierre to derive reduced-order models based on invariant manifold theory (S. Shaw and Pierre 1991; S. Shaw and Pierre 1993; S. W. Shaw et al. 1999), whereas the normal form approach has been proposed in (Jezequel and C. Lamarque 1991; Touzé et al. 2004; Touzé and Amabili 2006) in order to derive similar results. The exact link between the two approaches has been only recently uncovered by using the parametrisation method for invariant manifolds, which highlights that the two techniques are in fact two different styles of parametrisation that can be used to solve the invariance equation.

<footnote n="1">¹ IMSIA, CNRS, EDF, ENSTA Paris, Institut Polytechnique de Paris, Palaiseau, France</footnote>

<footnote n="2">² College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter, UK</footnote>
Whereas previous works on the subject focus on direct applications to large dimensional FE problems to underline the impressive gains in computing time that can be expected from the application of the method, this contribution aims to give more insights into the results that can be awaited from a broad use of the method to derive accurate, high-order, analytical results obtained thanks to symbolic computations on low-dimensional systems. To that purpose, a symbolic package which relies on the previous developments of the MORFE project (MORFE stands for Model Order Reduction for Finite Element problems) is released with the present article in order to help the analyst in producing automated solutions up to arbitrary order. This paper illustrates how such outputs can be used to analyse the high-order reduced dynamics in different scenarios, and using *ad hoc* assumptions. In particular, it will be clearly shown how the first-order solutions are completely equivalent to solutions from standard perturbation methods. When the solutions are analytically tractable, they will be symbolically analysed to show how one can easily get more accurate approximations that take into account the next orders. Finally, when the analytical solutions are too lengthy, it will be illustrated how one can then use numerical solutions to reach convergence to a reference solution, once the validity limit of the local expansions is not exceeded. The purpose of all these developments is thus to show the continuity that exists between low-order perturbative techniques and high-order numerical solutions, showing to the analyst how the method can be used in an integrated manner, providing both analytical approximations for the design and system understanding phase, and accurate numerical solutions.

The development and release of symbolic softwares for computing high-order normal form solutions, centre manifold approximations or automated perturbative techniques, is not new and has been largely documented in the past, see e.g. (Elphick et al. 1987; A. Roberts 1997; Leung and Q. Zhang 1998; Yagasaki 1998; Yagasaki 1999; Yagasaki and Ichikawa 1999; Huseyin and W. Zhang 2000; W. Zhang et al. 2000; Leung and Q. Zhang 2003) for different proposals. However, we believe that embedding the normal form development into the more general framework of the parametrisation method, and proposing a symbolic software to automatically derive such high-order solutions, including also the graph style parametrisation, is a necessary development to give more physical insights into the powerful results obtained when treating large dimensional systems. Moreover, we show in the course of the article that the results analysed are general and can be used to better understand the reduced dynamics of large-scale problems, as well as to propose arbitrary order analytical solutions for nonlinear vibration problems. Moreover, thanks to the original treatment of the non-autonomous forcing term recently proposed in (Vizzaccaro et al. 2024), the analysis can be enlarged to a high level of forcing and can treat any classical resonance scenario occurring in nonlinear vibrations: primary resonance, subharmonic and superharmonic resonance, parametric resonance, etc.

Another purpose of the analytical developments presented in this contribution is to give a detailed explanation of the different variants of normal forms that have been used in the literature on nonlinear vibration, by proposing a unified presentation that also explains their advantages and drawbacks. Indeed, while complex normal form (CNF) is generally in use in the mathematical literature, see e.g. (Haragus and Iooss 2011; Iooss and Adelmeyer 1999; Jezequel and C. Lamarque 1991; Gabale and Sinha 2009; A. J. Roberts 2014; Haller and Ponsioen 2016; Waswa and Redkar 2020), a real normal form (RNF) has been introduced in (Neild and Wagg 2011; Neild et al. 2015) and analyzed in (Vizzaccaro et al. 2022; Opreni et al. 2023) in the context of the parametrisation method. Furthermore, another variant of the normal form has been used in (Touzé et al. 2004; Touzé and Amabili 2006; Touzé 2014) in order to keep real oscillator-like equations throughout the process. This third variant is called the *oscillator normal form* (ONF) and has been fruitful in order to make a direct link with the calculations of invariant manifolds using the centre manifold technique as proposed by Shaw and Pierre, see e.g. (Touzé et al. 2021).

The software used to automatically produce the outputs of the high-order expansions presented in this paper is called *MORFE_Symbolic*. It is based on the implementation of the parametrisation method for nonlinear vibrating systems used in the MORFE project (Opreni et al. 2022). A key feature is also the treatment of the forcing term and the time dependence of the invariant manifolds, when harmonic excitation is considered (for the case of stochastic forcing, the reader is referred to (Wang and A. Roberts 2012)). This automatic treatment proposed in (Vizzaccaro et al. 2024), allows very general results that are not limited to a first-order expansion
in the forcing amplitude as proposed in (Breunung and Haller 2018; Ponsioen et al. 2018; Ponsioen et al. 2020; Jain and Haller 2022; Opreni et al. 2023). In particular, the nonlinear dependence upon the forcing amplitude can be taken into account with the method proposed in (Vizzaccaro et al. 2024), such that all types of resonance scenarios including the superharmonic can be analysed. The code MORFE_Symbolic also allows one to start from a differential algebraic equation (DAE), extending again the range of dynamical scenarios that can be analysed. This feature will be illustrated here in the case of the parametric excitation. A repository that contains different versions of the symbolic code is published in connection with this manuscript, where both Julia and Mathematica versions are released. The Julia version has the advantage of not using any proprietary software. However, today's symbolic computation capacities of the Julia package are not as efficient as those of Mathematica. Because of that, the more involved examples studied in this paper were treated with the Mathematica version, or using a mix of both codes, due to processing time.

As a final introductory remark, most of the results presented in the paper aim to analyse reduced-order dynamics by systematically understanding the effects of the terms produced by asymptotic expansions. The code MORFE_Symbolic, to which a short introduction is given in Appendix A, can also be used in a very general manner to produce many more different cases than the ones analysed here. We think that such a tool is a key requirement for the reinterpretation of the grammar of nonlinear oscillations. Another interesting feature of the results produced by the normal form approach is the creation of a dictionary of the simplest dynamical systems that contain a typical feature in nonlinear vibration. Such a dictionary can be very efficiently used in the realm of system identification with data-driven techniques. Thanks to the automated symbolic calculation, the method can be used to derive the dynamical systems that will produce the feature of e.g. a superharmonic resonance or a 1 : 2 internal resonance with the minimal number of monomials.

The outline of the paper is as follows. Section 2 is a pedagogical introduction to the different variants of normal form styles, which underlines some key features of the different techniques. Section 3 contains most of the analysis by focusing on cases where the reduced dynamics contains a single master mode. Typical scenarios in nonlinear vibrations are analysed: primary resonance, super- and sub-harmonic resonance and parametric excitation. Section 4 extends some of these results by accurately quantifying the effect that a slave mode can have on the reduced dynamics with a single master mode. In order to restrict the length of the paper, this section is deliberately shorter than the previous one and closes with an illustrative example.

## 2 Variants of normal form styles

This pedagogical introductory section aims to explain the different possible styles of normal forms that have been proposed in the literature for nonlinear vibratory systems. To that purpose and to make the presentation as simple as possible, the case of the unforced and undamped Duffing oscillator equation is considered, as it is sufficient to understand the different variants. The findings extend naturally to coupled nonlinear oscillators, which are studied in the next sections.

The starting point is the conservative unforced cubic Duffing oscillator

$$ \ddot{u} + \omega^2 u + hu^3 = 0. \qquad (1) $$

The three main variants that will be emphasised in this presentation are: the complex normal form (CNF), the real normal form (RNF) and the oscillator normal form (ONF). The complex normal form is recommended by mathematical textbooks (Iooss and Adelmeyer 1999; A. J. Roberts 2014), and has been used for example in (Jezequel and C. Lamarque 1991; Haller and Ponsioen 2016; Waswa and Redkar 2020). A pedagogical introduction with application to Duffing and Mathieu equations is also provided in (A. J. Roberts 2014). The starting point is to write Equation (1) at first-order, with a diagonalised linear part that contains the eigenvalues {$\pm i\omega$}. By doing so, complexification is enforced and the link to an oscillator-like equation is lost. This is the main reason why alternative procedures have been proposed in the vibration literature (namely RNF and ONF), in order to ease the realification and keep the link with oscillator equations.
To rewrite the system with a first-order diagonal linear part, the following linear change of coordinate can be applied

$$ \begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ i\omega & -i\omega \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix}, \qquad (2) $$

where $v = \dot{u}$ has been introduced. Note that other linear transforms with different normalisations can be applied here. Equation (1) then becomes

$$ \dot{y}_1 = i\omega y_1 + i \frac{h}{2\omega} (y_1^3 + 3y_1^2 y_2 + 3y_1 y_2^2 + y_2^3), \qquad (3a) $$

$$ \dot{y}_2 = -i\omega y_2 - i \frac{h}{2\omega} (y_1^3 + 3y_1^2 y_2 + 3y_1 y_2^2 + y_2^3). \qquad (3b) $$

In particular, one can see that, whereas only one monomial $u^3$ is present in Equation (1), eight monomials are in Equation (3). Note also that Equation (3b) is the complex conjugate of Equation (3a). This is the direct consequence of the fact that the initial problem, Equation (1) is second-order in time such that, when rewriting it as a first-order dynamical system, one of the two equations, namely $v = \dot{u}$, is tautological.

The normal form procedure can be unfolded on Equation (3). It starts with introducing a nonlinear change of coordinates between modal coordinates ($y_1, y_2$) and normal coordinates ($z_1, z_2$) as

$$ y_1 = z_1 + a_{11}z_1^3 + a_{12}z_1^2z_2 + a_{13}z_1z_2^2 + a_{14}z_2^3 \qquad (4a) $$

$$ y_2 = z_2 + a_{21}z_1^3 + a_{22}z_1^2z_2 + a_{23}z_1z_2^2 + a_{24}z_2^3. \qquad (4b) $$

In Equation (4), no quadratic monomials are present because the original system only contains cubic terms. The normal dynamics is also introduced as

$$ \dot{z}_1 = i\omega z_1 + f_{11}z_1^3 + f_{12}z_1^2z_2 + f_{13}z_1z_2^2 + f_{14}z_2^3, \qquad (5a) $$

$$ \dot{z}_2 = -i\omega z_2 + f_{21}z_1^3 + f_{22}z_1^2z_2 + f_{23}z_1z_2^2 + f_{24}z_2^3, \qquad (5b) $$

where the linear part is preserved because the nonlinear change of coordinates Equation (4) is identity-tangent.

The sixteen unknown coefficients {$a_{ij}$} and {$f_{ij}$}, $i \in \{1, 2\}$ and $j \in \{1, 2, 3, 4\}$, are found by plugging Equations (4) to (5) into Equation (3), and comparing coefficients for each monomial of the normal coordinates. This leads to the four following equations, related to Equation (3a), which are known in the literature as order-three homological equations (note that homological equations are then automatically derived for each higher order):

$$ \text{For } z_1^3: \quad (i\omega + i\omega + i\omega) a_{11} + f_{11} = i\omega a_{11} + i \frac{h}{2\omega}, \qquad (6a) $$

$$ \text{For } z_1^2 z_2: \quad (i\omega + i\omega - i\omega) a_{12} + f_{12} = i\omega a_{12} + i \frac{3h}{2\omega}, \qquad (6b) $$

$$ \text{For } z_1 z_2^2: \quad (i\omega - i\omega - i\omega) a_{13} + f_{13} = i\omega a_{13} + i \frac{3h}{2\omega}, \qquad (6c) $$

$$ \text{For } z_2^2: \quad (-i\omega - i\omega - i\omega) a_{14} + f_{14} = i\omega a_{14} + i \frac{h}{2\omega}, \qquad (6d) $$

while another four, not reported here for the sake of brevity, are derived from Equation (3b).

Equation (6) is an underdetermined system of four equations for eight unknowns, admitting an infinity of possible solutions. Equations (6a) to (6d) can be solved easily for $a_{11}$, $a_{13}$ and $a_{14}$, by imposing $f_{11} = f_{13} = f_{14} = 0$. This choice is the classical one to derive the normal form of the system, with the idea of simplifying as much as possible the normal dynamics Equation (5). Once the coefficients of the change of coordinates are determined, if Equation (4) is replaced into Equation (3), terms of order higher than 3 will appear in the expressions, such that another change of coordinates, similar to Equation (4) but with higher order terms, could be made,
generalizing the procedure up to arbitrary order. At each time, new homological equations are derived in order to solve for the coefficients of the system.

A *trivial* resonance relationship occurs in Equation (6b) since $a_{12}$ vanishes, such that one cannot select $f_{12} = 0$. This is a direct consequence of the eigenspectrum composed of a purely complex conjugate pair {$\pm i\omega$}. The three variants of normal form styles discussed herein depart on the choices made in order to solve Equation (6), which can also be viewed as a more or less stringent interpretation of the resonance relationship. This is detailed now for each of the three normal form styles.

## 2.1 The complex normal form style

The complex normal form (CNF) style is the classical treatment proposed in mathematical textbooks to deal with purely imaginary complex eigenspectrum, see e.g. (Iooss 1988; Haragus and Iooss 2011; Wiggins 2003; Jezequel and C. Lamarque 1991). Referring to the simple case of the Duffing equation, it amounts to cancel the three non-resonant monomials $z_1^3, z_1z_2^2$ and $z_2^3$ in Equations (6a) to (6d) with the choice $f_{11} = f_{13} = f_{14} = 0$. For the only resonant monomial $z_1^2z_2$, then one imposes $a_{12} = 0$ and $f_{12} = 3ih/(2\omega)$. This follows from a strict interpretation of the resonance relationship stemming from the homological equations at each order. The generalisation of the resonance relationship to arbitrary order can be found in many classical books, see e.g. (Poincaré 1892; Iooss 1988; Guckenheimer and Holmes 1983; Manneville 1990; Murdock 2003), and reads, for a dynamical systems of dimension $n$ with eigenvalues {$\lambda_1, \dots, \lambda_n$}, and for a monomial of order $p$:

$$ \lambda_k = \sum_{i=1}^{n} m_i \lambda_i, \quad \text{with } m_i \ge 0 \quad \text{and} \quad \sum_{i=1}^{n} m_i = p. \qquad (7) $$

For conservative mechanical systems with an eigenspectrum composed of pairs of purely imaginary numbers, trivial resonance relationships appear at each odd order in the normal form computation. Focusing on the simple case of the Duffing equation with cubic nonlinearity, an interesting feature of the CNF is that, for each odd order, only one resonant monomial stays in the normal form, as a consequence of these trivial resonances. Symbolic calculation of the normal form up to arbitrary order is possible thanks to the code `MORFE_Symbolic`, which has been used to develop the calculations shown in this paper. As an illustration, we give below the CNF for the Duffing equation up to order 11:

$$ \dot{z}_1 = i\omega z_1 + i \frac{3h}{2\omega} z_1^2 z_2 - i \frac{51h^2}{2^4\omega^3} z_1^3 z_2^2 + i \frac{1419h^3}{2^7\omega^5} z_1^4 z_2^3 - i \frac{47505h^4}{2^{10}\omega^7} z_1^5 z_2^4 + i \frac{438825h^5}{2^{11}\omega^9} z_1^6 z_2^5, \quad (8a) $$

$$ \dot{z}_2 = -i\omega z_2 - i \frac{3h}{2\omega} z_1 z_2^2 + i \frac{51h^2}{2^4\omega^3} z_1^2 z_2^3 - i \frac{1419h^3}{2^7\omega^5} z_1^3 z_2^4 + i \frac{47505h^4}{2^{10}\omega^7} z_1^4 z_2^5 - i \frac{438825h^5}{2^{11}\omega^9} z_1^6 z_2^5. \quad (8b) $$

As announced, one can observe that only one resonant monomial of the form $z_1^{p+1}z_2^p$ stays in the normal dynamics for each odd order $2p + 1$. Equation (8b) is the complex conjugate of Equation (8a), and one can also observe the change of sign in each successive odd order coefficients. For the sake of completeness, the nonlinear change of coordinates is reported in Appendix B.

An essential property of the CNF solution for conservative systems is that an analytic backbone curve is easily computed for any order. This property is known and has been, for example, already used in (Breunung and Haller 2018). Let us recall why this property is essentially linked to the CNF style. In order to derive physical characteristics of the original system, or even just to compute numerical continuation solution on the reduced dynamics, a realification procedure needs to be applied to the complex normal form as given e.g. in Equation (8), see for example (Haro et al. 2016; Vizzaccaro et al. 2022; Opreni et al. 2023) for general discussions. Different realification procedures can be used, choosing for instance cartesian or polar coordinates. A key feature of the CNF is to provide very simple expressions when realification with polar coordinates is employed. Let us introduce the polar form of the normal coordinates as

$$ z_1 = 1/2\rho e^{i\alpha}, \qquad (9a) $$
$$z_2 = 1/2\rho e^{-i\alpha}, \qquad (9b)$$

where it should be noted that $\rho$ and $\alpha$ are functions of time, even though this dependence is not made explicit to lighten the notation. Since only one resonant monomial ($z_1^{p+1}z_2^p$ for the first equation and $z_1^p z_2^{p+1}$ for the second equation), is present for each odd order, a simple calculation shows that, whatever the order of the expansion used for the normal form, the dynamical equation for the amplitude $\rho$ simply reads $\dot{\rho} = 0$. This is a direct consequence of the existence of the Lyapunov subcenter manifold (LSM), that is densely filled with a family of periodic orbits that are parametrised by their amplitude $\rho$. The equation for the phase then contains all the important dynamical terms, and reads, for example, up to order 11:

$$\dot{\alpha} = \omega + \frac{3h}{2^3\omega}\rho^2 - \frac{51h^2}{2^8\omega^3}\rho^4 + \frac{1419h^3}{2^{13}\omega^5}\rho^6 - \frac{47505h^4}{2^{18}\omega^7}\rho^8 + \frac{438825h^5}{2^{21}\omega^9}\rho^{10}. \qquad (10)$$

Note in this case that a simplification by $\rho$ is made in this last equation; this explains why an odd order $2p + 1$ in the normal form gives rise to a power $\rho^{2p}$ in the resulting equation for the phase. Since $\rho$ is constant then this equation can be directly integrated and reveals the frequency-amplitude relationship of the conservative problem, also known as the backbone curve. With the help of the symbolic processor MORFE_Symbolic, high-order backbone curves can be easily derived; an example is given below up to order 11:

$$\omega_{\text{NL}} = \omega \left( 1 + \frac{3h}{2^3\omega^2}\rho^2 - \frac{51h^2}{2^8\omega^4}\rho^4 + \frac{1419h^3}{2^{13}\omega^6}\rho^6 - \frac{47505h^4}{2^{18}\omega^8}\rho^8 + \frac{438825h^5}{2^{21}\omega^{10}}\rho^{10} \right). \qquad (11)$$

Note that this property is particularly important and meaningful since it applies to the solutions derived by the parametrisation method using CNF. Consequently, even when dealing with FE problems that involve a very large number of degrees-of-freedom (dofs), an analytical backbone curve can be derived when reducing the problem to a single NNM with CNF, meaning that no extra calculation (like a continuation procedure) is needed to obtain the amplitude-frequency relationship for the conservative problem (Breunung and Haller 2018).

It is important to highlight, however, that the above backbone curve is written as a function of the amplitude $\rho$ of the normal variables $z_1$ and $z_2$, and not that of the maximum value of the physical displacement, which we denote $u_{\text{max}}$ in what follows. In order to obtain the physical displacement amplitude, the normal coordinates' polar representation can be inserted into the nonlinear mapping displacement equation. Substitution of Equation (9) into Equation (B.1a), limited here at order 3 for the sake of brevity, yields

$$u = \rho \cos \alpha + \frac{h}{32\omega^2}\rho^3 \cos 3\alpha - \frac{3h}{16\omega^2}\rho^3 \cos \alpha. \qquad (12)$$

Since all cosines are in phase, the maximum displacement amplitude is simply given by the sum of their coefficients, obtained for $\alpha = 0$. This procedure is easily extendable to any order and is automated in MORFE_Symbolic. Up to order 11, the relationship between the maximum physical displacement amplitude $u_{\text{max}}$ and the amplitude of the normal coordinate $\rho$ is

$$u_{\text{max}} = \rho \left( 1 - \frac{5h}{2^5\omega^2}\rho^2 + \frac{25h^2}{2^8\omega^4}\rho^4 - \frac{2781h^3}{2^{15}\omega^6}\rho^6 + \frac{90493h^4}{2^{20}\omega^8}\rho^8 - \frac{3234957h^5}{2^{25}\omega^{10}}\rho^{10} \right). \qquad (13)$$

With Equations (11) to (13) it is possible to find coordinate pairs $(\omega_{\text{NL}}, u_{\text{max}})$ that yield the backbone of the systems parametrized by $\rho$. Despite the distinction between $u_{\text{max}}$ and $\rho$, it is interesting to note, nevertheless, that these two quantities coincide up to the second order, which implies that Equation (11) can be used as an approximation for the backbone for small displacement amplitudes.

As an illustration, let us show, for this simple example of the conservative cubic Duffing oscillator, the convergence of the backbone curve for different orders $o$ of normal form expansion. Figure 1 shows backbone curves obtained with CNF in the way described above up to order $o = 25$. They are compared with the backbone in terms of $u_{\text{max}}$ found using an exact solution based on Jacobi elliptic functions, available for example in (Salas and Castillo H. 2014).
<figure><img src="image_2.png" /><figcaption>Figure 1 Backbone curves for the Duffing oscillator up to order o = 25 with a complex normal form style and parameter values set as ω = 1, h = 1. The backbone curves are compared with an analytical solution based on Jacobi elliptic functions (Salas and Castillo H. 2014). (a) Comparison in terms of the maximum displacement u<sub>max</sub>. (b) Comparison as a function of the amplitude ρ of the normal variable.</figcaption></figure>

From Figure 1(a) it can be seen that there exists a maximum validity range for the approxima-
tions given by the normal form expansion, corresponding in this case to a displacement amplitude
of approximately u<sub>max</sub> = 0.8. This is in line with the fact that the normal form relies on a local
theory, and this upper bound for validity limit has been explored for example in (C.-H. Lamarque
et al. 2012). Figure 1(b) shows the same backbone as a function of the amplitude ρ of the normal
variable. Since Equation (13) gives u<sub>max</sub> as a function of ρ, and since the analytical solution
with Jacobi elliptic functions is given for the amplitude u<sub>max</sub>, Equation (13) needs to be inverted
to plot the reference solution. This has been done in Figure 1(b) together with a first-order
approximation u<sub>max</sub> = ρ which is frequently used in perturbation methods. This highlights
the difference between u<sub>max</sub> and ρ as amplitudes are increasing. In this case, the first-order
approximation u<sub>max</sub> = ρ can be used for amplitudes up to about ρ = 0.5, but higher orders cannot
be neglected thereafter. It should be noted that Equation (13) is no longer valid once a certain level
of amplitude is reached, as can be seen by the sudden deviation of the black curve on Figure 1(b)
from where it is graphically expected to go. This deviation occurs at approximately ρ = 0.9,
which corresponds to u<sub>max</sub> = 0.85, and is clearly beyond the validity limit of the approximation.

As a last point of discussion, the effect of a viscous damping term in the Duffing equation is
investigated. The damping is appended to Equation (1) with a modal damping factor ξ, and gives

$$
\ddot{u} + 2\xi\omega\dot{u} + \omega^2 u + hu^3 = 0. \tag{14}
$$

The damping affects the eigenvalues, which now become

$$
\lambda_{1,2} = -\xi\omega \pm i\omega\sqrt{1 - \xi^2} = \omega(-\xi \pm i\delta), \quad (15)
$$

where $\delta = \sqrt{1 - \xi^2}$ has been introduced in order to ease the analytical expressions. A small damping assumption is made such that (i) underdamped oscillations are considered ($\xi < 1$), and (ii) the resonance check in the normal form procedure is done on the imaginary (oscillatory) parts only, as proposed in (Touzé and Amabili 2006; Haller and Ponsioen 2016; Vizzaccaro et al. 2022) for nonlinear vibratory systems. Due to the increase in length we consider the CNF up to order 7
only:

$$ \dot{z}_1 = (i\delta\omega - \xi\omega)z_1 + i\frac{3h}{2\delta\omega}z_1^2z_2 - i\frac{3h^2(17 - 14\xi^2 + 4i\xi\delta)}{8\omega^3\delta^2(-2i\delta + \xi)}z_1^3z_2^2 \\ - i\frac{3h^3(378\xi^6 - 48\xi^4 - 1420\xi^2 - 531i\delta\xi + 378i\delta\xi^5 - 207i\delta\xi^3 + 946)}{64\omega^5\delta^3(9\xi^4 - 9\xi^2 - 4)}z_1^4z_2^3, \quad (16) $$

the second equation for $z_2$ being the complex conjugate of this one, i.e. $\dot{z}_2 = \dot{z}_1$. The linear term is due to the change in the eigenvalue. The same resonant monomials are present in the normal form since the resonance relationships are verified based on the assumption of small damping. Only the coefficients are slightly modified by the damping ratio $\xi$. Proceeding with realification as in the undamped case, it is possible to find the following two equations, by considering the real and imaginary parts of Equation (16) once $z_1$ and $z_2$ are substituted by their polar representation:

$$ \dot{\rho} = -\xi\omega\rho + \frac{3^2\xi h^2(2\xi^2 - 3)}{2^7\delta^2\omega^3(1 + 3\delta^2)}\rho^5 - \frac{3^3\xi h^3(42\xi^4 - 23\xi^2 - 59)}{2^{12}\delta^2\omega^5(1 + 3\delta^2)(1 + 3\xi^2)}\rho^7 \quad (17a) $$

$$ \dot{\alpha} = \delta\omega + \frac{3h}{2^3\delta\omega}\rho^2 + \frac{3h^2(12\xi^2 - 17)}{2^6\delta\omega^3(1 + 3\delta^2)}\rho^4 + \frac{3h^3(189\xi^6 - 24\xi^4 - 710\xi^2 + 473)}{2^{11}\delta\omega^5(1 + 3\delta^2)(1 + 3\xi^2)}\rho^6. \quad (17b) $$

Equation (17a) gives the decay rate of the oscillation amplitude. Note that up to the first order, the equation can be integrated to give the usual exponential decay predicted by the linear theory. A nonlinear damping ratio $\xi_{NL}$ can be defined by dividing Equation (17a) by $-\rho\omega$. A nonlinear instantaneous frequency $\omega_{NL}$ can also be derived from Equation (17b) through $\dot{\alpha} = \omega_{NL}$, which yields the two equations

$$ \xi_{NL} = \xi \left( 1 - \frac{3^2 h^2 (2\xi^2 - 3)}{2^7 \delta^2 \omega^4 (1 + 3\delta^2)} \rho^4 - \frac{3^3 h^3 (42\xi^4 - 23\xi^2 - 59)}{2^{12} \delta^2 \omega^6 (1 + 3\delta^2) (1 + 3\xi^2)} \rho^6 \right) \quad (18a) $$

$$ \omega_{NL} = \delta\omega \left( 1 + \frac{3h}{2^3\delta^2\omega^2}\rho^2 + \frac{3h^2(12\xi^2 - 17)}{2^6\delta^2\omega^4(1 + 3\delta^2)}\rho^4 + \frac{3h^3(189\xi^6 - 24\xi^4 - 710\xi^2 + 473)}{2^{11}\delta^2\omega^6(1 + 3\delta^2)(1 + 3\xi^2)}\rho^6 \right). \quad (18b) $$

Equation (18b) can be viewed as a damped backbone curve that smoothly perturbs from the undamped case when $\xi \neq 0$ (Touzé and Amabili 2006; Llave and Kogelbauer 2019). It gives an analytical formula for taking the damping into account in the backbone curve and quantitatively compares conservative and damped cases. Equation (18a) underlines that the linear viscous damping creates a nonlinear decay rate. However, its effect is only apparent at large amplitudes since the nonlinearities are of order $\rho^4$ or higher.

In order to better understand the evolution of the above quantities, the terms inside the brackets in Equations (18a) and (18b) are plotted as a function of $\rho$ for $\omega = 1$ and $h = 1$ in Figure 2(a) and Figure 2(b). From the figures, it is possible to notice that the inclusion of damping does not significantly alter the shape of the backbones up to $\rho = 0.8$, which corresponds to the validity limit of the asymptotic expansion.

It is also worth mentioning that, in contrast to the undamped case, obtaining the maximum physical displacement $u_{max}$ as a function of amplitude $\rho$ is not feasible analytically, since the coefficients of the nonlinear mappings are now complex. Both cosines and sines appear in the equation for $u(\rho)$ (analogous to Equation (12)), hence there is no simple analytical solution for $u_{max}$.

As a summary, the CNF shares a number of advantageous characteristics. It has a lot of symmetries, leads to the most parsimonious representation of the normal dynamics, and analytical backbone curves are directly attainable with polar realification. The only drawback is that the normal dynamics is expressed with complex coordinates. In order to obtain normal form calculations that stay in a real formulation, different variants have thus been proposed in the past.

## 2.2 The real normal form

The real normal form (RNF) has been first introduced by Neild and Wagg (2011); Neild et al. (2015); Liu and Wagg (2019); Nasir et al. (2021); Wagg (2022), and reformulated in the context of the parametrisation method in (Vizzaccaro et al. 2022; Opreni et al. 2023). Note that the denomination
<figure><img src="image_3.png" /><figcaption>Figure 2 Duffing oscillator with increasing values of the damping ratio $\xi$: (a) damped backbone curve and (b) nonlinear damping ratio. The results are obtained using a complex normal form style, for the parameter values $\omega = 1$ and $h = 1$.</figcaption></figure>

RNF is not used by Neild and Wagg who called the method DNF for direct normal form, see e.g. (Elliott et al. 2018). In their case, the adjective direct was used to mean that the normal form technique is directly applied to second-order problems in time that are under study for mechanical vibration. Also in their approach, the developments were slightly different from the one presented herein since the target was to obtain time-domain approximations of the solutions. On the other hand, the method has been named RNF in (Vizzaccaro et al. 2022) in order to keep the adjective direct to refer to calculations that affect the physical, and not the modal coordinates.

The main idea of the RNF consists of keeping two monomials as resonant in Equation (3a). In addition to the term $z_1^2 z_2$, which is trivially resonant and kept in the CNF, the monomial $z_1 z_2^2$ is also defined as a resonant one. This implies that in Equation (6) the choice to zero coefficient $f_{12}$ is not made, but rather coefficient $a_{12}$ is set to zero. In such way, $f_{12} = i3h/(2\omega)$, changing the final expressions for Equations (4) to (5). The origin of this choice lies in the fact that the goal of the calculation proposed in (Neild and Wagg 2011; Neild et al. 2015) was to apply the normal form computation to second-order problem, which gives additional constraint to the developments. Since the second derivative of the mapping with respect to time needs to be computed, the homological equations need to be written with these terms, and squares of the eigenfrequencies appear naturally in the resonance relationships. This can also be seen in the homological equations derived in (Vizzaccaro et al. 2022) that have been rewritten only for the displacement mapping, by using the relationship that exists at any order between the displacement and the velocity mappings. In this case, the ill-conditioning of the term in front of the mapping, which is due to resonance relationships, appears with squared values. Rewriting the resonance relationships Equation (7) with squares leads to

$$ (\lambda_k)^2 = \left( \sum_{i=1}^{n} m_i \lambda_i \right)^2, \quad \text{with } m_i \ge 0 \quad \text{and} \quad \sum_{i=1}^{n} m_i = p, \qquad (19) $$

which is indeed the resonance relationship used to derive the RNF. For the sake of illustration, the RNF of the Duffing equation is here given up to order 7, it reads

$$ \dot{z}_1 = i\omega z_1 + i \frac{3h}{2\omega} z_1^2 z_2 + i \frac{3h}{2\omega} z_1 z_2^2 - i \frac{15h^2}{2^4\omega^3} z_1^3 z_2^2 - i \frac{3h^2}{2^3\omega^3} z_1^2 z_2^3 + i \frac{267h^3}{2^7\omega^5} z_1^4 z_2^3 - i \frac{3h^3}{2^7\omega^5} z_1^3 z_2^4, \qquad (20a) $$

$$ \dot{z}_2 = -i\omega z_2 - i \frac{3h}{2\omega} z_1^2 z_2 - i \frac{3h}{2\omega} z_1 z_2^2 + i \frac{3h^2}{2^3\omega^3} z_1^3 z_2^2 + i \frac{15h^2}{2^4\omega^3} z_1^2 z_2^3 + i \frac{3h^3}{2^7\omega^5} z_1^4 z_2^3 - i \frac{267h^3}{2^7\omega^5} z_1^3 z_2^4, \qquad (20b) $$
while the nonlinear mapping is given in Appendix C. As a consequence of the choice retained for fulfilling the resonance relationship, the same monomials now appear on the two lines of the normal dynamics. Interestingly, the cubic order terms share the same coefficients. This property is useful in order to retrieve an oscillator equation when coming back to real coordinates using a cartesian representation. We define the cartesian real coordinates ($a_1, a_2$) as

$$a_1 = z_1 + z_2, \qquad (21a)$$

$$a_2 = \frac{z_1 - z_2}{i}. \qquad (21b)$$

By stopping the RNF developments in Equation (20) at order three, one can see that the dynamics for $a_1$ is simple and reads $\dot{a}_1 = -\omega a_2$. Consequently, an oscillator equation can be written for the cartesian coordinates ($a_1, a_2$) as

$$\ddot{a}_1 + \omega^2 a_1 + \frac{3h}{4}(a_1^3 + a_1 \frac{\dot{a}_1^2}{\omega^2}) = 0. \qquad (22)$$

Whereas the CNF was conveniently realified with polar coordinates, the RNF is better suited for realification using cartesian coordinates, thanks to the symmetry properties appearing in the coefficients of Equation (20). Unfortunately the nice property $\dot{a}_1 = -\omega a_2$ is completely linked to the fact that cubic coefficients are all equal. From order 5, this property is lost, but realification is still easy to manage since the two different coefficients sum up and hence can be factorized. Equations (20) to (21) yield

$$\dot{a}_1 = -\omega a_2 + \frac{9h^2}{256\omega^3} a_1^4 a_2 + \frac{9h^2}{128\omega^3} a_1^2 a_2^3 + \frac{9h^2}{256\omega^3} a_2^5 \\ - \frac{135h^3}{4096\omega^5} a_1^6 a_2 - \frac{405h^3}{4096\omega^5} a_1^4 a_2^3 - \frac{405h^3}{4096\omega^5} a_1^2 a_2^5 - \frac{135h^3}{4096\omega^5} a_2^7, \qquad (23a)$$

$$\dot{a}_2 = \omega a_1 + \frac{3h}{4\omega} a_1^3 + \frac{3h}{4\omega} a_1 a_2^2 - \frac{21h^2}{256\omega^3} a_1^5 - \frac{21h^2}{128\omega^3} a_1^3 a_2^2 - \frac{21h^2}{256\omega^3} a_1 a_2^4 \\ + \frac{33h^3}{1024\omega^5} a_1^7 + \frac{99h^3}{1024\omega^5} a_1^5 a_2^2 + \frac{99h^3}{1024\omega^5} a_1^3 a_2^4 + \frac{33h^3}{1024\omega^5} a_1 a_2^6. \qquad (23b)$$

From these two equations, one can try to recover an oscillator-like equation. Deriving the first equation with respect to time, replacing $\dot{a}_2$ by the second expression, and stopping the developments at order 5, one obtains:

$$\ddot{a}_1 + \omega^2 a_1 + \frac{3h}{4} (a_1^3 + a_1 a_2^2) - \frac{15h^2}{128\omega^2} (a_1^5 + 2a_1^3 a_2^2 + a_1 a_2^4) = 0. \qquad (24)$$

This underlines that, for order 5 and higher, a closed-form expression involving only $a_1$ is not possible anymore such that two equations need to be kept with both $a_1$ and $a_2$. Equation (23) or Equation (24) can then be solved with a perturbative approach or a numerical continuation method to derive the backbone curve.

All these results show that the RNF can be advantageously used with Cartesian coordinates for realification. It can be easily automated since the choice of resonant monomials derives from a broader interpretation of the resonance relationship. Arbitrary order solutions are at hand and can be computed thanks to MORFE_Symbolic. However, the resulting equations are not oscillator-like from order 5. Besides, it seems that no simple and exact solution allows giving an analytical backbone curve at arbitrary order, as was the case for the CNF. To draw out such a solution a few more assumptions need to be made, following for example the approximations used in (Neild and Wagg 2011; Neild et al. 2015). This is illustrated in Appendix C.

## 2.3 The oscillator normal form

The oscillator normal form (ONF) has been first introduced in (Touzé 2003; Touzé et al. 2004; Touzé and Amabili 2006; Touzé 2014), with the main idea of keeping oscillator equations without using any complex formulation. To that purpose, the linear part is not made diagonal with complex
entries, but stays under its anti-diagonal formulation, see (Touzé et al. 2004). To better understand how the ONF can be interpreted from the previous example where complexification is used, the key point is to understand that all complex monomials appearing due to the complexification of a real one, need to be kept in the analysis, in order to make possible the come back to oscillator-like equations. For the Duffing equation, the real monomial $u^3$ in Equation (1) gives rise to 8 monomials for $y_1$ and $y_2$ as shown in Equation (3). Since all these are mandatorily needed to be able to reconstruct $u^3$, it means that no terms in Equation (6) will be cancelled. The choice $a_{11} = a_{12} = a_{13} = a_{14} = 0$ is selected. As a consequence, the Duffing Equation (1) is under its oscillator normal form. In ONF, the trivially resonant monomial is $u^3$ and cannot be cancelled, but many other terms, which are not linked to trivial resonances, can be cancelled in the process. In particular, all quadratic terms are not resonant and can be eliminated thanks to a nonlinear change of coordinate, see (Touzé et al. 2004) for general discussions and (Touzé 2014) for examples and classification of nonlinear terms thanks to this interpretation of the resonance relationship.

One of the main advantages of the ONF is thus to keep oscillator-like equations throughout the process. The nonlinear change of coordinate is given between two real coordinates that are homogeneous to a displacement and a velocity, whereas this interpretation is lost when using complex formulations. This choice came along with other advantages. For example, only the ONF allows drawing out a term-by-term comparison of the NNM calculation using either the centre manifold technique as proposed by Shaw and Pierre, or the normal form approach, see e.g. (Touzé et al. 2004; Touzé et al. 2021) for such discussions. Thanks to the ONF, a direct comparison with the quadratic manifold method with modal derivatives has also been made possible. Indeed, since real coordinates are used in the two approaches, term-by-term comparisons are at hand, which allows understanding that the quadratic manifold is a simplification of the general formula given by ONF, see (Vizzaccaro, Salles, et al. 2021; Shen et al. 2021; Touzé et al. 2021) for more details. A final advantage of the ONF is that it can be rewritten from physical coordinates, which allows deriving a non-intrusive version of the reduction technique using the normal form, which has been named DNF for direct normal form, see (Vizzaccaro, Shen, et al. 2021; Opreni et al. 2021).

However, numerous drawbacks are linked to this formulation. First, it is difficult to translate the choice on the resonant monomial as a broader algebraic interpretation of the resonance relationships, as it has been possible for the RNF with Equation (19). As a consequence, it appears very difficult (and maybe not possible) to generalise the ONF to arbitrary order and automate its computation. As an illustration of these complications, one can refer to (Shami et al. 2022) to see how the ONF with resonant quadratic terms can be computed up to order three. As a consequence, the ONF will not be much commented on in the rest of this paper, and is not included in MORFE_Symbolic.

# 3 High-order solutions for single-degree-of-freedom systems with forcing and damping

This section extends the analysis using normal form expansions to illustrate how this formalism allows an understanding of the main features of nonlinear oscillations. Primary and secondary resonances will be analysed. The first case under study is a Duffing oscillator with quadratic and cubic nonlinearities, derived to explain how the even order nonlinear terms can be cancelled and how the high-order terms on the backbone can be analysed in terms of hardening/softening behaviour.

## 3.1 Duffing oscillator with quadratic and cubic nonlinearities

In this section, the results of the previous section are generalised to a Duffing equation with quadratic and cubic nonlinearities:

$$ \ddot{u} + \omega^2 u + gu^2 + hu^3 = 0. \qquad (25) $$

Using the automated symbolic development provided by MORFE_Symbolic, the complex normal form (CNF) can be written up to arbitrary order. The first step consists in diagonalising the linear part using Equation (2), then the nonlinear mappings and the reduced dynamics are
computed, by applying the procedure that is illustrated by Equations (4) to (6). For the sake of illustration, the normal dynamics up to order 7 is shown here:

$$ \dot{z}_1 = i\omega z_1 + i \frac{-10g^2 + 9h\omega^2}{6\omega^3} z_1^2 z_2 + i \frac{-3140g^4 + 8388g^2h\omega^2 - 1377h^2\omega^4}{432\omega^7} z_1^3 z_2^2 + i \frac{-523960g^6 + 2186724g^4h\omega^2 - 1913274g^2h^2\omega^4 + 114939h^3\omega^6}{10368\omega^{11}} z_1^4 z_2^3. \quad (26) $$

while the associated nonlinear mapping is reported in Appendix D. For the sake of brevity, only the equation for $z_1$ is shown in Equation (26), the second equation for $z_2$ being simply its complex conjugate. Rewriting the coefficients appearing in Equation (26) as $if^{(p)}$, for odd p corresponding to the odd monomial remaining:

$$ \dot{z}_1 = if^{(1)}z_1 + if^{(3)}z_1^2 z_2 + if^{(5)}z_1^3 z_2^2 + \dots, \quad (27) $$

the backbone is analytic and reads, using a polar representation for ($z_1, z_2$) as in Section 2.1, Equation (9):

$$ \omega_{\text{NL}} = f^{(1)} + f^{(3)} \left(\frac{\rho}{2}\right)^2 + f^{(5)} \left(\frac{\rho}{2}\right)^4 + \dots \quad (28) $$

Given the expressions of the coefficients $f^{(p)}$ shown in Equation (26), one can see that, assuming $h \ge 0$, the quadratic and cubic nonlinearities play opposite roles in defining the hardening/softening behaviour. This has been used for a long time to predict, only from the sign of $f^{(3)}$, the type of nonlinearity for structures, using for example the oscillator normal form, see e.g. (Touzé et al. 2004; Touzé and Thomas 2006; Touzé and Amabili 2006; Touzé et al. 2008). Arbitrary order expansions allow a finer understanding of the hard/soft transition, and have already been used for rotating beams in (Martin et al. 2023).

Figure 3 shows the behaviour of the coefficients $f^{(p)}$ for $p = 3, 5, 7$ and $9$. $f^{(3)}$ changes sign only once (for the case under study, with $h = 1$ fixed and varying $g$), meaning that the transition from hardening to softening behaviour occurs at $g = 3\sqrt{h}/10$. On the other hand, higher-order coefficients are polynomials of higher degrees in the coefficients $g$ and $h$ and thus have numerous zeros, see Figures 3(a) and 3(b). Figure 3(c) shows the backbone curves obtained just before and after the transition of the cubic coefficient (while the next orders do not change sign), namely for $g = 0.85$ ($f^{(3)} > 0$, and $f^{(5)} > 0$, $f^{(7)} < 0$, $f^{(9)} < 0$), and $g = 1$ ($f^{(3)} < 0$, and $f^{(5)} > 0$, $f^{(7)} < 0$, $f^{(9)} < 0$). A reference solution obtained numerically by continuation is compared to two truncations, respectively to orders 7 and 9, of the analytical backbone curve. For $g = 0.85$, the negative signs of $f^{(7)}$ and $f^{(9)}$ change the high amplitude behaviour of orders 7 and 9 truncations that depart from reference, needing orders higher than 9 to achieve convergence up to the selected amplitudes. For $g = 1$, the positive coefficient $f^{(5)}$ plays the major role since $f^{(3)}$ is close to zero.

Two other cases where the dominant behaviour is softening are shown in Figure 3d for $g = 1.5$ and $g = 1.65$, around the change of sign of $f^{(5)}$. Again, in such a situation, the mixed signs of the different truncations ask for high-order development to reach convergence. Nevertheless, the qualitative change of behaviours is captured. Also, the validity limit of the normal form development is probably impacted by the different values of quadratic and cubic coefficients.

For the sake of completeness, the RNF and ONF analysis are reported in Appendix E. In these two cases, automatic analysis of the backbone curve up to high order is more difficult.

## 3.2 Effect of forcing and damping

In this section, an analysis of the forced and damped Duffing oscillator with high-order normal form expansions is developed. The goal is to show how normal form approach relates to and extends the established procedure in nonlinear vibration theory, based on perturbation methods. All these results will then be helpful to understand and analyse reduced-order models with a single nonlinear normal mode assumption. The starting point is thus a forced-damped Duffing oscillator

$$ \ddot{u} + \omega^2 u + 2\xi\omega\dot{u} + hu^3 = \frac{\kappa}{2} (e^{+i\Omega t} + e^{-i\Omega t}), \quad (29) $$
<figure><img src="image_4.png" /><figcaption>**Figure 3** Hardening/softening behaviour for the Duffing oscillator with quadratic and cubic nonlinearities. A complex normal form style is used with parameter values set as $\omega = 1$, $h = 1$. (a) Behaviour of the coefficients $f^{(p)}$ as a function of $g$. (b) Zoom on the area highlighted by the black dashed box in (a). The vertical dashed lines indicate the values of $g$ for which the backbones have been calculated. (c-d) Backbone curves for different values of $g$. A reference solution obtained by numerical continuation is compared to the backbone curves obtained with CNF up to orders $o = 7$ and $o = 9$.</figcaption></figure>

where the forcing has already been rewritten with complex notation, while a viscous damping is introduced with the damping ratio $\xi$. We start by considering the case of a forcing frequency $\Omega$ which is far from the primary and secondary resonances.

### 3.2.1 Non-resonant excitation

The case of an out-of-resonance forcing frequency $\Omega$ leads to a minimal number of resonant monomials in the normal form, since only the trivial resonances are considered. Thanks to the efficient treatment of the non-autonomous forcing term presented in (Vizzaccaro et al. 2024), an arbitrary order expansion in terms of the power of the non-autonomous term is also implemented in MORFE_Symbolic. This key feature allows computing and analysing the high-order terms produced by the forcing in the normal form, a case that is not considered for example in (Breunung and Haller 2018; Opreni et al. 2023; Jain and Haller 2022), where only linear terms of the forcing are included.

The method proposed in (Vizzaccaro et al. 2024) to deal with the forcing term is to make the system autonomous by adding two additional coordinates, namely $z_3 = e^{i\Omega t}$ and $z_4 = e^{-i\Omega t}$, augmenting the size of the original system but keeping the added coordinates unchanged. More precisely, the state variables of the augmented system become $\mathbf{y} = [u, v, z_3, z_4]^T$, while the augmented normal coordinates are $\mathbf{z} = [z_1, z_2, z_3, z_4]^T$, such that only the first two physical coordinates need to be sought as a parametrisation of the normal ones, i.e.

$$ \begin{bmatrix} u \\ v \end{bmatrix} = W(z), \qquad (30) $$

with $W(z)$ denoting the nonlinear mappings. Additionally, because of the definition of the added variables, their reduced dynamics is trivial, and already known:

$$ \dot{z}_3 = \Omega z_3, \qquad (31a) $$

$$ \dot{z}_4 = -\Omega z_4. \qquad (31b) $$

This leads to an efficient reformulation of the parametrisation method asking for slight modifications in the algorithm as compared to the autonomous case. When a single master mode (NNM)
is selected, the reduced dynamics depends on four coordinates ($z_1, z_2, z_3, z_4$), which can be easily reduced to ($z_1, z_2$) only by replacing $z_3 = e^{+i\Omega t}$ and $z_4 = e^{-i\Omega t}$.

Another key feature of the method proposed in (Opreni et al. 2023; Vizzaccaro et al. 2024), to handle the external periodic forcing, is to compute the parametrisation for only one value of the external forcing. Rigorously speaking, the parametrisation needs to be computed for each forcing frequency, as proposed, for example, in (Jiang et al. 2005), to produce exact solutions. However, to alleviate the associated computational burden, the parametrisation can be computed for a single excitation frequency $\Omega_p$, and this ROM can be used to draw out rapidly frequency response curves (FRC), assuming that the dependence upon $\Omega$ is small. This strategy is also adopted here to propose analytical solutions obtained with symbolic computations that lend themselves well to analysis.

The complex normal form (CNF) up to order three for the forced, undamped ($\xi = 0$) Duffing oscillator Equation (29), and for this out-of-resonance case, is computed at a single excitation frequency $\Omega_p = \Omega$. It reads

$$ \dot{z}_1 = i\omega z_1 + i \frac{3h}{2\omega} z_1^2 z_2 + i \frac{3h\kappa^2}{4\omega(\Omega^2 - \omega^2)^2} z_1 z_3 z_4, \quad (32) $$

where again only the first equation is shown, the second equation for $z_2$ being the complex conjugate. At this order, only two resonant monomials are present in the CNF. The first one, $z_1^2 z_2$, has already been commented on in Section 2.1 in Equation (8). The second one, $z_1 z_3 z_4$, depends on the forcing, and scales with the amplitude of the forcing squared, $\kappa^2$. This monomial is very interesting since it appears only due to the complete treatment of the forcing shown in (Vizzaccaro et al. 2024). Since $z_3 = e^{+i\Omega t}$ and $z_4 = e^{-i\Omega t}$, the monomial can be interpreted as having a direct consequence on the nonlinear oscillation frequency. In particular, it shows the dependence of the free oscillation frequencies upon the forcing amplitude. For such conservative dynamics, the system's response is characterised by a quasi-periodic behaviour, stemming from the contributions of the free oscillations and of the forced response, whose effect is completely embedded into the nonlinear mapping equations. As commented next, the free oscillation term is generally related to the transient and is damped out when losses are taken into account. Interestingly, the monomial $z_1 z_3 z_4$ corresponds to a trivial resonance and will thus be present in all normal form dynamics that are considered in the next sections (primary and secondary resonance). It is thus further analysed, in particular for secondary resonances, since in this case the forcing is not assumed to be small such that the dependence of the oscillation frequencies upon forcing amplitude is not negligible.

The nonlinear mapping up to order 3, that corresponds to Equation (32), reads:

$$ u = (z_1 + z_2) + \left( \frac{h}{8\omega^2} z_1^3 - \frac{3h}{4\omega^2} z_1^2 z_2 - \frac{3h}{4\omega^2} z_1 z_2^2 + \frac{h}{8\omega^2} z_2^3 \right) - \left( \frac{1}{(\Omega^2 - \omega^2)^2} \frac{\kappa}{2} (z_3 + z_4) \right) \\ - \left( \frac{3h\kappa}{2(\Omega^2 - \omega^2)(\Omega + \omega)(\Omega + 3\omega)} z_1^2 z_3 + \frac{3h\kappa}{2(\Omega^2 - \omega^2)(\Omega - \omega)(\Omega - 3\omega)} z_1^2 z_4 \right. \\ \left. + \frac{3h\kappa}{(\Omega^2 - \omega^2)^2} z_1 z_2 (z_3 + z_4) \right) \quad (33) $$

In the above equation, note that the sums and subtractions of variables $z_3$ and $z_4$ can be written as cosines and sines. The velocity mapping is not shown here for the sake of brevity, and also because it can be retrieved from the derivation of Equation (33) with respect to time; see the remark in Appendix B. Note also that Equation (33) has been truncated to $O(\kappa^2)$ terms only to make it more concise.

Equation (33) shows that the autonomous deformation of the associated invariant manifolds, already present in Equation (B.1), is recovered but there are also additional non-autonomous deformations, driven by the added coordinates $z_3$ and $z_4$. The linear terms in $z_3$ and $z_4$ in the first line of Equation (33) account for a rigid-body rotation of the invariant manifold, already commented on in (Opreni et al. 2023). Then the terms of the next two lines show that together with this rigid-body motion, the manifold shows deformations along the phase of the forcing.
The symbolic code MORFE_Symbolic can be used to derive higher-order approximations of the normal form and nonlinear change of coordinates in this case of non-resonant excitation. Interestingly, since only trivially resonant monomials are in the normal form, the dynamics will be composed of two different terms. First, the autonomous terms with resonant monomials will strictly follow those appearing in the unforced case, see Equation (8) for an example up to order 9. In addition to these, all the trivially resonant monomials will reappear, multiplied by a factor of the form $z_3^p z_4^q$ with $p \ge 1$. Higher order results are not shown here but can be automatically produced via MORFE_Symbolic.

Finally, the effect of the viscous damping term is investigated to close this case of non-resonant excitation. With the damping, the eigenvalues of Equation (29) are $\lambda_{1,2} = -\xi\omega \pm i\omega\sqrt{1 - \xi^2} = \omega(-\xi \pm i\delta)$, see Equation (15). The CNF up to order three for the non-resonant response is then

$$ \dot{z}_1 = (i\delta\omega - \xi\omega) z_1 + i \frac{3h}{2\delta\omega} z_1^2 z_2 + \frac{3ih\kappa^2}{4\delta\omega ((\Omega^2 - \omega^2)^2 + 4\xi^2\omega^2\Omega^2)} z_1 z_3 z_4, \quad (34) $$

the second equation for $z_2$ being the complex conjugate of this one. The real part in the first monomial shows that these oscillations are damped out. Consequently, as commented in (Opreni et al. 2023), the reduced dynamics is trivial and all the effects are embedded in the nonlinear mapping, where only the terms related to $z_3$ and $z_4$ remain.

### 3.2.2 Primary resonance

The case of a resonant forcing with $\Omega \approx \omega$ is discussed next. In this case of primary resonance, we compute the parametrisation at the value of the damped oscillation frequency $\delta\omega$ (imaginary part of the eigenvalue), such that $\Omega_p = \delta\omega$. With this choice and taking the damping into account, the complex normal form (CNF) at primary resonance up to order 3 becomes

$$ \dot{z}_1 = f_1 z_1 + f_2 z_1^2 z_2 + f_3 z_3 + f_4 z_1^2 z_4 + f_5 z_1 z_2 z_3 + f_6 z_1 z_3 z_4 + f_7 z_2 z_3^2 + f_8 z_3^2 z_4, \quad (35) $$

where only the first equation is given, the second on $z_2$ being the complex conjugate. The coefficients $f_j, j = 1, \dots, 8$, read:

$$ f_1 = \lambda_1, \quad f_2 = i \frac{3h}{2\delta\omega}, \quad f_3 = -i \frac{\kappa}{4\delta\omega}, \quad f_4 = i \frac{3h\kappa}{8\delta^2\omega^3 (2\delta + i\xi)} \quad (36a) $$

$$ f_5 = i \frac{3h\kappa}{4\delta^2\omega^3 (2\delta - i\xi)}, \quad f_6 = i \frac{3h\kappa^2}{16\delta^3\omega^5 (4\delta^2 + \xi^2)}, \quad (36b) $$

$$ f_7 = i \frac{3h\kappa^2}{32\delta^3\omega^5 (2\delta - i\xi)^2}, \quad f_8 = i \frac{3h\kappa^3}{128\delta^4\omega^7 (2\delta - i\xi)^2 (2\delta + i\xi)}. \quad (36c) $$

Additionally, the nonlinear mappings for this case are given in Appendix F. Note that, as opposed to the non-resonant case, the coefficients of the normal form are not a function of the forcing frequency $\Omega$. This is a direct consequence of the choice $\Omega_p = \delta\omega$. A different choice could be introduced here, for example by selecting $\Omega_p = \omega$, which appears closer to the assumption made in perturbative techniques. Indeed, these two choices are almost indistinguishable for most practical cases where the small damping assumption is met such that $\delta \approx 1$. Since the analysis developed here does not necessarily assume small damping, it appeared more coherent to expand around $\Omega_p = \delta\omega$, which also leads to simpler expressions for the coefficients. For comparison purposes, the coefficients obtained with the expansion around $\Omega_p = \omega$ are given in Appendix G.

The first two monomials in Equation (35) are respectively the linear term and the cubic nonlinear stiffness that were already present without forcing. The third term corresponds to the direct resonant forcing. The next two monomials, $z_1^2 z_4$ and $z_1 z_2 z_3$, correspond to nonlinear parametric-like excitation terms. In the sequence, $z_1 z_3 z_4$, is the trivially resonant monomial already commented on in the previous section. Its coefficient $f_6$ scales as $O(\kappa^2)$, like the next monomial that introduces a second harmonic forcing through $z_3^2$. Finally, the last coefficient $f_8$ scales as $O(\kappa^3)$ and is a higher-order effect of the direct forcing.

These reduced dynamics are automatically derived in the context of the parametrisation method up to high-order and are used to analyse the primary resonance of finite element models,
see e.g. (Vizzaccaro et al. 2022; Opreni et al. 2023; Vizzaccaro et al. 2024; Jain and Haller 2022; Li et al. 2022). The purpose here is to give some insights into these results, which are obtained using high-order solutions that are realified with polar or cartesian representations, by analysing the wealth of all the terms involved in the solution. Note that for obtaining the normal form given in Equation (35), no specific assumptions on damping or forcing have been made and the solutions are not stopped at first order for these two terms, as commonly assumed in perturbative solutions. Hence the solution given by Equation (35) contains a priori more information.

To give more insights to Equation (35), polar coordinates are introduced as

$$z_{1,2} = \frac{\rho}{2} e^{\pm i\theta}, \quad z_{3,4} = e^{\pm i\Omega t}. \qquad (37)$$

An autonomous system can be derived for the amplitude $\rho$ and the phase $\psi$, defined as $\psi = \theta - \phi$, where $\phi = \Omega t$ is introduced to make the system autonomous with $\dot{\phi} = \Omega$. The resulting system reads

$$\dot{\rho} = A_0 + A_1^c \cos \psi + A_1^s \sin \psi + A_2^c \cos (2\psi) + A_2^s \sin (2\psi), \qquad (38a)$$

$$\rho \dot{\psi} = B_0 + B_1^c \cos \psi + B_1^s \sin \psi + B_2^c \cos (2\psi) + B_2^s \sin (2\psi). \qquad (38b)$$

The coefficients in Equation (38) have explicit expressions as functions of the coefficients $f_j$ in Equation (35), which have been split into to their real part $f_j^R$ and imaginary part $f_j^I$. The expressions read:

$$\begin{align}
A_0 &= \rho f_1^R + \frac{\rho^3}{4} f_2^R + \rho f_6^R, & B_0 &= \rho f_1^I + \frac{\rho^3}{4} f_2^I + \rho f_6^I - \rho \Omega, \nonumber \\
A_1^c &= 2 f_3^R + \frac{\rho^2}{2} f_4^R + \frac{\rho^2}{2} f_5^R + 2 f_8^R, & A_1^s &= 2 f_3^I - \frac{\rho^2}{2} f_4^I + \frac{\rho^2}{2} f_5^I + 2 f_8^I, \nonumber \\
B_1^c &= 2 f_3^I + \frac{\rho^2}{2} f_4^I + \frac{\rho^2}{2} f_5^I + 2 f_8^I, & B_1^s &= -2 f_3^R + \frac{\rho^2}{2} f_4^R - \frac{\rho^2}{2} f_5^R - 2 f_8^R, \tag{39} \\
A_2^c &= \rho f_7^R, & A_2^s &= \rho f_7^I, \nonumber \\
B_2^c &= \rho f_7^I, & B_2^s &= -\rho f_7^R. \nonumber
\end{align}$$

Deriving a closed-form expression for the FRC from the fixed points of Equation (38) is difficult because of the presence of the second harmonic excitation, which creates the terms with arguments $2\psi$. From the analysis of the monomials in Equation (35), the second harmonic terms are only created by $z_3^2$, such that the coefficients $A_2^c, A_2^s, B_2^c$ and $B_2^s$ are directly proportional to $f_7$. To simplify the analysis, the usual assumption that a small forcing is enough to lead to large amplitude solutions in case of primary resonance can be made, such that one could neglect the terms in $\kappa^2$ and $\kappa^3$ in Equation (35). With this assumption, which leads to disregarding $f_6, f_7$ and $f_8$, and using $\sin^2 \psi + \cos^2 \psi = 1$, one can obtain the following relationship as an expression for the frequency response curve:

$$\left( \frac{A_0}{A_1^s} + \frac{A_1^c}{A_1^s} \frac{A_0 B_1^s - A_1^s B_0}{A_1^s B_1^c - A_1^c B_1^s} \right)^2 + \left( \frac{A_0 B_1^s - A_1^s B_0}{A_1^s B_1^c - A_1^c B_1^s} \right)^2 = 1. \qquad (40)$$

With this equation, it is still too cumbersome to derive an explicit expression for the amplitude $\rho$ as a function of the excitation frequency $\Omega$ for direct comparisons with a perturbative solution. To proceed, we assume that the damping $\xi$ is small, as it is usually done in a perturbative scheme. Neglecting all high-order terms of the damping in the expressions of the coefficients $f_j$, $j = 1, \dots, 5$; leads to these expressions:

$$f_1^R = -\xi\omega, \quad f_1^I = \omega, \quad f_2^R = 0, \quad f_2^I = \frac{3h}{2\omega}, \quad f_3^R = 0, \qquad (41a)$$

$$f_3^I = -\frac{\kappa}{4\omega}, \quad f_4^R = 0, \quad f_4^I = \frac{3h\kappa}{16\omega^3}, \quad f_5^R = 0, \quad f_5^I = \frac{3h\kappa}{8\omega^3}, \qquad (41b)$$

such that an explicit expression for the FRC, with a small damping assumption, reads

$$\rho \Omega = \rho \omega + \rho^3 \frac{3h}{8\omega} \pm \frac{9h\rho^2 - 16\omega^2}{3h\rho^2 - 16\omega^2} \sqrt{\left(\frac{3h\rho^2 - 16\omega^2}{32\omega^3}\kappa\right)^2 - \rho^2 \xi^2 \omega^2}. \qquad (42)$$
This expression is still more complete than the one given by a first-order perturbative solution like the method of multiple scales (MMS). In particular, it involves $f_4$ and $f_5$ coefficients, which are linked to the monomials $z_1^2 z_4$ and $z_1 z_2 z_3$ in Equation (35). Note that these two monomials refer to a nonlinear parametric excitation that is generally overlooked. Neglecting these two terms leads to considering only the first three monomials in Equation (35), and to the classical first-order perturbative solution of the FRC, see e.g. (Nayfeh and Mook 1995), page 166. This can be justified on the coefficients by assuming the same hypotheses as in the perturbative solution. Namely that, additionally to $\kappa$ and $\xi$, the nonlinear coefficient $h$ is also small, such that when expanding Equation (42) in power series and considering only the leading order terms, the classical result is found:

$$ \Omega = \omega + \rho^2 \frac{3h}{8\omega} \pm \sqrt{\frac{\kappa^2}{4\rho^2\omega^2} - \xi^2\omega^2}. \quad (43) $$

This development shows the wealth of the automated high-order solution provided by the complex normal form, which can be simplified in order to retrieve the results of the first-order multiple scales expansion. Of course, higher-order perturbative expansions should give back the same term, but the change of paradigm in the approach proposed here is to give automatically high-order expansions that can then be analysed by introducing asymptotic and ordering between the different terms, leading to explicit analytical expressions in the simplest cases.

Figure 4 illustrates the previous derivations by comparing analytical and numerical results. Figure 4(a) first compares analytical expressions based on the different approximations. The first-order multiple scales solution, Equation (43), is compared to the solution given by Equation (42), for classical parameter values corresponding to usual assumptions: $\omega = 1.5$, $h = 1$, small damping with $\xi = 0.02$, and $\kappa = 0.1$, leading to a maximal vibration amplitude $\rho$ close to 1, in order to see the effect of moderate amplitudes. It shows that neglecting the two monomials $z_1^2 z_4$ and $z_1 z_2 z_3$ in Equation (35) has an important effect on the prediction of the maximum amplitude, underlining that, when $\rho$ is close to 1, a first-order perturbative solution is not accurate enough. The explicit expression Equation (42) is also compared to the implicit analytical expression given by Equation (40) and with an implicit solution for the fixed points of Equation (38), computed symbolically after replacing numerical values for the system parameters, showing that, at this level of amplitude and for these parameter values, the assumptions leading to Equation (42) are accurate enough.

<figure><img src="image_5.png" /><figcaption>Figure 4 FRCs for the primary resonance of the damped cubic Duffing oscillator. Parameter values are $\omega = 1.5$, $h = 1$, $\xi = 0.02$ and $\kappa = 0.1$. (a) Comparison of the analytical FRCs given by Equations (40), (42) and (43) and one obtained by symbolically solving for the fixed points of Equation (38) after replacing the numerical values of the parameters. (b) Comparison of numerically computed FRCs, obtained from MORFE_Symbolic using Matcont (Dhooge et al. 2004), with a reference solution, also obtained by continuation, and with the first-order multiple scale solution given by Equation (43), denoted as MMS.</figcaption></figure>

Figure 4(b) compares numerical solutions to show the convergence of the high-order normal form, as well as the effect of the nonlinear mapping. For the sake of illustration, the first-order
multiple scales solution is also shown. Whereas the amplitude reported in Figure 4(a) is that of the normal coordinate $\rho$, Figure 4(b) is given for the maximum displacement of the physical coordinate $u$, denoted as $u_{\max}$. Note that, as compared to the calculations for the backbone curves shown in Section 2.1, it is not possible here to recover an easy analytic expression relating $\rho$ to the physical displacement through the nonlinear mapping, because the forcing terms are not in phase with the autonomous ones. It should be highlighted, however, that the analytical expression (first-order perturbative solution) given by Equation (43) assumes (at first order), $u_{\max} = \rho$. The numerical solutions are computed numerically with a continuation procedure embedded in the package Matcont (Dhooge et al. 2004). Four different truncations of the normal form expansion are shown. Following (Vizzaccaro et al. 2024), the selected truncation order is denoted as $O(z^p, \epsilon^q)$, meaning a maximal order $p$ in the $z = (z_1, z_2, z_3, z_4)$ coordinate and $q$ for the non-autonomous $(z_3, z_4)$ variables. In this example, the case $p = q$ is selected for convenience, and $p = q = 3$ corresponds to Equation (35). These results show that, for such level of amplitude, the first-order multiple scales solution is far from the full-order solution. The order 3 expansion is also not accurate enough, while convergence is almost achieved at order 9. Besides, the effect of the nonlinear mapping is very important, as can be seen by comparing the order 3 solutions when parameterised in amplitude by $\rho$ or the maximum displacement $u_{\max}$.

All this development underlines the wealth of the symbolic solution, which can be analysed in order to produce, whenever possible, analytical expressions that can be more complete than first-order perturbative solution. This also highlights the fact that, when amplitudes are close to 1, one needs to resort to numerical approximations since low-order approximations are not accurate enough.

### 3.2.3 Superharmonic resonance

In this section, the 3 : 1 superharmonic resonance of the Duffing oscillator is investigated thanks to high-order normal form approximations and analytical expressions. Equation (29) is considered as starting point and the forcing frequency $\Omega$ is such that $\Omega \approx \omega/3$. The complex normal form (CNF) up to order three is obtained from MORFE_Symbolic, and the expansion point for computing the parametrisation has been selected as $\Omega_p = \omega/3$. The choice $\Omega_p = \delta\omega/3$ has not been selected in this case since it does not yield substantial simplifications to the coefficients. It reads:

$$ \dot{z}_1 = \lambda_1 z_1 + i \frac{3h}{2\delta\omega} z_1^2 z_2 + \frac{243i h\kappa^2}{16\delta\omega^5 (9\xi^2 + 16)} z_1 z_3 z_4 - \frac{729h\kappa^3}{128\delta\omega^7 (3\xi - 4i)^3} z_3^3, \quad (44) $$

where the notations of the previous sections are used, and the second equation for $z_2$ has been omitted since it is the complex conjugate. The associated nonlinear mapping is given in Appendix H for the sake of completeness. As compared to the primary resonance, Equation (44) contains fewer terms. On the other hand, only one more resonant monomial as compared to the non-resonant case shown in Equation (32), is present. The added monomial is $z_3^3$, which is indeed the resonant term due to 3 : 1 superharmonic resonance, since $z_3 = e^{i\Omega t}$ and $\Omega \approx \omega/3$.

At this order three of the development, no terms in the reduced dynamics involve powers of the forcing coordinates, such as $z_3^2$ that was for example in Equation (35) for the primary resonance. Consequently, no higher harmonics of the forcing will appear in the solution, such that analytical expressions for the frequency response curve could be easily derived at this order. Using the polar coordinates, redefining the phase as $\psi = \theta - 3\phi$ in order to render the system autonomous, and searching for fixed points such that $\dot{\rho} = \dot{\psi} = 0$, the following relationship is derived:

$$ B_r + A_i \sin(\psi) + A_r \cos(\psi) = 0, \quad (45a) $$

$$ B_i - A_r \sin(\psi) + A_i \cos(\psi) = 0, \quad (45b) $$

where $A_r$, $A_i$, $B_r$ and $B_i$ have explicit expressions as functions of the coefficients of the normal
dynamics Equation (44):

$$ 
\begin{aligned}
A_r &= \frac{6561h\kappa^3\xi(16 - 3\xi^2)}{64\delta(9\xi^2 + 16)^3\omega^7}, & A_i &= -\frac{729h\kappa^3(27 - 6\delta^2)}{16\delta(9\xi^2 + 16)^3\omega^7} \\
B_r &= -\rho\xi\omega, & B_i &= \rho\delta\omega - 3\rho\Omega + \rho^3\frac{3h}{8\delta\omega} + \rho\frac{243h\kappa^2}{16\delta(9\xi^2 + 16)\omega^5}.
\end{aligned}
\quad (46) $$

Squaring each line of Equation (45), summing and using $\sin^2 \psi + \cos^2 \psi = 1$ in order to eliminate the angle $\psi$ leads to the following relationship that gives the frequency response function, i.e. the amplitude $\rho$ as a function of the forcing frequency $\Omega$:

$$ A_i^2 + A_r^2 = B_i^2 + B_r^2. \quad (47) $$

From this last equation, an explicit expression for the frequency response curve is derived, and the coefficients from Equation (46) are substituted to obtain:

$$ 3\rho\Omega = \rho\delta\omega + \rho^3\frac{3h}{2^3\delta\omega} + \rho\frac{3^5h\kappa^2}{2^4\delta(9\xi^2 + 16)\omega^5} \pm \sqrt{\frac{3^{12}h^2\kappa^6}{2^{12}\delta^2\omega^{14}(9\xi^2 + 16)^3} - \rho^2\xi^2\omega^2}. \quad (48) $$

This expression for the FRC of the superharmonic case is very close to the one obtained using a first-order multiple scales expansion. It is however more general since no small damping assumption has been made yet. Using a first-order expansion on the damping term $\xi$, following the guideline used in the previous section, gives:

$$ 3\rho\Omega = \rho\omega + \rho^3\frac{3h}{2^3\omega} + \rho\frac{3^5h\kappa^2}{2^8\omega^5} \pm \sqrt{\frac{3^{12}h^2\kappa^6}{2^{24}\omega^{14}} - \rho^2\xi^2\omega^2}, \quad (49) $$

this last expression being exactly equivalent to the one reported in (Nayfeh and Mook 1995) using the first-order method of multiple scales (MMS). At this level of the asymptotic expansion (order 3), the normal form of the superharmonic solution is thus equivalent to a first-order perturbative solution. This is logical because in the present situation, the amplitude $\kappa$ can no longer be considered small, an effect that is accounted for in the multiple scales solution by assuming that the forcing appears at the order zero of the solution, since a secondary resonance is considered.

The terms which are outside the square root in Equation (48), represent a shift of the traditional backbone curve with forcing amplitude. Interestingly, this dependence only comes from the monomial $z_1z_3z_4$. The expression that describes this curve is

$$ 3\omega_{\text{NL}} = \delta\omega + \rho^2\frac{3h}{2^3\delta\omega} + \frac{3^5h\kappa^2}{2^4\delta(9\xi^2 + 16)\omega^5}. \quad (50) $$

When $\kappa = 0$, this last equation recovers the backbone curve of the primary resonance (see Equation (18b) truncated at the third-order), which is simply shifted and scaled to one-third of the eigenfrequency. For a fixed value of $\kappa \neq 0$, it aligns to a backbone passing through the maximum of the FRC, as illustrated in Figure 5(a).

Additionally, it is possible to derive an expression for the curve that connects the maxima of the superharmonic FRCs, which can be viewed as a generalised backbone curve for this resonance scenario. It is obtained as a curve parameterised by $\kappa$: for each forcing amplitude, there is a value for $\rho$ such that the frequency in Equation (48) is single-valued, and a corresponding value of $\Omega$:

$$ 
\begin{cases}
\rho(\kappa) = \sqrt{\frac{3^{12}h^2\kappa^6}{2^{12}\delta^2\omega^{16}\xi^2(9\xi^2 + 16)^3}}, \\
\Omega(\kappa) = \frac{\delta\omega}{3} + \frac{3^{12}h^3\kappa^6}{2^{15}\delta^3\omega^{17}\xi^2(9\xi^2 + 16)^3} + \frac{3^3h\kappa^2}{2^3\delta(9\xi^2 + 16)\omega^5}.
\end{cases}
\quad (51) $$

Equations (48) to (49) are used to plot the FRCs in Figure 5(a) for three different values of forcing amplitude: $\kappa \in \{0.5, 0.7, 0.9\}$. The other parameters are fixed as $\omega = 1$, $\xi = 0.2$, $h = 3$. The figure also shows different shifted backbones for the problem. The generalized backbone
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5 FRCs for the 3:1 superharmonic resonance of the damped cubic Duffing oscillator. The results are obtained using a complex normal form style, and parameter values are set as $\omega = 1$, $h = 3$, $\xi = 0.2$. (a) FRCs in the frequency-amplitude of the normal coordinate plane. Note how the forcing-dependent backbones follow the shift of the curves with increasing $\kappa$. In the legend, the grey lines indicate which equation has been used to calculate each curve, while different colors are associated with different $\kappa$ values. For $\kappa = 0$ there is no FRC, only the forcing-dependent backbone. (b) Three-dimensional view of the FRCs and of the frequency response manifold. The curve that unites the peaks of the FRCs is also shown as a red dash-dotted line.</figcaption></figure>

curve, following the peaks of the FRCs and given by Equation (51), is shown in red dash-dotted line. In Figure 5(a), the two expressions for the FRC, respectively without assumption on the damping, Equation (48), or assuming a first-order, Equation (49), are also shown, underlining the minimal difference that is brought about by considering accurately the damping ratio as compared to the first-order perturbative solution. Note that, for illustrative purpose, values of $\rho$ largely exceeding 1 have been selected. Furthermore, the range of frequencies shown is large and exceeds the vicinity of the superharmonic resonance, with values approaching the region where the primary resonance influences the response. Therefore, the illustrative character of Figure 5(a) should be kept in mind, where parameter values such that a clear shift of the FRCs could be observed were chosen purposefully.

The curves are also plotted in 3D space in Figure 5(b), by adding the forcing amplitude $\kappa$ as an additional coordinate. A frequency response manifold, given by Equation (47) is represented in order to highlight the set on which the solutions are. The generalised backbone curve given by Equation (51) is represented as the red dash-dotted line joining all the maxima of this surface.

Figure 6 completes the analysis for this case, presenting the FRCs in terms of physical displacement amplitude, $u_{max}$, and for higher orders of parametrisation. The curves therein are obtained, such as in the primary resonance case, by numerical continuation with the package Matcont (Dhooge et al. 2004), since it is not possible to find an analytical expression relating $u_{max}$ to $\rho$. The figure also shows the multiple scales solution given by Nayfeh and Mook (1995, p. 176). It should be noted that, in contrast with the case of primary resonance, this solution does not simply assume that $u_{max} = \rho$. Indeed, the time series of the displacement is given as

$$u = \rho \cos (3\Omega t - \gamma) + \frac{\kappa}{\omega^2 - \Omega^2} \cos \Omega t, \qquad (52)$$

with $\gamma$ being a phase defined in the reference, in the same page mentioned before. The curve corresponding to this solution was found by calculating the maxima of this expression over a period of oscillation for different values of $\rho$. Note that, since the damping is small, the difference between this curve and the one of the $O(z^3, \varepsilon^3)$ parametrisation is only due to the nonlinear mappings relating $u$ to $\rho$. Consequently, even though the same FRC in terms of $\rho$ is found for both of them, the approximation given by the normal form computation up to order 3 is more accurate in this case, thanks to the correction brought about by the mapping. Nevertheless, an
order three expansion is still not enough when compared to the response of the full system, i.e. higher order terms need to be taken into account to converge to the exact solution, which happens at order seven.

<figure><img src="image_8.png" /><figcaption>**Figure 6** FRCs for the 3:1 superharmonic resonance of the damped cubic Duffing oscillator. Parameter values are $\omega = 1.5$, $h = 1$, $\xi = 0.002$ and $\kappa = 0.5$. Comparison of numerically computed FRCs, obtained from MORFE_Symbolic and solved using numerical continuation, with a reference solution, also obtained by continuation but on the original equations, and with the first-order multiple scales solution (MMS).</figcaption></figure>

At this point it is interesting to highlight that the solution obtained by the CNF converges to the reference one only in the vicinity of the superharmonic resonance, but starts to show slight discrepancies outside this region, see e.g. for $\Omega/\omega > 0.35$ in Figure 6. This is inherent to the local nature of the normal form solutions, that are developed for specific resonances, and which should not be expected to converge for the whole range of forcing frequencies. Specifically, the reference solution starts to be affected by the primary resonance when $\Omega/\omega > 0.35$. This effect is not captured by the normal form solution, no matter the maximal order chosen for the parametrisation.

### 3.2.4 Subharmonic resonance

In this section, the 1:3 subharmonic resonance, where the excitation frequency is in the vicinity of three times the natural frequency $\omega$, is investigated, for the Duffing equation with cubic nonlinearity, Equation (29). The complex normal form (CNF) up to the third order for both autonomous and non-autonomous variables, with the expansion point for the parametrisation being $\Omega_p = 3\omega$ is

$$ \dot{z}_1 = \lambda_1 z_1 + \frac{3ih}{2\delta\omega} z_1^2 z_2 + \frac{3ih\kappa^2}{16\delta\omega^5 (16 + 9\xi^2)} z_1 z_3 z_4 - \frac{3h\kappa(3\xi - 4i)}{8\delta\omega^3 (16 + 9\xi^2)} z_2^2 z_3. \quad (53) $$

Again, the equation for $z_2$ is not reported for the sake of brevity, being the complex conjugate. The nonlinear mapping up to order 3 is reported in Appendix I. The normal form dynamics looks similar to the superharmonic case, only the last monomial being changed, since now the resonant monomial with the assumption $\Omega \approx 3\omega$ is $z_2^2 z_3$. The term $z_1^2 z_2$ refers to the cubic nonlinearity, while the second monomial $z_1 z_3 z_4$ is the trivially resonant term, scaling as the square of the forcing amplitude $\kappa^2$, and making clearly appear the hard non-resonant excitation that is key for secondary resonances.

Proceeding similarly as in the superharmonic case to find an analytic expression for the FRC, the system to solve can be put in the following form, this time with $\psi = 3\theta - \phi$:

$$ B_r + A_i \sin(\psi) + A_r \cos(\psi) = 0, \quad (54) $$

$$ B_i - A_r \sin(\psi) + A_i \cos(\psi) = 0, \quad (55) $$

with:

$$ A_r = -\rho^2 \frac{3h\kappa\xi}{2^4\delta\omega^3 (16 + 9\xi^2)}, \quad A_i = \rho^2 \frac{3h\kappa}{2^2\delta\omega^3 (16 + 9\xi^2)}, \quad (56a) $$
$$B_r = -\rho\xi\omega, \quad B_i = \rho\delta\omega - \rho\frac{\Omega}{3} + \rho^3\frac{3h}{2^3\delta\omega} + \rho\frac{3h\kappa^2}{2^4\delta\omega^5(16+9\xi^2)}. \quad (56b)$$

Solving the system, it is possible to find the following expression for the FRC:

$$A_i^2 + A_r^2 = B_i^2 + B_r^2. \quad (57)$$

Substituting coefficients $A_r$, $A_i$, $B_r$ and $B_i$ into this expression yields a biquadratic equation in $\rho$. A peculiarity of the subharmonic resonance is to give rise to detached solution branches or isola, that are not connected to the main solution branch (Nayfeh and Mook 1995). Consequently, deriving the existence condition for such solutions is important in this context. The derivation of this condition from the biquadratic equation is reported in Appendix J. Only the result is shown here, the boundary for the region where solutions can exist being given by

$$\Omega = 3\delta\omega + \frac{63h\kappa^2}{2^7\delta\omega^5(16+9\xi^2)} + \frac{2^5\delta\omega\xi^2(16+9\xi^2)}{h\kappa^2}. \quad (58)$$

From Equation (57) it is also possible to derive an explicit expression for the FRC:

$$\frac{\Omega}{3} = \delta\omega + \frac{3h}{2^3\delta\omega}\rho^2 + \frac{3h\kappa^2}{2^4\delta\omega^5(16+9\xi^2)} \pm \sqrt{\frac{3^2h^2\kappa^2}{2^8\delta^2\omega^6(16+9\xi^2)}\rho^2 - \xi^2\omega^2}. \quad (59)$$

As in the superharmonic case, the terms outside the square root in Equation (59) can be interpreted to understand the shift of the primary resonance backbone curve for this subharmonic scenario, including the effect of the monomial $z_1z_3z_4$. Its expression is given by

$$\frac{\omega_{\text{NL}}}{3} = \delta\omega + \frac{3h}{2^3\delta\omega}\rho^2 + \frac{3h\kappa^2}{2^4\delta\omega^5(16+9\xi^2)}. \quad (60)$$

Also, by requiring Equation (59) to be single-valued, a generalised backbone curve that connects the minima of the FRCs' isolas can be obtained as:

$$\left\{ \begin{aligned} \rho(\kappa) &= \frac{2^4\omega^4\delta\xi\sqrt{16+9\xi^2}}{3h\kappa}, \\ \Omega(\kappa) &= 3\delta\omega + \frac{2^5\omega^7\delta\xi^2(16+9\xi^2)}{h\kappa^2} + \frac{3^2h\kappa^2}{2^4\delta\omega^5(16+9\xi^2)}. \end{aligned} \right. \quad (61)$$

The FRC can be further simplified by adopting a small damping hypothesis. Similarly to the previous superharmonic case, this assumption allows us to recover exactly the solution given by the first-order multiple scales development (Nayfeh and Mook 1995):

$$\frac{\Omega}{3} = \omega + \frac{3h}{2^3\omega}\rho^2 + \frac{3h\kappa^2}{2^8\omega^5} \pm \sqrt{\frac{3^2h^2\kappa^2}{2^{12}\delta^2\omega^6}\rho^2 - \xi^2\omega^2}. \quad (62)$$

The small damping assumption allows for an additional simplification in the expression of the boundary of the existence region, Equation (58), that also results in the same formula as the one found with a first-order multiple scales solution:

$$\Omega = 3\omega + \frac{63h\kappa^2}{2^{11}\omega^5} + \frac{2^9\omega^7\xi^2}{h\kappa^2}. \quad (63)$$

It should be noted that Equation (63) is quite close to the second expression in Equation (61) once it is simplified to retain only first-order terms in the damping:

$$\Omega = 3\omega + \frac{3^2h\kappa^2}{2^8\delta\omega^5} + \frac{2^9\omega^7\xi^2}{h\kappa^2}. \quad (64)$$

The two curves have different interpretations, however. While Equation (63) gives the possible combinations of parameters $\Omega$ and $\kappa$ such that non-trivial solutions might exist, Equation (64) only gives the frequency value at which the minimum of the FRC occurs as a function of $\kappa$. Both
<figure><img src="image_9.png" /><figcaption>**Figure 7** 1:3 subharmonic resonance of the damped cubic Duffing oscillator. Results obtained using a complex normal form style, with parameter values as $\omega = 1.5$, $h = 1$ and $\xi = 0.02$. (a) Comparison of the existence condition and the curve following the FRCs minima in $(\Omega, \kappa^2)$ space. (b) FRCs for $\kappa = 5, 10, 15$. The figure also shows the shifted backbones and generalised backbone joining the FRCs minima. In the legend, the grey lines indicate which equation has been used to calculate each curve, while different colors are associated with different $\kappa$ values. For $\kappa = 0$ there is no FRC, only the forcing-dependent backbone.</figcaption></figure>

expressions are plotted in Figure 7(a), with parameter values fixed as $\omega = 1.5$, $h = 1$ and $\xi = 0.02$. It is worth mentioning that the curve defining the FRC minima lies inside the existence region given by Equation (63), as should be expected.

Equations (59) to (62) are used to plot the FRCs in Figure 7(b). For the figure, the forcing amplitude assumed three values: $\kappa \in \{5, 10, 15\}$, each of which is associated with blue, dark green and light green curves, respectively, while the other parameters remained the same as in the previous plot. The figure also shows the shifted backbones given by Equation (60) for the three different forcing values; and for $\kappa = 0$. The generalised backbone that connects the minima of the FRCs, whose expression is found in Equation (61), is also shown, in red dash-dotted line.

Numerical results are also presented for this situation in Figure 8, where the FRCs in terms of physical displacement amplitude $u_{max}$ are plotted. Once again, the numerical continuation package Matcont (Dhooge et al. 2004) is employed to find the FRCs stemming from MORFE_Symbolic and that of the full system. Additionally, a first-order multiple scales solution is also presented. Its displacement time series is obtained from

$$u = \rho \cos \frac{1}{3}(\Omega t - \gamma) + \frac{\kappa}{\omega^2 - \Omega^2} \cos \Omega t, \qquad (65)$$

with $\gamma$ defined in (Nayfeh and Mook 1995), page 180. The difference between the multiple scales solution and the one from the $O(z^3, \epsilon^3)$ parametrisation is even more pronounced in this case, as compared to the superharmonic resonance scenario. Again, the nonlinear mappings arising from the parametrisation method prove essential to obtain accurate approximations of the full system. Comparing the curves from the Figure 8, it can be seen that convergence is more difficult to reach. However, an order 9 parametrisation seems to be an acceptable approximation to the reference curve in the displayed frequency range.

### 3.3 Parametric excitation

In this section, the case of a parametric excitation is considered. The starting point is a damped cubic Duffing equation with a forcing on the right-hand side that depends on the displacement, reading:

$$\ddot{u} + 2\xi\omega\dot{u} + \omega^2 u + hu^3 = \kappa u \cos \Omega t. \qquad (66)$$

There exist several parametric instabilities for different values of the excitation frequency $\Omega$, that are generally reported in a so-called Strutt-Ince stability diagram, see e.g. (Nayfeh and Mook 1995;
<figure><img src="image_10.png" /><figcaption>**Figure 8** FRCs for the 1:3 subharmonic resonance of the damped cubic Duffing oscillator. Parameter values are $\omega = 1.5$, $h = 1$, $\xi = 0.002$ and $\kappa = 1$. Comparison of numerically computed FRCs, obtained from MORFE_Symbolic with higher orders and solved by continuation, with a reference solution, also obtained by continuation, and with the first-order multiple scales solution (MMS).</figcaption></figure>

Grandi et al. 2021). Here, the analysis is restricted to the primary parametric resonance when the excitation frequency is close to twice the eigenfrequency, $\Omega \approx 2\omega$. The package MORFE_Symbolic relies on the generic treatment of non-autonomous terms proposed in (Vizzaccaro et al. 2024), which can handle a first-order differential-algebraic equation (DAE). Also, the starting point assumes that a quadratic recast of the equations of motion has been performed, by adding new variables to express the initial problem with only quadratic nonlinearities (Guillot et al. 2019). Equation (66) is thus rewritten as follows to be automatically treated by MORFE_Symbolic:

$$ \dot{u} = v, \qquad (67a) $$

$$ \dot{v} = -2\xi\omega v - \omega^2 u - hur_1 + ur_2, \qquad (67b) $$

$$ 0 = r_1 - u^2, \qquad (67c) $$

$$ 0 = r_2 - \frac{\kappa}{2}(z_3 + z_4), \qquad (67d) $$

$$ \dot{z}_3 = i\Omega z_3, \qquad (67e) $$

$$ \dot{z}_4 = -i\Omega z_4. \qquad (67f) $$

The complex normal form (CNF), up to the third order, is given in the equation that follows. The expansion point for the parametrisation is selected as $\Omega_p = 2\omega$.

$$ \dot{z}_1 = f_1 z_1 + f_2 z_1^2 z_2 + f_3 z_2 z_3 + f_4 z_1 z_3 z_4. \qquad (68) $$

Again, note that, as compared to the out-of-resonance scenario shown in Equation (34), the only added monomial to be considered for the parametric excitation is $z_2z_3$, which is indeed resonant since $\Omega \approx 2\omega$. The coefficients $f_j, j = 1, \dots, 4$ in Equation (68) are

$$ f_1 = \lambda_1 = -\xi\omega + i\omega\delta, \qquad f_2 = i\frac{3h}{2\delta\omega}, \qquad (69a) $$

$$ f_3 = -i\frac{\kappa}{2^2\delta\omega}, \qquad f_4 = i\frac{\kappa^2(\delta - 1)}{2^5\omega^3\xi^2\delta^2}. \qquad (69b) $$

Note in particular that in the present case, $f_2, f_3$ and $f_4$ are purely imaginary, which is a consequence of the starting point Equation (66). However, in a more general context of parametrically excited systems, these coefficients can have a non-vanishing real part. This is observed for example in the case of continuous structures where an external forcing leads to a parametric excitation, as the case studied for example in (Opreni et al. 2023) (beam with in-plane forcing where the parametric excitation leads to transverse vibrations), and in (Frangi et al. 2025) (parametric excitation due to electro-mechanical coupling). Hence, for the sake of generality, all the introduced coefficients are
considered to be complex, with $f_j = f_j^R + if_j^I$ to distinguish real and imaginary parts. This point is further addressed in Section 4.2 where a parametrically excited two-dofs system is considered. Note finally that $f_4$ scales with $\kappa^2$, it is thus a second-order term with respect to the forcing, which is neglected with first-order assumptions used in (Breunung and Haller 2018; Jain and Haller 2022; Opreni et al. 2023), and also neglected in first-order perturbative solutions, such as the MMS and the method of varying amplitude, as can be seen in (Benacchio et al. 2022).

Substituting $z_1 = \rho/2e^{i\theta}$ and $z_2 = \rho/2e^{-i\theta}$, defining $\psi = 2\theta - \Omega t = 2\theta - \phi$ to make the system autonomous, we obtain a first-order autonomous dynamical system

$$ \dot{\rho} = B_r + A_r \cos \psi + A_i \sin \psi, \quad (70a) $$

$$ \frac{\rho}{2} \dot{\psi} = B_i - A_r \sin \psi + A_i \cos \psi. \quad (70b) $$

where the introduced coefficients are written as a function of the real and imaginary parts of the $f_j$ (respectively denoted as $f_j^R$ and $f_j^I$) as:

$$ A_r = \rho f_3^R, \quad B_r = \rho f_1^R + \frac{\rho^3}{4} f_2^R + \rho f_4^R, \quad (71a) $$

$$ A_i = \rho f_3^I, \quad B_i = \rho f_1^I - \rho \frac{\Omega}{2} + \frac{\rho^3}{4} f_2^I + \rho f_4^I. \quad (71b) $$

Following the same procedure as in the previous superharmonic and subharmonic cases, an implicit expression for the branches of solutions (FRC) is

$$ A_r^2 + A_i^2 = B_r^2 + B_i^2. \quad (72) $$

Once the values of $A_r$, $A_i$, $B_r$ and $B_i$ are substituted, it is possible to factor out $\rho$ from all the terms, resulting in an equation of the form $\rho^2 f(\Omega, \rho) = 0$. Thus, $\rho = 0$ is always a fixed point of the system, which is a known result in parametrically excited systems (Nayfeh and Mook 1995). The non-zero solution for the amplitude $\rho$, that yields an explicit expression for the bifurcated solution branches, can be found by solving for $f(\Omega, \rho) = 0$, which gives

$$ \frac{\Omega^2}{4} - g^I(\rho)\Omega + |g(\rho)|^2 - |f_3|^2 = 0, \quad (73) $$

where the complex $g(\rho) = g^R(\rho) + ig^I(\rho)$ has been introduced as

$$ g(\rho) = f_1^R + \frac{\rho^2}{4} f_2^R + f_4^R + i \left( f_1^I + \frac{\rho^2}{4} f_2^I + f_4^I \right). \quad (74) $$

The solutions to Equation (73) are

$$ \Omega = 2 \left( f_1^I + \frac{\rho^2}{4} f_2^I + f_4^I \right) \pm 2 \sqrt{|f_3|^2 - \left( f_1^R + \frac{\rho^2}{4} f_2^R + f_4^R \right)^2}, \quad (75) $$

which is an explicit expression for the FRC of the parametric excitation. In the present case for the coefficients $f_j$ given in Equation (69), this equation is

$$ \Omega = 2\omega\delta \left( 1 + \frac{3h}{2^3\delta^2\omega^2}\rho^2 + \frac{\kappa^2(\delta - 1)}{2^5\omega^4\delta^3\xi^2} \right) \pm 2\sqrt{\frac{\kappa^2}{2^4\omega^2\delta^2} - \xi^2\omega^2}. \quad (76) $$

To simplify Equation (76), a first level of approximation can be found by neglecting high order terms on the forcing, i.e. the third term that scales with $\kappa^2$ in the first bracket disappears. This term comes from the monomial $z_1z_3z_4$ and is a second-order term for the forcing, which is neglected in usual perturbative solutions. With this assumption, the expression for the FRC simplifies to

$$ \Omega = 2\omega\delta \left( 1 + \frac{3h}{2^3\delta^2\omega^2}\rho^2 \right) \pm 2\sqrt{\frac{\kappa^2}{2^4\omega^2\delta^2} - \xi^2\omega^2}. \quad (77) $$
<figure><img src="image_11.png" /><figcaption>**Figure 9** Parametric excitation with parameter values set as $\omega = 1$ and $h = 1$, and a complex normal form analysis. (a) Marginal stability curves or Floquet tongues for $\xi = 0, 0.05, 0.1, 0.2$. (b) FRCs for different levels of approximation. The damping factor is set as $\xi = 0.2$ and the forcing amplitude as $\kappa = 1$.</figcaption></figure>

Following that, the solution given by a first-order multiple scales development can be retrieved by further imposing a small damping assumption, where only the leading order term in $\xi$ is kept and $\delta = 1$ is set. The FRC is then

$$ \Omega = 2\omega \left( 1 + \frac{3h}{2^3 \omega^2} \rho^2 \right) \pm 2 \sqrt{\frac{\kappa^2}{2^4 \omega^2} - \xi^2 \omega^2}. \quad (78) $$

The bifurcation points from which the non-zero solution emerges are found by letting $\rho = 0$ in Equation (76). Let us call $\Omega_{a,b}$ these two period-doubling points, where a Hopf bifurcation occurs since a fixed point loses stability and a limit cycle emerges. They read, for the simple case of Equation (66)

$$ \Omega_{a,b} = 2\omega\delta \left( 1 + \frac{\kappa^2(\delta - 1)}{2^5\omega^4\delta^3\xi^2} \right) \pm 2\sqrt{\frac{\kappa^2}{2^4\omega^2\delta^2} - \xi^2\omega^2}. \quad (79) $$

These two equations can be represented as a function of the forcing amplitude $\kappa$ and give the marginal stability curves (or Floquet tongues) where the parametric excitation gives rise to non-zero solutions. They are presented in Figure 9(a) for different values of $\xi$, and retrieve the usual result where the minimum of these tongues increases with the damping ratio. Interestingly, the curves are not symmetric, a feature already reported with MMS solutions, see e.g. (Thomsen 2003), such that their shapes in the low-frequency range lead to crossing points, see Figure 9(a). About the bifurcated branches, their stability can be analyzed through a linear stability analysis. The detailed calculations for this case are presented in Appendix K.

A further level of refinement of the solution can be obtained by pushing the normal form development up to order four. In this case, the expression for the FRC becomes more complicated, but has the same form as Equation (40). For the sake of conciseness, the explicit expressions are not reported here, but are compared to the solutions given by Equations (76) to (78) in Figure 9(b), where the stability calculations developed in Appendix K are also included by indicating the unstable sections of the FRC with dashed lines. In the Figure 9(b), the influence of different terms on the FRC expression can be appreciated. While passing from the multiple scales solution to the complete order three normal form development, the effect of not considering small damping ($\delta^2 \approx 1$) is responsible for displacing the bifurcation points, see Figure 9(b). Additionally, the effect of the $z_1z_3z_4$ monomial is given by a shift in the midpoint between the two bifurcation points, generated by the term proportional to $\kappa^2$ inside the bracket in Equation (78). Interestingly, the two bifurcated branches of the FRC do not close with a third-order approximation in the normal form. This can be easily seen from the equations since the term under the square root is constant, see also Figure 9(b). This result is however particular to the starting point used in this study and the fact that the $f_j$ coefficients reported in Equation (69) are purely imaginary. A case
<figure><img src="image_12.png" /><figcaption>**Figure 10** FRCs for the parametric resonance of the damped cubic, parametrically excited Duffing oscillator. Common parameters for both figures are $\omega = 1$ and $h = 1$. Comparison of numerically computed FRCs, obtained from MORFE_Symbolic, with a reference solution obtained by numerical continuation. (a) Backbones for small levels of forcing and damping, with $\xi = 0.01$ and $\kappa = 0.043$. (b) Backbones for higher levels of forcing and damping, with $\xi = 0.2$ and $\kappa = 0.9$.</figcaption></figure>

with non-vanishing real parts is shown in Section 4.2 to extend the analysis. Note that, in the presented first-order perturbative solution, the solution branches are closing by using the actual forcing frequency $\Omega$ in the equations, as done for example in (Thomsen 2003) for the MMS and in (Benacchio et al. 2022) for the method of varying amplitude, a particular feature of these methods. In the present formulation, where the expansion is centered around fixed frequency $\Omega_p$, a fourth-order normal form needs to be computed in order to recover the expected solution with closing bifurcated branches.

This resonance scenario is also studied numerically in Figure 10. The FRCs for the physical maximal displacement $u_{max}$ are obtained from the output of MORFE_Symbolic, which is then solved numerically with a continuation procedure using the package Matcont (Dhooge et al. 2004). They are compared to the full system solutions, taken as reference, and obtained by numerical continuation directly on the initial problem. Two sets of parameter values are selected to highlight different features of the solutions. In Figure 10(a), parameter values are set as $\omega = 1$, $h = 1$, $\xi = 0.02$ and $\kappa = 0.043$. In this case, the order 3 solution is omitted, since its bifurcated branches do not close. Full convergence is obtained for an order 11 development. In Figure 10(b), the values of both the forcing and damping are severely increased to $\kappa = 0.9$ and $\xi = 0.2$, while keeping $\omega = 1$ and $h = 1$. This choice has been selected in order to better underline their effect on the bifurcation point and illustrate how these particular features are retrieved by the CNF solution. The numerical results show that the prediction of the bifurcation point given by the CNF is confirmed. MMS solutions and order three developments completely overpredict the maximum amplitude, such that higher orders are mandatorily needed. Note, however, that, with this choice of parameters, the CNF does not converge to the exact solution. Even though all the features are correctly retrieved, an underestimation in the maximum amplitude persists, underlining that the validity limits of the method have been reached. Also, a small shift in the bifurcation points that could not be captured by the CNF is visible.

# 4 Theoretical results, high-order developments on two-dof systems

This section aims to illustrate, with simple examples involving a two-degrees-of-freedom system, how the presence of a slave mode can modify the dynamics of the master mode. Three different cases with a master reduced dynamics considering a single NNM are considered. First, a generic two-dof system without internal resonance is considered in order to highlight the effects of the linear and nonlinear characteristics of the slave mode on the dynamics of the master mode. Then, a parametrically excited two-dof system is considered. Finally, a mechanical system that serves as
an illustrative example is used to highlight the previous findings.

## 4.1 Case without internal resonance

A conservative two-dof system with quadratic and cubic nonlinearities is considered as the starting point:

$$ \ddot{u}_1 + \omega_1^2 u_1 + g_{11}^1 u_1^2 + g_{12}^1 u_1 u_2 + g_{22}^1 u_2^2 + h_{111}^1 u_1^3 + h_{112}^1 u_1^2 u_2 + h_{122}^1 u_1 u_2^2 + h_{222}^1 u_2^3 = 0, \quad (80a) $$

$$ \ddot{u}_2 + \omega_2^2 u_2 + g_{11}^2 u_1^2 + g_{12}^2 u_1 u_2 + g_{22}^2 u_2^2 + h_{111}^2 u_1^3 + h_{112}^2 u_1^2 u_2 + h_{122}^2 u_1 u_2^2 + h_{222}^2 u_2^3 = 0. \quad (80b) $$

Note that, when needed, modal damping of the form $2\xi_j\omega_j\dot{u}_j$, for $j = 1, 2$, can be added to take losses into account. Also, external forcing can be appended to the right-hand sides. Since the internal forces of mechanical systems are generally derived from a potential energy, some symmetry relationships exist between the quadratic and cubic coefficients, see e.g. Appendix A in (Touzé et al. 2021) for a general discussion. In the present case, the quadratic coefficients satisfy

$$ g_{12}^1 = 2g_{11}^2, \quad \text{and} \quad g_{12}^2 = 2g_{22}^1, \quad (81) $$

meaning that only four free coefficients are at hand for quadratic terms. Of the 8 possible cubic nonlinear coefficients, only 5 are free since the following three relationships are fulfilled:

$$ h_{112}^1 = 3h_{111}^2, \quad h_{122}^1 = h_{112}^2, \quad \text{and} \quad h_{122}^2 = 3h_{222}^1. \quad (82) $$

Assume that mode 1 is the master mode, and that no internal resonance exists between $\omega_1$ and $\omega_2$. In that case, the reduced dynamics contains a single pair of master coordinates, ($z_1, z_2$). Since the analytical expressions become too lengthy for an easy interpretation, we begin with the case of a system with only cubic nonlinearities, such that all $g_{ij}^p$ coefficients vanish, for $p, i, j = 1, 2$. The reduced dynamics, up to order 5 is then

$$ \dot{z}_1 = i\omega_1 z_1 + i \frac{3h_{111}^1}{2\omega_1} z_1^2 z_2 + 3i \left( \frac{(41\omega_1^2 - 5\omega_2^2)(h_{111}^1)^2}{\omega_1(9\omega_1^2 - \omega_2^2)(\omega_1^2 - \omega_2^2)} - \frac{17(h_{111}^1)^2}{16\omega_1^3} \right) z_1^3 z_2^2. \quad (83) $$

This result has been obtained with MORFE_Symbolic using the complex normal form style, and the assumptions given by Equation (82) have been taken into account. The nonlinear mappings that relate the original coordinates ($u_1, v_1, u_2, v_2$) to the normal complex master variables ($z_1, z_2$), are reported in Appendix L.

It is interesting to observe that only two cubic coefficients appear in the reduced dynamics: $h_{111}^1$ and $h_{111}^2$. While $h_{111}^1$ refers only to the cubic term in $u_1$ and logically appears in Equation (83), the other cubic coefficient which plays a role is $h_{111}^2$. Because of the symmetry relationship Equation (82), this coefficient is related to two different monomials: the invariant-breaking term $u_1^3$ in the second equation, and $u_1^2 u_2$ in the first equation. The invariant-breaking term is the only one (in such two-dofs scenario) responsible for the loss of invariance of the linear eigenspace. Consequently, only $h_{111}^2$ is responsible for the curvatures of the invariant manifold in phase space. This is also evident in the equations governing the geometry of the manifolds, which are contained in the nonlinear mappings reported in Appendix L. Notice that enforcing $h_{111}^2 = 0$ in Equation (L.1) leads to $u_2 = v_2 = 0$, showing that the manifold has no extra curvatures in the direction of the second mode. Finally, looking in detail into the analytical equations of the nonlinear mappings shown in Appendix L, one can see that these expressions only depend on $h_{111}^1, h_{111}^2$, and $h_{122}^1$. Again, it appears logical that $h_{122}^1$ is involved since it is linked in the first equation to the second trivially resonant monomial $u_1 u_2^2$ (but also to monomial $u_1^2 u_2$ in the second equation due to symmetry), following the terminology introduced in (Touzé et al. 2004; Touzé 2014; Touzé et al. 2021). Interestingly, the present calculations underline that up to order 5 this term has an effect only on the nonlinear mappings, but not on the reduced dynamics. However, once higher-order developments are pursued, the other missing coefficients start to appear, first in the nonlinear mappings and then in the reduced dynamics equations. As a final remark on higher orders, we highlight that the monomials present in the reduced dynamics are
the same as the ones examined for the Duffing equation, due to the simple rule of construction of the resonance. In the present case, only terms of the form $z_1^{p+1}z_2^p$, with $p \ge 1$, will stay in the reduced dynamics, while the coefficients are modified by the presence of the second oscillator.

Thanks to the property of the CNF already underlined in Section 2.1, the backbone curve is analytic up to the desired order. Since the coefficients are lengthy for orders higher than 5, only the fifth-order is shown here, but of course, closed-form expressions are readily available for high-order terms:

$$ \frac{\omega_{\text{NL}}}{\omega_1} = 1 + \frac{3}{8} \frac{h_{111}^1}{\omega_1^2} \rho^2 + \frac{3}{16} \left( \frac{41\omega_1^2 - 5\omega_2^2}{(9\omega_1^2 - \omega_2^2)(\omega_1^2 - \omega_2^2)} \left( \frac{h_{111}^2}{\omega_1^2} \right)^2 - \frac{17}{16} \left( \frac{h_{111}^1}{\omega_1^2} \right)^2 \right) \rho^4. \quad (84) $$

Interestingly, the coefficient $h_{111}^2$ which appears in the cubic term might have an effect on the hardening/softening behaviour at large amplitude. If the first-order term is unequivocally dictated by the value of $h_{111}^1$, one can see that the next order depends in a complicated manner also on $h_{111}^2$, $\omega_1$ and $\omega_2$. This is illustrated in Figure 11(a) where three different scenarios have been tested. In each case, a high-order backbone curve is computed from the expressions obtained with the parametrisation method up to order 25, and is compared to a reference solution obtained by a numerical continuation technique.

<figure><img src="image_13.png" /><figcaption>Figure 11 Backbones up to order 25 for the 2-dofs system without internal resonance. The curves were obtained using a complex normal form style, and parameter values were fixed as $\omega_1 = 1$, $h_{111}^1 = 1$, $h_{111}^2$, $g_{11}^1$ and $g_{11}^2$ varying, and all other nonlinear stiffness coefficients zero (other than the ones fixed by symmetry conditions). The backbones given by the parametrisation method are compared with reference solutions obtained using the numerical continuation package Matcont (Dhooge et al. 2004). (a) Oscillator with only cubic coefficients. (b) Oscillator with quadratic and cubic coefficients. Additional parameters fixed: $\omega_2 = 1.57$, $h_{111}^2 = 1$.</figcaption></figure>

The first case is where the invariant-breaking term vanishes, i.e. $h_{111}^2 = 0$. In such a case, the backbone follows that of a single-dof Duffing equation. Then two other cases are selected, using $h_{111}^2 = 1$, $\omega_1 = 1$, and either $\omega_2 = 1.57$, or $\omega_2 = 0.637$. The slave mode's eigenfrequencies were purposefully selected such that no low-order resonance relationships are verified. Figure 11(a) shows the obtained results. When $\omega_2 = 1.57$, some terms in the quintic coefficient shown in Equation (84) are negative, such that the hardening behaviour turns back to softening at higher amplitudes. This effect is correctly captured by the asymptotic solution to order 25 which is close to the numerical solution up to amplitudes around 0.65. In this case, the validity limit of the
asymptotic approach seems to decrease as compared to the case $h_{111}^2 = 0$. On the other hand, when the slave mode has the lowest eigenfrequency with $\omega_2 = 0.637$, the hardening behaviour is enhanced by the higher-order terms. In this case, the validity limit seems to be even smaller since the order 25 solution departs from the reference from amplitudes around 0.45.

We now consider the effect of the quadratic terms by considering all quadratic and cubic coefficients in Equation (80). Since the full analytical expressions begin to be lengthy, only the third-order reduced dynamics with CNF is reported to see how the quadratic coefficients aggregate to form the third-order coefficient in the normal form that dictates the hardening/softening behaviour. It reads:

$$ \dot{z}_1 = i\omega_1 z_1 + i\omega_1 \left( \frac{3h_{111}^1}{2\omega_1^2} - \frac{5}{3}\left(\frac{g_{11}^1}{\omega_1^2}\right)^2 - \frac{8\omega_1^2 - 3\omega_2^2}{4\omega_1^2 - \omega_2^2}\left(\frac{g_{11}^2}{\omega_1\omega_2}\right)^2 \right) z_1^2 z_2. \quad (85) $$

Notice that, at this order, two quadratic coefficients are of special interest. First, the “self-quadratic” one, $g_{11}^1$, which is present in a single-dof Duffing analysis, see Section 3.1. Second, the invariant-breaking quadratic coefficient $g_{11}^2$. Also note that the denominator of the invariant breaking term vanishes when a 2:1 internal resonance such that $\omega_2 = 2\omega_1$ is met, a feature already discussed in (Touzé et al. 2004; Touzé and Thomas 2006).

Figure 11(b) illustrates this effect of the quadratic terms by considering the increase of both $g_{11}^1$ and $g_{11}^2$ with the same trend, by imposing $g_{11}^1 = g_{11}^2$, and increasing them from 0.5 to 1.5. As expected, the quadratic nonlinearity favours the softening behaviour, and one can see that the types of behaviour ranged by this choice are various and generally well reproduced by the normal form analysis, up to a certain amplitude that should correspond to the validity limit of the asymptotic development.

## 4.2 Parametrically excited system

This section considers a simplified two-dofs system which is representative of a flat structure that is externally excited with an in-plane force, which leads to transverse vibrations through a parametric excitation. In such a case, studied numerically with a cantilever beam in (Opreni et al. 2023), the external forcing transforms to a parametric excitation in the reduced dynamics. The simplest system that can reproduce such an effect needs to contain a quadratic coupling between master and slave modes, as is the case in beam structures, for example, (Vizzaccaro et al. 2020; Givois et al. 2019). The model introduced in (Vizzaccaro, Salles, et al. 2021) is thus selected here, i.e.

$$ \ddot{u}_1 + 2\xi\dot{u}_1 + u_1 + 2gu_1u_2 + hu_1^3 = 0, \quad (86a) $$

$$ \ddot{u}_2 + 2\xi\omega_2\dot{u}_2 + \omega_2^2 u_2 + gu_1^2 = \kappa \cos \Omega t, \quad (86b) $$

where the forcing frequency $\Omega$ is selected in the vicinity of twice the eigenfrequency of the master mode (normalized at 1 here), $\Omega \approx 2$, in order to activate the principal parametric resonance. Additionally, a small damping assumption is introduced to simplify the coefficients. In this case, the reduced dynamics up to the third order, with a parametrisation computed with the choice $\Omega_p = 2$, is

$$ \dot{z}_1 = f_1 z_1 + f_2 z_1^2 z_2 + f_3 z_2 z_3 + f_4 z_1 z_3 z_4, \quad (87) $$

with the coefficients

$$ \begin{aligned} f_1 &= i - \xi, & f_3 &= \frac{g\kappa(i - \xi)}{2(i\xi + 1)(\omega_2 - 2)(4i\xi + \omega_2 + 2)}, & f_4 &= -\frac{g^2\kappa^2(2i - 7\xi)}{32(\omega_2^2 - 4)^2} \\ f_2 &= \frac{(i - \xi)(-16g^2i\xi\omega_2 + 32g^2i\xi - 6g^2\omega_2^2 + 16g^2 - 12hi\xi\omega_2^3 + 24hi\xi\omega_2^2 - 3h\omega_2^4 + 12h\omega_2^2)}{2\omega_2^2(i\xi + 1)(\omega_2 - 2)(4i\xi + \omega_2 + 2)}. \end{aligned} \quad (88) $$

The reduced dynamics equation is, as expected, the same as the one shown in Equation (68). In this situation, however, the coefficients $f_j$, and in particular $f_2$, have non-vanishing real parts.
<figure><img src="image_14.png" /><figcaption>**Figure 12** Externally forced 2-dofs system leading to a parametric resonance, Equation (86). Comparison of the results obtained with a third-order complex normal form style and the method of multiple scales (MMS). Parameter values set as $\omega_2 = 1.57$, $\xi = 0.2$, $h = 0.5$, $g = 5$ and $\kappa = 0.2$.</figcaption></figure>

As a consequence, the term under the square root in Equation (75) depends on $\rho$ and the FRC closes at order three. This situation is illustrated in Figure 12.

From the figure, note that for the selected parameter values the bifurcated branches of the FRC display a softening behaviour. Inspecting Equation (75) shows that the hardening/softening behaviour is governed by the imaginary part $f_2^I$ of the coefficient $f_2$. In the present case, it can be written explicitly as:

$$f_2^I = g^2 \frac{3\omega_2^2 - 8}{\omega_2^2(4 - \omega_2^2)} - \frac{3h}{2}. \qquad (89)$$

Consequently, the hardening/softening transition occurs for this specific case when $f_2^I$ vanishes, and a hardening behaviour is obtained for $h < g^2(6\omega_2^2 - 16)/(3\omega_2^2(4 - \omega_2^2))$. Again, the wealth of the proposed ROMs can be here underlined since they provide physical insights and predictive interpretations that can be directly checked on the numerical results. Otherwise, the same comments made in Section 3.3 also apply here. The slight shift of the bifurcation points between the multiple scales solution and the order 3 normal form development is also retrieved and is illustrated in Figure 12. A stability analysis similar to the one given in Appendix K could also be applied here, but is not done for concision, and therefore not reported in the figure.

## 4.3 Illustrative example: a mass connected to two nonlinear springs

This last section aims at illustrating some of the results presented in the paper to the case of a 2 dof system consisting of a mass connected to two elastic nonlinear springs. This example has been introduced in (Touzé et al. 2004) and then used in several articles as a benchmark study (Touzé and Amabili 2006; C.-H. Lamarque et al. 2012; Breunung and Haller 2018; Liu and Wagg 2019). The equations of motion read:

$$\ddot{u}_1 + \omega_1^2 u_1 + \frac{\omega_1^2}{2}(3u_1^2 + u_2^2) + \omega_2^2 u_1 u_2 + \frac{\omega_1^2 + \omega_2^2}{2} u_1 (u_1^2 + u_2^2) = 0 \qquad (90a)$$

$$\ddot{u}_2 + \omega_2^2 u_2 + \frac{\omega_2^2}{2}(3u_2^2 + u_1^2) + \omega_1^2 u_1 u_2 + \frac{\omega_1^2 + \omega_2^2}{2} u_2 (u_1^2 + u_2^2) = 0. \qquad (90b)$$

It should be highlighted that the nonlinear stiffness coefficients can be written as a function of the eigenfrequencies of the system, leading to simplified expressions for the coefficients that appear in the high-order normal forms. Up to order 5, using again the CNF and showing only the equation for $z_1$, we obtain

$$\dot{z}_1 = i\omega_1 z_1 + i \frac{4\omega_1 (-3\omega_1^2 + \omega_2^2)}{4\omega_1^2 - \omega_2^2} z_1^2$$
$$ + i \frac{\omega_1 (-9072\omega_1^{10} + 27180\omega_1^8\omega_2^2 - 19624\omega_1^6\omega_2^4 + 5835\omega_1^4\omega_2^6 - 754\omega_1^2\omega_2^8 + 35\omega_2^{10})}{576\omega_1^{10} - 1072\omega_1^8\omega_2^2 + 652\omega_1^6\omega_2^4 - 177\omega_1^4\omega_2^6 + 22\omega_1^2\omega_2^8 - \omega_2^{10}} z_1^3 z_2^2. \quad (91) $$

Furthermore, using polar realification, an analytic expression for the backbone can be found as

$$ \omega_{\text{NL}} = \omega_1 + \frac{-3\omega_1^3 + \omega_1\omega_2^2}{4\omega_1^2 - \omega_2^2}\rho_1^2 + \frac{-9072\omega_1^{11} + 27180\omega_1^9\omega_2^2 - 19624\omega_1^7\omega_2^4 + 5835\omega_1^5\omega_2^6 - 754\omega_1^3\omega_2^8 + 35\omega_1\omega_2^{10}}{16(4\omega_1^2 - \omega_2^2)^3(9\omega_1^2 - \omega_2^2)(\omega_1^2 - \omega_2^2)}\rho_1^4. \quad (92) $$

We denote as $f^{(p)}$ the coefficient of the monomial of order $p$ in the normal form. The third-order term, which governs the first curvature of the backbone curve and the hardening/softening behaviour, has already been analysed in (Touzé et al. 2004) with the oscillator normal form, leading to the same coefficient as $f^{(3)}$ reported here with the CNF. Thanks to the high-order expansions, the analysis can be pushed further by considering the variations of $f^{(5)}$, whose expression is stated in Equation (91), and $f^{(7)}$, not shown here for the sake of brevity, whose variations are represented in Figure 13 for $\omega_1 = 1$ and varying $\omega_2$.

<figure><img src="image_15.png" /><figcaption>**Figure 13** Behaviour of the normal form coefficients $f^{(3)}$, $f^{(5)}$ and $f^{(7)}$ for the two-dofs system given in Equation (90) as a function of $\omega_2$, with $\omega_1 = 1$.</figcaption></figure>

From the figure, it is straightforward to predict the type of nonlinearity for the order 7 backbone curve. Taking for instance $\omega_2 = 1.5$, the system shows a softening behaviour, since all three coefficients $f^{(3)}$, $f^{(5)}$, and $f^{(7)}$, are negative. On the contrary, for $\omega_2 = 2.5$, since $f^{(3)}$ and $f^{(5)}$ are positive, the backbone will first display a hardening behaviour, turning to softening for large amplitudes since $f^{(7)}$ is highly negative and will dominate.

Figure 13 also illustrates that the coefficients have a singular behaviour when internal resonances are crossed, a typical feature resulting from the small denominator problem. The third-order coefficient $f^{(3)}$ displays a singularity when a 1 : 2 internal resonance $\omega_2 = 2\omega_1$ appears, a behaviour already analysed e.g. in (Rega et al. 2000; Touzé et al. 2004; Touzé and Thomas 2006; Arafat and Nayfeh 2003). Interestingly, the order five coefficient $f^{(5)}$ shows in addition a singularity for the first third-order resonances, namely the 1 : 1 and the 1 : 3 resonances defined by $\omega_2 = \omega_1$ and $\omega_2 = 3\omega_1$. This is a direct consequence that the order five coefficient is built from the elimination of the cubic terms. Continuing further, the order seven coefficient $f^{(7)}$ has an additional singularity at the 1 : 4 resonance $\omega_2 = 4\omega_1$, following the cancellation of the quartic terms in the normal form. Finally, it is also possible to notice that for $\omega_2 > 4$ the values of the coefficients remain virtually constant, and the backbone behaviour is hardening at low amplitudes followed by softening at high amplitudes.

Whereas most of the reported analyses have been focused on FRCs, this section closes with numerous illustrative examples highlighting the behaviour of the invariant manifold serving as a reduced-order subspace in the method. First, the convergence in terms of the geometry and the effect of higher-order terms on the computed curvatures are illustrated.
<figure><img src="image_16.png" /><figcaption>**Figure 14** Phase space representation for the two-dofs system given by Equation (90). (a) Comparison of periodic orbits calculated with an order 20 parametrisation and the invariant manifold obtained by numerical continuation for $\omega_1 = 1$ and $\omega_2 = 1.57$. (b) Comparison of the manifolds of order 3 and 15, in blue and orange respectively, obtained by the parametrisation method with the one found by continuation, in pink, for $\omega_1 = \sqrt{0.5}$ and $\omega_2 = \sqrt{6}$. The black point represents the origin. (c-d) Cross-sections at $v_1 = 0.1$ and $u_1 = -0.2$, respectively, obtained from case (b). Solutions of order 3, 5, 7 and 15 are compared to the numerical continuation reference.</figcaption></figure>

Figure 14(a) compares the exact shape of the invariant manifold, which has been computed by numerical continuation, to periodic orbits calculated using an order 20 parametrisation, with $\omega_1 = 1$ and $\omega_2 = 1.57$. It can be seen that the orbits show an almost perfect agreement with the manifold in a region surrounding the fixed point located at the origin. In the same spirit, Figure 14(b) compares the numerical manifold obtained by continuation for $\omega_1 = \sqrt{0.5}$ and $\omega_2 = \sqrt{6}$, depicted in blue, to different orders of the asymptotic expansions: order 3 in orange, and order 15, in pink. It clearly highlights that, with increasing orders, the asymptotically computed manifolds can retrieve the complex curvatures of the exact solution. To aid the visualisation, section cuts of Figure 14(b) are presented in Figures 14(c) and 14(d), where orders 5 and 7 solutions are also included to better appreciate the convergence. Note that Figure 14(d) shows a cross-section at $u_1 = -0.2$ which is quite far from the origin, hence explaining the observed discrepancies.

The next case illustrates how the linear viscous damping terms might affect the shape of the invariant manifolds, thus highlighting the difference between a conservative and a damped manifold. Viscous damping of the form $2\xi_i\omega_i$, $i \in \{1, 2\}$, is appended to Equation (90). Figure 15 shows the invariant manifolds computed by the parametrisation method with CNF, with system's parameters fixed as $\omega_1 = \sqrt{0.5}$, $\omega_2 = \sqrt{6}$ and $\xi_1 = \xi_2 = 0.1$. In Figure 15(a) the damped manifold, in blue, is compared with the one obtained without damping, in yellow. It is possible to notice
<figure><img src="image_17.png" /><figcaption>Figure 15 Influence of damping on the invariant manifolds for the mass connected to two springs. (a) Comparison of the order 15 damped (blue) and undamped (yellow) manifolds obtained by the parametrisation method with CNF style, for $\omega_1 = \sqrt{0.5}$, $\omega_2 = \sqrt{6}$ and $\xi_1 = \xi_2 = 0.1$. Two periodic orbits, in red, are also shown for the full conservative system as well as a damped orbit of the full system, in blue. The black point represents the origin. (b-c) Cuts at $v_1 = 0$ and $u_1 = -0.2$, respectively, obtained from (a).</figcaption></figure>

that in the dissipative case, the curvatures of the manifold are importantly attenuated, which can be better appreciated by inspecting Figures 15(b) to 15(c), which show cross-sections of the manifolds.

Lastly, the effect of the external forcing on the manifolds is examined in Figure 16, to illustrate the behaviour of time-dependent invariant manifolds in the context of model order reduction. Two different orders of truncations relative to the forcing term are considered. A first case where a first-order assumption on the forcing is adopted, following (Opreni et al. 2023; Jain and Haller 2022) (first column, Figures 16(a) to 16(e), $O(z^3, \epsilon^1)$ truncation). In the second case, an $O(z^3, \epsilon^3)$ truncation is shown, highlighting the higher-order effects of the non-autonomous terms computed thanks to the development shown in (Vizzaccaro et al. 2024) (second column, Figures 16(b) to 16(f)).

Additionally, three cases relative to the forcing configuration are examined. In these situations, a non-autonomous right-hand side is added to the equations of motion:

$$ \ddot{u}_1 + 2\xi_1\omega_1\dot{u}_1 + \omega_1^2 u_1 + \frac{\omega_1^2}{2}(3u_1^2 + u_2^2) + \omega_2^2 u_1 u_2 + \frac{\omega_1^2 + \omega_2^2}{2}u_1(u_1^2 + u_2^2) = \kappa_1 \cos \Omega t \quad (93a) $$

$$ \ddot{u}_2 + 2\xi_2\omega_2\dot{u}_2 + \omega_2^2 u_2 + \frac{\omega_2^2}{2}(3u_2^2 + u_1^2) + \omega_1^2 u_1 u_2 + \frac{\omega_1^2 + \omega_2^2}{2}u_2(u_1^2 + u_2^2) = \kappa_2 \cos \Omega t. \quad (93b) $$

In all scenarios, the forcing is considered resonant to the first (master) mode, such that $\Omega \approx \omega_1$.
<figure><img src="image_18.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(f)</figcaption></figure>

<caption>**Figure 16** Influence of forcing on the invariant manifolds for the mass connected to two springs. In all figures, parameter values are fixed as $\omega_1 = 1$, $\omega_2 = 1.57$ and $\xi_1 = \xi_2 = 0.01$, and a primary resonance situation is considered. For figures (a), (c) and (e), $\kappa = 1$ and an $O(z^3, \epsilon^1)$ truncation is chosen, while for figures (b), (d) and (f), $\kappa = 0.2$ and an $O(z^3, \epsilon^3)$ is adopted. The path of the fixed point in phase space over time is also represented: (a-b) Forcing perfectly aligned with the master mode. (c-d) Forcing perfectly orthogonal to the master mode. (e-f) Forcing in both master and slave modes.</caption>
In the first of them, corresponding to Figures 16(a) to 16(b), the forcing is perfectly aligned with the master eigenspace. This case is obtained by setting a harmonic forcing only on Equation (93a), i.e. $κ₁ = κ, κ₂ = 0$. In this situation, when the development of the nonautonomous terms stops at order 1, the manifolds exhibit a rigid body rotation along a path that remains on the master eigenplane, as shown in (Opreni et al. 2023) and illustrated in Figure 16(a). Once high-order forcing terms are considered, however, the rotation of the fixed point shows slight deviations from the linear eigenplane due to the higher-order effects. Additionally, the motion of the manifold is no longer a rigid body one, and shows important deformations, as can be seen in Figure 16(b).

The second case corresponds to a forcing that is now orthogonal to the master eigenspace, obtained by setting a harmonic forcing only on Equation (93b), i.e. $κ₁ = 0, κ₂ = κ$. The first-order development for the non-autonomous part predicts a rigid body motion of the manifold along a circular motion on the slave eigenplane (Opreni et al. 2023), which when projected into the adopted representation space corresponds to a straight vertical line, as seen in Figure 16(c). With the inclusion of the high-order forcing terms once again the manifold deforms, and the motion path along the phase space is no longer circular. Its projection to the representation space, nevertheless, still remains a straight line, as can be seen in Figure 16(d), which indicates that it is entirely contained in the slave eigenspace.

Finally, Figures 16(e) to 16(f) showcase the general case where the forcing is neither orthogonal nor parallel to the master eigenspace. In this situation we chose $κ₁ = κ₂ = κ$ and, as can be observed from the figures, the motion of the manifolds is a combination of the two previous situations.

# 5 Conclusion

High-order, automated solutions for nonlinear vibrations, have been derived and analysed thanks to a symbolic version of the parametrisation method for invariant manifold, implemented in the package `MORFE_Symbolic`. The contribution aims to show how this general method can be used efficiently to derive several useful results to understand the reduced dynamics of large dimensional systems. It also underlines the existing continuity between low-order analytical perturbative solutions and high-order numerical solutions. When the solutions are analytically tractable, we have shown how they can easily reduce to known results obtained with perturbative methods, and how they extend by considering more effects with fewer assumptions. Indeed, a key feature of the parametrisation method is to offer high-order expansions with a single assumption on the smallness of the amplitudes. This is in contrast with perturbative solutions where a scaling of the different terms (nonlinearity, forcing, damping) needs to be introduced with an $ε$ parameter at the first step. This is not needed here and comes automatically from the analysis.

Most of the analyses we perform here use the complex normal form (CNF), after having introduced the differences between three variants of normal forms that have been derived in the nonlinear vibration literature. The symbolic package also allows derivations of solutions with the so-called graph style, a feature that has not been investigated here but which can be automatically computed to compare with normal form style. Analysis on primary, secondary and parametric resonance have been exhibited, and the scope has been limited to cases where the master dynamics contains a single master mode. We highlight that the results can be extended and take into account many more scenarios.

In the course of the developments, the wealth of the computed solutions has been regularly commented on, showing how they generalize previously known results limited to first orders, and how they can converge to the exact solution by considering the effects of the higher orders. We believe that systematically using the complex normal form, in combination with the parametrisation method, defines a new powerful grammar for solving nonlinear vibration problems, from low-order analytical solutions to high-order numerical solutions. Besides, normal forms can be used for system identification since providing the skeleton of the dynamics and the simplest dynamical system reproducing a given resonance scenario.

The only limitation of the approach is the validity limit in terms of amplitude, which has been commented on throughout the text while analysing the results. Future work will thus consider extending the presented results to other resonance scenarios, to build a dictionary of reference
solutions. As an example, general solutions for 1 : 2 internal resonance could be derived to generalise the results shown in (Gobat et al. 2021), and resonant phase lags could also be obtained with high-order solutions to generalise the results shown in (Volvert and Kerschen 2022). Finally, the question of accurately estimating the validity limits of the method in terms of amplitudes should be tackled to provide the analyst with a priori and a posteriori amplitude bounds.

## A Short presentation of MORFE_Symbolic.jl

This appendix gives a brief introduction to the use and output possibilities of MORFE_Symbolic. The code is available as a repository on version control platform github, under the MORFE project, and is accessible via the link github.com/MORFEproject/MORFE_Symbolic. Instructions on how to download and run the code are given there. Additionally, the code presents descriptions to all its functions, making their inputs and outputs explicit for the interested user.

MORFE_Symbolic treats dynamical systems written in first-order format, and deals only with quadratic nonlinearities. In such a way, quadratic recast (Guillot et al. 2019) has to be applied as an intermediate step in case arbitrary nonlinear functions should be considered, such as is exemplified in Section 3.3 for the parametric excitation case. At present, the code only supports cases with mono-harmonic excitation.

Since most of the examples of interest to the authors are mechanical systems, auxiliary functions allowing for an easier treatment in such cases, i.e. the construction of the system of equations from stiffness, mass and damping matrices and the quick assembly of such matrices for specific cases, have been devised. The reader is referred to the code repository, where the examples treated in this contribution are available, in the folder test, and showcase all of the present code capabilities.

Other than defining the equations of motion, inputs concerning the parametrisation method also need to be given to the code, such as the number of autonomous and non-autonomous variables considered in the parametrisation (i.e. the number of master modes), the maximum expansion order and a symbolic list of algebraic relationships between the master modes' eigenvalues, used to determine the resonance relationships. As an example, consider the case of an autonomous lightly damped mechanical system with two master modes in a 1 : 2 internal resonance. If the eigenvalues are {$\lambda_1, \lambda_2, \lambda_3, \lambda_4$}, and such that $\lambda_1$ and $\lambda_2$ refer to the first mode and $\lambda_3$ and $\lambda_4$ to the second mode, the relationships $\lambda_1 \approx -\lambda_2, \lambda_3 \approx -\lambda_4$ and $\lambda_1 \approx 2\lambda_3$ need to be given in a specific list of resonance conditions as an input.

The code also presents functions capable of performing polar realification, restricted to some specific cases where it is advantageous, and cartesian realification. Furthermore, dedicated functions to the output of the reduced dynamics, nonlinear mappings, backbone curve, nonlinear damping and physical amplitudes in LaTeX form are available. For cases where more powerful symbolic capabilities than those of the Julia language are required, a function that writes output files for use in Mathematica environment is also present, making it possible to export all of the necessary variables for post-treatment in the software.

## B Nonlinear change of coordinates for the CNF of the Duffing equation

In this appendix, the nonlinear mappings used to compute the complex normal form of the conservative Duffing equation with only cubic nonlinearity, are given, for the sake of completeness. The normal form up to order 11 is given in the main text as Equation (8). The nonlinear mapping allowing one to get from Equation (3) to the normal form, up to order order 11, reads:

$$ u = (z_1 + z_2) + \left( \frac{h}{8\omega^2} z_1^3 - \frac{3h}{4\omega^2} z_1^2 z_2 - \frac{3h}{4\omega^2} z_1 z_2^2 + \frac{h}{8\omega^2} z_2^3 \right) + \left( \frac{h^2}{64\omega^4} z_1^5 - \frac{39h^2}{64\omega^4} z_1^4 z_2 + \frac{69h^2}{32\omega^4} z_1^3 z_2^2 \right. \\ \left. + \frac{69h^2}{32\omega^4} z_1^2 z_2^3 - \frac{39h^2}{64\omega^4} z_1 z_2^4 + \frac{h^2}{64\omega^4} z_2^5 \right) + \left( \frac{h^3}{512\omega^6} z_1^7 - \frac{73h^3}{512\omega^6} z_1^6 z_2 + \frac{1569h^3}{512\omega^6} z_1^5 z_2^2 \right. \\ \left. - \frac{2139h^3}{256\omega^6} z_1^4 z_2^3 - \frac{2139h^3}{256\omega^6} z_1^3 z_2^4 + \frac{1569h^3}{512\omega^6} z_1^2 z_2^5 - \frac{73h^3}{512\omega^6} z_1 z_2^6 + \frac{h^3}{512\omega^6} z_2^7 \right) + \left( \frac{h^4}{4096\omega^8} z_1^9 \right) $$
$$ 
\begin{aligned}
& -\frac{107h^4}{4096\omega^8}z_1^8z_2 + \frac{1049h^4}{1024\omega^8}z_1^7z_2^2 - \frac{65379h^4}{4096\omega^8}z_1^6z_2^3 + \frac{75891h^4}{2048\omega^8}z_1^5z_2^4 + \frac{75891h^4}{2048\omega^8}z_1^4z_2^5 \\
& -\frac{65379h^4}{4096\omega^8}z_1^3z_2^6 + \frac{1049h^4}{1024\omega^8}z_1^2z_2^7 - \frac{107h^4}{4096\omega^8}z_1z_2^8 + \frac{h^4}{4096\omega^8}z_2^9 \\
& + \left( \frac{h^5}{32768\omega^{10}}z_1^{11} - \frac{141h^5}{32768\omega^{10}}z_1^{10}z_2 + \frac{7979h^5}{32768\omega^{10}}z_1^9z_2^2 - \frac{222369h^5}{32768\omega^{10}}z_1^8z_2^3 + \frac{2798925h^5}{32768\omega^{10}}z_1^7z_2^4 - \frac{727419h^5}{4096\omega^{10}}z_1^6z_2^5 \\
& -\frac{727419h^5}{4096\omega^{10}}z_1^5z_2^6 + \frac{2798925h^5}{32768\omega^{10}}z_1^4z_2^7 - \frac{222369h^5}{32768\omega^{10}}z_1^3z_2^8 + \frac{7979h^5}{32768\omega^{10}}z_1^2z_2^9 - \frac{141h^5}{32768\omega^{10}}z_1z_2^{10} \\
& +\frac{h^5}{32768\omega^{10}}z_2^{11} \right)
\end{aligned}
\quad (B.1a) $$

$$ 
\begin{aligned}
v = (i\omega z_1 - i\omega z_2) + \left( i\frac{3h}{8\omega}z_1^3 + i\frac{3h}{4\omega}z_1^2z_2 - i\frac{3h}{4\omega}z_1z_2^2 - i\frac{3h}{8\omega}z_2^3 \right) + \left( i\frac{5h^2}{64\omega^3}z_1^5 - i\frac{81h^2}{64\omega^3}z_1^4z_2 \right. \\
\left. -i\frac{69h^2}{32\omega^3}z_1^3z_2^2 + i\frac{69h^2}{32\omega^3}z_1^2z_2^3 + i\frac{81h^2}{64\omega^3}z_1^2z_2^4 - i\frac{5h^2}{64\omega^3}z_2^5 \right) + \left( i\frac{7h^3}{512\omega^5}z_1^7 - i\frac{305h^3}{512\omega^5}z_1^6z_2 \right. \\
\left. + i\frac{2691h^3}{512\omega^5}z_1^5z_2^2 + i\frac{2139h^3}{256\omega^5}z_1^4z_2^3 - i\frac{2139h^3}{256\omega^5}z_1^3z_2^4 - i\frac{2691h^3}{512\omega^5}z_1^2z_2^5 + i\frac{305h^3}{512\omega^5}z_1z_2^6 \right. \\
\left. -i\frac{7h^3}{512\omega^5}z_2^7 \right) + \left( i\frac{9h^4}{4096\omega^7}z_1^9 - i\frac{665h^4}{4096\omega^7}z_1^8z_2 + i\frac{3895h^4}{1024\omega^7}z_1^7z_2^2 - i\frac{98757h^4}{4096\omega^7}z_1^6z_2^3 \right. \\
\left. -i\frac{75891h^4}{2048\omega^7}z_1^5z_2^4 + i\frac{75891h^4}{2048\omega^7}z_1^4z_2^5 + i\frac{98757h^4}{4096\omega^7}z_1^3z_2^6 - i\frac{3895h^4}{1024\omega^7}z_1^2z_2^7 + i\frac{665h^4}{4096\omega^7}z_1z_2^8 \right. \\
\left. -i\frac{9h^4}{4096\omega^7}z_2^9 \right) + \left( i\frac{11h^5}{32768\omega^9}z_1^{11} - i\frac{1161h^5}{32768\omega^9}z_1^{10}z_2 + i\frac{45437h^5}{32768\omega^9}z_1^9z_2^2 - i\frac{757245h^5}{32768\omega^9}z_1^8z_2^3 \right. \\
\left. + i\frac{3848751h^5}{32768\omega^9}z_1^7z_2^4 + i\frac{727419h^5}{4096\omega^9}z_1^6z_2^5 - i\frac{727419h^5}{4096\omega^9}z_1^5z_2^6 - i\frac{3848751h^5}{32768\omega^9}z_1^4z_2^7 + i\frac{757245h^5}{32768\omega^9}z_1^3z_2^8 \right. \\
\left. -i\frac{45437h^5}{32768\omega^9}z_1^2z_2^9 + i\frac{1161h^5}{32768\omega^9}z_1z_2^{10} - i\frac{11h^5}{32768\omega^9}z_2^{11} \right)
\end{aligned}
\quad (B.1b) $$

Note that the relationship $\dot{u} = v$ is preserved by the nonlinear change of coordinates. Indeed, taking the derivative of Equation (B.1a) with respect to time and eliminating the time dependence with the reduced dynamics given by Equation (8), yields Equation (B.1b). This holds in general and can be verified in each of the treated examples, independent of the normal form style. It is a direct consequence of the fact that the initial problem is second-order in time, such that the velocity mapping can be expressed as a function of the displacement mapping, see e.g. (Vizzaccaro et al. 2022; Opreni et al. 2023; Vizzaccaro et al. 2024) for general discussions related to this point.

It is also possible to obtain the nonlinear mappings relating normal and modal coordinates, which are here given for the sake of completeness:

$$ 
\begin{aligned}
y_1 = z_1 + \left( \frac{h}{4\omega^2}z_1^3 - \frac{3h}{4\omega^2}z_1z_2^2 - \frac{h}{8\omega^2}z_2^3 \right) + \left( \frac{3h^2}{64\omega^4}z_1^5 - \frac{15h^2}{16\omega^4}z_1^4z_2 + \frac{69h^2}{32\omega^4}z_1^2z_2^3 + \frac{21h^2}{64\omega^4}z_1z_2^4 \right. \\
\left. -\frac{h^2}{32\omega^4}z_2^5 \right) + \left( \frac{h^3}{128\omega^6}z_1^7 - \frac{189h^3}{512\omega^6}z_1^6z_2 + \frac{1065h^3}{256\omega^6}z_1^5z_2^2 - \frac{2139h^3}{256\omega^6}z_1^4z_2^3 - \frac{561h^3}{512\omega^6}z_1^3z_2^4 - \frac{561h^3}{512\omega^6}z_1^2z_2^5 \right. \\
\left. +\frac{29h^3}{128\omega^6}z_1z_2^6 - \frac{3h^3}{512\omega^6}z_2^7 \right) + \left( \frac{5h^4}{4096\omega^8}z_1^9 - \frac{193h^4}{2048\omega^8}z_1^8z_2 + \frac{309h^4}{128\omega^8}z_1^7z_2^2 - \frac{20517h^4}{1024\omega^8}z_1^6z_2^3 \right. \\
\left. +\frac{75891h^4}{2048\omega^8}z_1^4z_2^5 + \frac{16689h^4}{4096\omega^8}z_1^3z_2^6 - \frac{1423h^4}{1024\omega^8}z_1^2z_2^7 + \frac{279h^4}{4096\omega^8}z_1z_2^8 - \frac{h^4}{1024\omega^8}z_2^9 \right) \\
+ \left( \frac{3h^5}{16384\omega^{10}}z_1^{11} - \frac{651h^5}{32768\omega^{10}}z_1^{10}z_2 + \frac{6677h^5}{8192\omega^{10}}z_1^9z_2^2 - \frac{489807h^5}{32768\omega^{10}}z_1^8z_2^3 + \frac{1661919h^5}{16384\omega^{10}}z_1^7z_2^4 \right. \\
\left. -\frac{727419h^5}{4096\omega^{10}}z_1^5z_2^6 - \frac{524913h^5}{32768\omega^{10}}z_1^4z_2^7 + \frac{133719h^5}{16384\omega^{10}}z_1^3z_2^8 - \frac{18729h^5}{32768\omega^{10}}z_1^2z_2^9 + \frac{255h^5}{16384\omega^{10}}z_1z_2^{10} \right. \\
\left. -\frac{5h^5}{32768\omega^{10}}z_2^{11} \right)
\end{aligned}
\quad (B.2a) $$

$$ y_2 = z_2 - \left( \frac{h}{8\omega^2}z_1^3 - \frac{3h}{4\omega^2}z_1^2z_2 + \frac{h}{4\omega^2}z_2^3 \right) - \left( \frac{h^2}{32\omega^4}z_1^5 + \frac{21h^2}{64\omega^4}z_1^4z_2 + \frac{69h^2}{32\omega^4}z_1^3z_2^2 - \frac{15h^2}{16\omega^4}z_1z_2^4 \right) $$
$$ 
\begin{aligned}
& +\frac{3h^2}{64\omega^4}z_2^5 - \left(\frac{3h^3}{512\omega^6}z_1^7 + \frac{29h^3}{128\omega^6}z_1^6z_2 - \frac{561h^3}{512\omega^6}z_1^5z_2^2 - \frac{2139h^3}{256\omega^6}z_1^4z_2^3 + \frac{1065h^3}{256\omega^6}z_1^2z_2^5\right. \\
& \left. -\frac{189h^3}{512\omega^6}z_1z_2^6 + \frac{h^3}{128\omega^6}z_2^7\right) - \left(\frac{h^4}{1024\omega^8}z_1^9 + \frac{279h^4}{4096\omega^8}z_1^8z_2 - \frac{1423h^4}{1024\omega^8}z_1^7z_2^2 + \frac{16689h^4}{4096\omega^8}z_1^6z_2^3\right. \\
& \left. +\frac{75891h^4}{2048\omega^8}z_1^5z_2^4 - \frac{20517h^4}{1024\omega^8}z_1^3z_2^6 + \frac{309h^4}{128\omega^8}z_1^2z_2^7 - \frac{193h^4}{2048\omega^8}z_1z_2^8 + \frac{5h^4}{4096\omega^8}z_2^9\right) \\
& \left(\frac{5h^5}{32768\omega^{10}}z_1^{11} + \frac{255h^5}{16384\omega^{10}}z_1^10z_2 - \frac{18729h^5}{32768\omega^{10}}z_1^9z_2^2 + \frac{133719h^5}{16384\omega^{10}}z_1^8z_2^3 - \frac{524913h^5}{32768\omega^{10}}z_1^7z_2^4\right. \\
& \left. -\frac{727419h^5}{4096\omega^{10}}z_1^6z_2^5 + \frac{1661919h^5}{16384\omega^{10}}z_1^4z_2^7 - \frac{489807h^5}{32768\omega^{10}}z_1^3z_2^8 + \frac{6677h^5}{8192\omega^{10}}z_1^2z_2^9 - \frac{651h^5}{32768\omega^{10}}z_1z_2^{10}\right. \\
& \left. +\frac{3h^5}{16384\omega^{10}}z_2^{11}\right)
\end{aligned} 
\quad (B.2b) $$

Note that the symbolic code MORFE_Symbolic allows producing these equations up to an arbitrary order, shown here up to order 11.

## C Additional results for the real normal form of the Duffing oscillator

In this appendix, more details on the real normal form (RNF) are given. First, for the sake of completeness and in order to draw out comparisons with CNF, the nonlinear mapping between the original and normal coordinates is given up to order 5:

$$ u = (z_1 + z_2) + \left( \frac{h}{8\omega^2} z_1^3 + \frac{h}{8\omega^2} z_2^3 \right) + \left( \frac{h^2}{64\omega^4} z_1^5 - \frac{21h^2}{64\omega^4} z_1^4 z_2 - \frac{21h^2}{64\omega^4} z_1 z_2^4 + \frac{h^2}{64\omega^4} z_2^5 \right) \quad (C.1a) $$

$$ v = (i\omega z_1 - i\omega z_2) + \left( i\frac{3h}{8\omega}z_1^3 - i\frac{3h}{8\omega}z_2^3 \right) + \left( i\frac{5h^2}{64\omega^3}z_1^5 - i\frac{27h^2}{64\omega^3}z_1^4z_2 + i\frac{27h^2}{64\omega^3}z_1z_2^4 - i\frac{5h^2}{64\omega^3}z_2^5 \right) \quad (C.1b) $$

These equations can be directly compared to those of the CNF in Equation (B.1). Since more monomials have been considered resonant and kept in the normal form for the RNF as compared to the CNF, the immediate consequence is that the nonlinear mapping for the CNF contains fewer terms.

Next, we present some calculations and approximations in order to derive an analytical backbone curve. To that purpose, the idea is to use assumptions and calculation procedures used in (Neild and Wagg 2011; Neild et al. 2015), that have been reinterpreted and reworked from the formalism derived in the present paper. The two main ideas used in (Neild and Wagg 2011; Neild et al. 2015) consist of using a mixed formulation with the initial coordinates (u, v) and the normal coordinates (z₁, z₂), and then to introduce a small bookkeeping parameter ε and asymptotic expansions to solve for the backbone curve order by order. Another key point is also to re-introduce oscillator-like equations with second-order derivatives in time. To that purpose, Equation (C.1b) up to order 5 can be differentiated with respect to time, yielding

$$ 
\begin{aligned}
\dot{v} = \ddot{u} = i\omega(z_1 - \dot{z}_2) + i\frac{9h}{8\omega}(z_1^2\dot{z}_1 - z_2^2\dot{z}_2) &+ i\frac{h^2}{64\omega^3}(25z_1^4\dot{z}_1 - 108z_1^3z_2\dot{z}_1 - 27z_1^4\dot{z}_2 + 27z_2^4\dot{z}_2 + 108z_1z_2^3\dot{z}_2 - 25z_2^4\dot{z}_2). 
\end{aligned} 
\quad (C.2) $$

Then, considering that the reduced dynamics up to order 5 is given by

$$ \dot{z}_1 = i\omega z_1 + i\frac{3h}{2\omega}(z_1^2 z_2 + z_1 z_2^2) - i\frac{3h^2}{16\omega^3}(5z_1^3 z_2^2 + 2z_1^2 z_2^3) \quad (C.3a) $$

$$ \dot{z}_2 = -i\omega z_2 - i\frac{3h}{2\omega}(z_1^2 z_2 + z_1 z_2^2) + i\frac{3h^2}{16\omega^3}(2z_1^3 z_2^2 + 5z_1^2 z_2^3), \quad (C.3b) $$

and introducing it on the previous equation while neglecting terms of order higher than 5, the following relationship can be found:
$$ \ddot{u} = -\omega^2 \left[ z_1 + z_2 + \frac{h}{8\omega^2} (z_1^3 + z_2^3) + \frac{h^3}{64\omega^4} (z_1^5 - 21z_1^4z_2 - 21z_1z_2^4 + z_2^5) \right] \\ - h \left[ z_1^3 + 3z_1^2z_2 + 3z_1z_2^3 + z_2^3 + \frac{3h}{8\omega^2} (z_1^5 + 2z_1^4z_2 + z_1z_2^4 + z_1^2z_2^3 + 2z_1z_2^4 + z_2^5) \right]. \quad (C.4) $$

In this last equation, the first term in bracket is exactly equal to $u$ while the second to $u^3$, meaning that this equation is a rewriting of the original problem (Duffing equation) where the nonlinear mapping between the initial and normal coordinates has been made explicit. In order to derive an analytical backbone curve, polar coordinates can be introduced as $z_1 = \rho/2e^{i\omega_{NL}t}$ and $z_2 = \rho/2e^{-i\omega_{NL}t}$, with $\rho$ and $\omega_{NL}$ constant in time, and where the phase of the harmonics is taken to be zero. Additionally, in the spirit of perturbation methods, a bookkeeping parameter $\epsilon$ is introduced, in order to scale different orders of magnitude. Specifically, it is assumed that $\rho$ is of $O(\epsilon)$, while $\omega_{NL}^2$ can be decomposed in the spirit of a perturbative solution as

$$ \omega_{NL}^2 = \omega_{1,0}^2 + \epsilon\omega_{1,1}^2 + \epsilon^2\omega_{1,2}^2 + \epsilon^3\omega_{1,3}^2 + \epsilon^4\omega_{1,4}^2 + O(\epsilon^5). \quad (C.5) $$

Inserting these assumptions into Equation (C.1a) and taking two time derivatives yields

$$ u = \epsilon \frac{\rho}{2} (e^{i\omega_{NL}t} + e^{-i\omega_{NL}t}) + \epsilon^3 \frac{h\rho^3}{64\omega^2} (e^{3i\omega_{NL}t} + e^{-3i\omega_{NL}t}) \\ + \epsilon^5 \frac{h^2\rho^5}{2048\omega^4} (e^{5i\omega_{NL}t} - 21e^{3i\omega_{NL}t} - 21e^{-3i\omega_{NL}t} + e^{-5i\omega_{NL}t}), \quad (C.6) $$

and

$$ \ddot{u} = -\epsilon \frac{\rho}{2} \omega_{1,0}^2 (e^{i\omega_{NL}t} + e^{-i\omega_{NL}t}) - \epsilon^2 \frac{\rho}{2} \omega_{1,1}^2 (e^{i\omega_{NL}t} + e^{-i\omega_{NL}t}) \\ - \epsilon^3 \left[ \frac{\rho}{2} \omega_{1,2}^2 (e^{i\omega_{NL}t} + e^{-i\omega_{NL}t}) + \frac{9h\rho^3}{64\omega^2} \omega_{1,0}^2 (e^{3i\omega_{NL}t} + e^{-3i\omega_{NL}t}) \right] \\ - \epsilon^4 \left[ \frac{\rho}{2} \omega_{1,3}^2 (e^{i\omega_{NL}t} + e^{-i\omega_{NL}t}) + \frac{9h\rho^3}{64\omega^2} \omega_{1,1}^2 (e^{3i\omega_{NL}t} + e^{-3i\omega_{NL}t}) \right] \\ - \epsilon^5 \left[ \frac{\rho}{2} \omega_{1,4}^2 (e^{i\omega_{NL}t} + e^{-i\omega_{NL}t}) + \left( \frac{9h^2\rho^3}{64\omega^2} \omega_{1,2}^2 + \frac{189h^2\rho^5}{2048\omega^4} \omega_{1,0}^2 \right) \right. \\ \left. \times (e^{3i\omega_{NL}t} + e^{-3i\omega_{NL}t}) + \frac{25h^2\rho^5}{2048\omega^4} \omega_{1,0}^2 (e^{5i\omega_{NL}t} + e^{-5i\omega_{NL}t}) \right]. \quad (C.7) $$

These expressions can be inserted, together with the polar representation of the normal variables, into Equation (C.4). Equating the different orders of $\epsilon$ gives

$$ \epsilon : \quad \omega_{1,0}^2 = \omega^2 \quad (C.8) $$

$$ \epsilon^2 : \quad \omega_{1,1}^2 = 0 \quad (C.9) $$

$$ \epsilon^3 : \quad \omega_{1,2}^2 = \frac{3h}{4} \rho^2 \quad (C.10) $$

$$ \epsilon^4 : \quad \omega_{1,3}^2 = 0 \quad (C.11) $$

$$ \epsilon^5 : \quad \omega_{1,4}^2 = \frac{3h^2}{128\omega^2} \rho^4 \quad (C.12) $$

A direct solution to these equations yields the following backbone curve:

$$ \omega_{NL}^2 = \omega^2 + \frac{3h}{4} \rho^2 + \frac{3h^2}{128\omega^2} \rho^4. \quad (C.13) $$

This expression coincides with the one given in (Neild et al. 2015). Interestingly, this solving procedure yields a backbone curve that is directly expressed with the square of the radian eigenfrequencies, a distinctive feature from other perturbative techniques, which appears to stem directly from the treatment of the second-order time derivative and the choice of the RNF. Also, in comparison to the backbone curve given by the CNF shown in Equation (11), where each new
order appears with alternate signs, here only positive coefficients are present, meaning that another approximation seems to be given. Note, however, that the procedure to compute the backbone with the RNF is tedious and difficult to automatise, as compared to the analytical solution provided with the CNF. Hence a complete comparison of both methods at arbitrary order is not pushed further here.

## D Nonlinear mapping, complex normal form for the Duffing oscillator with quadratic and cubic terms

For the sake of completeness, the nonlinear mapping for the Duffing equation with quadratic and cubic nonlinearity, Equation (25), is here reported, up to order 5. Higher order terms can be derived automatically and symbolically with MORFE_Symbolic.

$$ u = (z_1 + z_2) + \left( \frac{g}{3\omega^2} z_1^2 - \frac{2g}{\omega^2} z_1 z_2 + \frac{g}{3\omega^2} z_2^2 \right) + \left( \frac{2g^2 + 3h\omega^2}{24\omega^4} z_1^3 + \frac{10g^2 - 9h\omega^2}{12\omega^4} z_1^2 z_2 \right. \\ \left. + \frac{10g^2 - 9h\omega^2}{12\omega^4} z_1 z_2^2 + \frac{2g^2 + 3h\omega^2}{24\omega^4} z_2^3 \right) + \left( \frac{2g^3 + 9gh\omega^2}{108\omega^6} z_1^4 + \frac{178g^3 - 333gh\omega^2}{108\omega^6} z_1^3 z_2 \right. \\ \left. + \frac{-68g^3 + 117gh\omega^2}{9\omega^6} z_1^2 z_2^2 + \frac{178g^3 - 333gh\omega^2}{108\omega^6} z_1 z_2^3 + \frac{2g^3 + 9gh\omega^2}{108\omega^6} z_2^4 \right) \\ + \left( \frac{20g^4 + 180g^2h\omega^2 + 81h^2\omega^4}{5184\omega^8} z_1^5 + \frac{436g^4 - 444g^2h\omega^2 - 351h^2\omega^4}{576\omega^8} z_1^4 z_2 \right. \\ \left. + \frac{3740g^4 - 9468g^2h\omega^2 + 1863h^2\omega^4}{864\omega^8} z_1^3 z_2^2 + \frac{3740g^4 - 9468g^2h\omega^2 + 1863h^2\omega^4}{864\omega^8} z_1^2 z_2^3 \right. \\ \left. + \frac{436g^4 - 444g^2h\omega^2 - 351h^2\omega^4}{576\omega^8} z_1 z_2^4 + \frac{20g^4 + 180g^2h\omega^2 + 81h^2\omega^4}{5184\omega^8} z_2^5 \right) \quad (D.1a) $$

$$ v = (i\omega z_1 - i\omega z_2) + \left( i \frac{2g}{3\omega} z_1^2 - i \frac{2g}{3\omega} z_2^2 \right) + \left( \frac{2ig^2 + 3i\omega^2}{8\omega^3} z_1^3 + \frac{-10ig^2 + 9i\omega^2}{12\omega^3} z_1^2 z_2 \right. \\ \left. + \frac{10ig^2 - 9i\omega^2}{12\omega^3} z_1 z_2^2 + \frac{-2ig^2 - 3i\omega^2}{8\omega^3} z_2^3 \right) + \left( \frac{2ig^3 + 9igh\omega^2}{27\omega^5} z_1^4 + \frac{118ig^3 - 279igh\omega^2}{54\omega^5} z_1^3 z_2 \right. \\ \left. + \frac{-118ig^3 + 279igh\omega^2}{54\omega^5} z_1 z_2^3 + \frac{-2ig^3 - 9igh\omega^2}{27\omega^5} z_2^4 \right) + \left( \frac{100ig^4 + 900ig^2h\omega^2 + 405ih^2\omega^4}{5184\omega^7} z_1^5 \right. \\ \left. + \frac{356ig^4 - 492ig^2h\omega^2 - 243ih^2\omega^4}{192\omega^7} z_1^4 z_2 + \frac{-3740ig^4 + 9468ig^2h\omega^2 - 1863ih^2\omega^4}{864\omega^7} z_1^3 z_2^2 \right. \\ \left. + \frac{3740ig^4 - 9468ig^2h\omega^2 + 1863ih^2\omega^4}{864\omega^7} z_1^2 z_2^3 + \frac{-356ig^4 + 492ig^2h\omega^2 + 243ih^2\omega^4}{192\omega^7} z_1 z_2^4 \right. \\ \left. + \frac{-100ig^4 - 900ig^2h\omega^2 - 405ih^2\omega^4}{5184\omega^7} z_2^5 \right) \quad (D.1b) $$

Comparing these equations to Equation (B.1), one can see how the quadratic term complexifies the result in terms of the nonlinear mapping, which now contains all even and odd powers. Also, letting $g = 0$ in Equation (D.1), one directly retrieves Equation (B.1) as expected.

## E RNF and ONF developments for the Duffing oscillator with quadratic and cubic terms

This appendix completes the analysis as of the Duffing oscillator with quadratic and cubic terms, Equation (25), by providing the results for both real normal form (RNF) and oscillator normal form (ONF). The case of the RNF is considered first, and the reduced dynamics up to order 5 reads:

$$ \dot{z}_1 = i\omega z_1 + i \frac{-10g^2 + 9h\omega^2}{6\omega^3} z_1^2 z_2 + i \frac{-10g^2 + 9h\omega^2}{6\omega^3} z_1 z_2^2 + i \frac{-1940g^4 + 6228g^2h\omega^2 - 405h^2\omega^4}{432\omega^7} z_1^3 z_2^2 $$
$$ + i \frac{-1060g^4 + 3060g^2h\omega^2 - 81h^2\omega^4}{216\omega^7} z_1^2 z_2^3, \qquad (E.1) $$

with the second one being its complex conjugate. Additionally, the nonlinear mapping for the displacement is

$$ \begin{aligned} u = (z_1 + z_2) &+ \left( \frac{g}{3\omega^2} z_1^2 - \frac{2g}{\omega^2} z_1 z_2 + \frac{g}{3\omega^2} z_2^2 \right) + \left( \frac{2g^2 + 3h\omega^2}{24\omega^4} z_1^3 + \frac{2g^2 + 3h\omega^2}{24\omega^4} z_2^3 \right) \\ &+ \left( \frac{2g^3 + 9gh\omega^2}{108\omega^6} z_1^4 \frac{358g^3 - 495gh\omega^2}{108\omega^6} z_1^3 z_2 + \frac{-26g^3 + 42gh\omega^2}{3\omega^6} z_1^2 z_2^2 + \frac{358g^3 - 495gh\omega^2}{108\omega^6} z_1 z_2^3 \right) \\ &+ \left( \frac{2g^3 + 9gh\omega^2}{108\omega^6} z_2^4 \right) + \left( \frac{20g^4 + 180g^2h\omega^2 + 81h^2\omega^4}{5184\omega^8} z_1^5 + \frac{212g^4 - 268g^2h\omega^2 - 63h^2\omega^4}{192\omega^8} z_1^4 z_2 \right) \\ &+ \left( \frac{212g^4 - 268g^2h\omega^2 - 63h^2\omega^4}{192\omega^8} z_1 z_2^4 + \frac{20g^4 + 180g^2h\omega^2 + 81h^2\omega^4}{5184\omega^8} z_2^5 \right), \end{aligned} \qquad (E.2) $$

while the one for the velocity is simply its time derivative.

As compared to the results given by the CNF reported in Section 3.1, one can observe that the cubic coefficient for both CNF and RNF has the same expression, such that the prediction of the transition hardening/softening behaviour shall happen for $g^2 = 9h\omega^2/10$. However, for the next order, note that in Equation (E.1) the two quintic coefficients in front of $z_1^3 z_2^2$ and $z_1^2 z_2^3$ are not the same. Combined with the fact that deriving an explicit analytical backbone curve for the RNF for arbitrary order requires further assumptions, see the analysis reported in Appendix C, we conclude that it is more difficult to analyse the successive sign change of the different coefficients to infer the trend of the backbone curve in the RNF as compared to CNF.

For the ONF, the result reported here is limited to the third-order because its automation appears cumbersome and has therefore not been implemented in MORFE_Symbolic. Indeed, attempts to compute higher orders even in simple cases have shown that the processing of the term is more complex, see e.g. (Shami et al. 2022) for an example of processing the cubic terms with a second-order internal resonance. On the other hand, calculations up to order three have global analytical expressions reported in (Touzé et al. 2004; Touzé and Amabili 2006). Using these general formulae, we obtain the nonlinear change of coordinates

$$ u = U - g/(3\omega^2)U^2 - 2g/(3\omega^4)V^2, \qquad (E.3) $$

for Equation (25), while the second equation of the nonlinear mapping for $v$ is not reported for the sake of brevity, since it can be easily recovered using $v = \dot{u}$. The dynamics in the normal coordinates $(U, V)$ up to order three is

$$ \ddot{U} + \omega^2 U + (h - 2g^2/(3\omega^2)) U^3 - 4g^2/(3\omega^4) UV^2 = 0. \qquad (E.4) $$

Since oscillator equations are enforced by the method, two monomials are present in Equation (E.4), $U^3$ and $UV^2$. A backbone curve for Equation (E.4) can be derived with a first-order perturbative method, yielding the nonlinear oscillation frequency $\omega_{NL}$ as a function of the amplitude $a$ as (Touzé et al. 2004):

$$ \omega_{NL} = \omega(1 + \Gamma a^2), \qquad (E.5) $$

where $\Gamma$ combines the two coefficients of the cubic monomials and reads:

$$ \Gamma = (-10g^2 + 9h\omega^2)/(24\omega^4). \qquad (E.6) $$

The prediction of hardening/softening behaviour transition point for the first-order term of the backbone curve is thus also predicted with ONF to occur at $g^2 = 9h\omega^2/10$.

## F Nonlinear mapping, complex normal form for the forced-damped cubic Duffing oscillator: primary resonance

This appendix gives the nonlinear mappings for the forced cubic Duffing oscillator in a situation of primary resonance. The parametrisation is pushed to order 3 only, but higher order
parametrisations can be easily obtained with MORFE_Symbolic.

$$
\begin{aligned}
u &= \left( z_1 + z_2 + \frac{i\kappa}{4\xi\delta\omega^2 + 8i\delta^2\omega^2}z_3 - \frac{i\kappa}{4\xi\delta\omega^2 - 8i\delta^2\omega^2}z_4 \right) + \left( -\frac{h}{4(\xi^2\omega^2 - 2\delta^2\omega^2 - 3i\xi\delta\omega^2)}z_1^3 \right. \\
  &\quad + \frac{3ih}{4\xi\delta\omega^2 - 4i\delta^2\omega^2}z_1^2z_2 + \frac{3h\kappa(10\delta\omega + i\xi\omega)}{16\delta\omega(-\xi^4\omega^4 + 16\delta^4\omega^4 + 20i\xi\delta^3\omega^4 + 5i\xi^3\delta\omega^4)}z_1^2z_3 \\
  &\quad + \frac{3h\kappa\xi}{16\delta^2\omega(\xi\omega - i\delta\omega)(\xi\omega - 2i\delta\omega)^2}z_1^2z_4 - \frac{3ih}{4\xi\delta\omega^2 + 4i\delta^2\omega^2}z_1z_2^2 \\
  &\quad - \frac{3h\kappa(\xi\omega - 4i\delta\omega)}{8\delta^2\omega^2(\xi\omega - 2i\delta\omega)(\xi^2\omega^2 + 2\delta^2\omega^2 + i\xi\delta\omega^2)}z_1z_2z_3 \\
  &\quad - \frac{3h\kappa(\xi\omega + 4i\delta\omega)}{8\delta^2\omega^2(\xi\omega + i\delta\omega)(\xi\omega - 2i\delta\omega)(\xi\omega + 2i\delta\omega)}z_1z_2z_4 \\
  &\quad + \frac{9h\kappa^2(\xi\omega - 10i\delta\omega)}{64\delta^2\omega^2(\xi\omega - i\delta\omega)(\xi\omega + 2i\delta\omega)^2(\xi^2\omega^2 - 8\delta^2\omega^2 - 6i\xi\delta\omega^2)}z_1z_3^2 \\
  &\quad + \frac{3ih\kappa^2}{32\delta^3\omega^3(\xi\omega - i\delta\omega)(\xi\omega - 2i\delta\omega)^2}z_1z_3z_4 \\
  &\quad + \frac{3ih\kappa^2(\xi\omega + 10i\delta\omega)}{64\delta^3\omega^3(\xi\omega + i\delta\omega)(\xi\omega - 2i\delta\omega)^2(\xi\omega + 2i\delta\omega)}z_1z_4^2 - \frac{h}{4(\xi^2\omega^2 - 2\delta^2\omega^2 + 3i\xi\delta\omega^2)}z_2^3 \\
  &\quad + \frac{3h\kappa\xi}{16\delta^2\omega(\xi\omega + i\delta\omega)(\xi\omega + 2i\delta\omega)^2}z_2^2z_3 \\
  &\quad + \frac{3ih\kappa(\xi\omega + 10i\delta\omega)}{16\delta\omega(\xi^4\omega^4 - 16\delta^4\omega^4 + 20i\xi\delta^3\omega^4 + 5i\xi^3\delta\omega^4)}z_2^2z_4 \\
  &\quad - \frac{3ih\kappa^2(\xi\omega - 10i\delta\omega)}{64\delta^3\omega^3(\xi\omega - i\delta\omega)(\xi\omega - 2i\delta\omega)(\xi\omega + 2i\delta\omega)^2}z_2z_3^2 \\
  &\quad - \frac{3ih\kappa^2}{32\delta^3\omega^3(\xi\omega + i\delta\omega)(\xi\omega + 2i\delta\omega)^2}z_2z_3z_4 \\
  &\quad + \frac{9h\kappa^2(\xi\omega + 10i\delta\omega)}{64\delta^2\omega^2(\xi\omega + i\delta\omega)(\xi\omega - 2i\delta\omega)^2(\xi^2\omega^2 - 8\delta^2\omega^2 + 6i\xi\delta\omega^2)}z_2z_4^2 \\
  &\quad + \frac{h\kappa^3(7i\xi^3\omega^3 - 872\delta^3\omega^3 + 340i\xi\delta^2\omega^3 + 10\xi^2\delta\omega^3)}{256\delta^3\omega^3(\xi\omega - i\delta\omega)(\xi\omega + 2i\delta\omega)^4(\xi\omega + 4i\delta\omega)(\xi^2\omega^2 - 8\delta^2\omega^2 - 6i\xi\delta\omega^2)}z_3^3 \\
  &\quad - \frac{3h\kappa^3(\xi^2\omega^2 - 24\delta^2\omega^2 + 2i\xi\delta\omega^2)}{256\delta^4\omega^4(\xi\omega - i\delta\omega)(\xi\omega - 2i\delta\omega)^2(\xi\omega + 2i\delta\omega)^3}z_3^2z_4 \\
  &\quad - \frac{3h\kappa^3(\xi^2\omega^2 - 24\delta^2\omega^2 - 2i\xi\delta\omega^2)}{256\delta^4\omega^4(\xi\omega + i\delta\omega)(\xi\omega - 2i\delta\omega)^3(\xi\omega + 2i\delta\omega)^2}z_3z_4^2 \\
  &\quad + \frac{h\kappa^3(-7i\xi^3\omega^3 - 872\delta^3\omega^3 - 340i\xi\delta^2\omega^3 + 10\xi^2\delta\omega^3)}{256\delta^3\omega^3(\xi\omega + i\delta\omega)(\xi\omega - 2i\delta\omega)^4(\xi\omega - 4i\delta\omega)(\xi^2\omega^2 - 8\delta^2\omega^2 + 6i\xi\delta\omega^2)}z_4^3
\end{aligned}
\quad (F.1)
$$

$$
\begin{aligned}
v &= \left( (-\xi\omega + i\delta\omega)z_1 + (-\xi\omega - i\delta\omega)z_2 + \frac{\kappa(\delta\omega - i\xi\omega)}{4\delta\omega(\xi\omega + 2i\delta\omega)}z_3 + \frac{\kappa\delta\omega + i\kappa\xi\omega}{4\xi\delta\omega^2 - 8i\delta^2\omega^2}z_4 \right) \\
  &\quad + \left( \frac{3h}{4\xi\omega - 8i\delta\omega}z_1^3 + \frac{3h(\xi\omega + i\delta\omega)}{4\delta\omega(\delta\omega + i\xi\omega)}z_1^2z_2 \right. \\
  &\quad + \frac{3h\kappa(22\delta\omega + 3i\xi\omega)}{16\delta\omega(\xi\omega + 2i\delta\omega)(\xi^2\omega^2 - 8\delta^2\omega^2 - 6i\xi\delta\omega^2)}z_1^2z_3 \\
  &\quad - \frac{3h\kappa\xi(\xi\omega + i\delta\omega)}{16\delta^2\omega(\xi\omega - i\delta\omega)(\xi\omega - 2i\delta\omega)^2}z_1^2z_4 + \frac{3h(\delta\omega + i\xi\omega)}{4\delta\omega(\xi\omega + i\delta\omega)}z_1z_2^2 \\
  &\quad + \frac{3h\kappa(\xi\omega + i\delta\omega)(\xi\omega - 4i\delta\omega)}{8\delta^2\omega^2(\xi\omega - 2i\delta\omega)(\xi^2\omega^2 + 2\delta^2\omega^2 + i\xi\delta\omega^2)}z_1z_2z_3 \\
  &\quad + \frac{3h\kappa(\xi\omega - i\delta\omega)(\xi\omega + 4i\delta\omega)}{8\delta^2\omega^2(\xi\omega + i\delta\omega)(\xi\omega - 2i\delta\omega)(\xi\omega + 2i\delta\omega)}z_1z_2z_4
\end{aligned}
$$
$$
\begin{aligned}
& -\frac{15h\kappa^2 (\xi\omega - 10i\delta\omega)}{64\delta^2\omega^2 (\xi\omega + 2i\delta\omega)^2 (\xi^2\omega^2 - 8\delta^2\omega^2 - 6i\xi\delta\omega^2)} z_1 z_3^2 \\
& +\frac{3h\kappa^2 (\delta\omega - i\xi\omega)}{32\delta^4\omega^3 (\xi\omega - i\delta\omega) (\xi\omega - 2i\delta\omega)^2} z_1 z_3 z_4 \\
& +\frac{3h\kappa^2 (\xi\omega - i\delta\omega) (10\delta\omega - i\xi\omega)}{64\delta^4\omega^3 (\xi\omega + i\delta\omega) (\xi\omega - 2i\delta\omega)^2 (\xi\omega + 2i\delta\omega)} z_1 z_4^2 + \frac{3h}{4\xi\omega + 8i\delta\omega} z_2^3 \\
& -\frac{3h\kappa\xi (\xi\omega - i\delta\omega)}{16\delta^2\omega (\xi\omega + i\delta\omega) (\xi\omega + 2i\delta\omega)^2} z_2^2 z_3 \\
& +\frac{3h\kappa (22\delta\omega - 3i\xi\omega)}{16\delta\omega (\xi\omega - 2i\delta\omega) (\xi^2\omega^2 - 8\delta^2\omega^2 + 6i\xi\delta\omega^2)} z_2^2 z_4 \\
& +\frac{3h\kappa^2 (\xi\omega + i\delta\omega) (10\delta\omega + i\xi\omega)}{64\delta^4\omega^3 (\xi\omega - i\delta\omega) (\xi\omega - 2i\delta\omega) (\xi\omega + 2i\delta\omega)^2} z_2 z_3^2 \\
& +\frac{3h\kappa^2 (\delta\omega + i\xi\omega)}{32\delta^4\omega^3 (\xi\omega + i\delta\omega) (\xi\omega + 2i\delta\omega)^2} z_2 z_3 z_4 \\
& -\frac{15h\kappa^2 (\xi\omega + 10i\delta\omega)}{64\delta^2\omega^2 (\xi\omega - 2i\delta\omega)^2 (\xi^2\omega^2 - 8\delta^2\omega^2 + 6i\xi\delta\omega^2)} z_2 z_4^2 \\
& +\frac{3h\kappa^3 (-3i\xi^3\omega^3 + 392\delta^4\omega^3 - 180i\xi\delta^2\omega^3 - 10\xi^2\delta\omega^3)}{256\delta^4\omega^3 (\xi\omega + 2i\delta\omega)^4 (\xi\omega + 4i\delta\omega) (\xi^2\omega^2 - 8\delta^2\omega^2 - 6i\xi\delta\omega^2)} z_3^3 \\
& +\frac{3h\kappa^3 (\xi\omega + i\delta\omega) (\xi^2\omega^2 - 24\delta^2\omega^2 + 2i\xi\delta\omega^2)}{256\delta^4\omega^4 (\xi\omega - i\delta\omega) (\xi\omega - 2i\delta\omega)^2 (\xi\omega + 2i\delta\omega)^3} z_3^2 z_4 \\
& +\frac{3h\kappa^3 (\xi\omega - i\delta\omega) (\xi^2\omega^2 - 24\delta^2\omega^2 - 2i\xi\delta\omega^2)}{256\delta^4\omega^4 (\xi\omega + i\delta\omega) (\xi\omega - 2i\delta\omega)^3 (\xi\omega + 2i\delta\omega)^2} z_3 z_4^2 \\
& +\frac{3h\kappa^3 (3i\xi^3\omega^3 + 392\delta^4\omega^3 + 180i\xi\delta^2\omega^3 - 10\xi^2\delta\omega^3)}{256\delta^4\omega^3 (\xi\omega - 2i\delta\omega)^4 (\xi\omega - 4i\delta\omega) (\xi^2\omega^2 - 8\delta^2\omega^2 + 6i\xi\delta\omega^2)} z_4^3
\end{aligned}
\quad (F.2)
$$

G Reduced dynamics coefficients, complex normal form for the forced-damped cubic Duffing oscillator at primary resonance computed for $\Omega_p = \omega$

This appendix gives coefficients of the complex normal form (CNF) of the forced-damped cubic Duffing oscillator in a primary resonance scenario. Differently from the main text where the parametrisation has been computed at the value $\Omega_p = \delta\omega$, here we show the analytical obtained result when the parametrisation is computed for $\Omega_p = \omega$. The eight coefficients of the resonant monomials shown in Equation (35) are then

$$
\begin{aligned}
f_1 &= \lambda_1, \quad f_2 = i \frac{3h}{2\delta\omega}, \quad f_3 = -i \frac{\kappa}{4\delta\omega}, \quad f_4 = \frac{3h\kappa (1 - \delta + i\xi)}{16\xi\delta^2\omega^3} \\
f_5 &= -\frac{3h\kappa (1 - \delta - i\xi)}{8\xi\delta^2\omega^3}, \quad f_6 = \frac{3ih\kappa^2 (1 - \delta)}{32\xi^2\delta^3\omega^5}, \quad (G.1) \\
f_7 &= -\frac{3ih\kappa^2 (-\xi^2 + i\delta\xi - \delta - i\xi + 1)}{64\xi^2\delta^3\omega^5}, \quad f_8 = \frac{3h\kappa^3 (\xi^2 - i\delta\xi + 2\delta + i\xi - 2)}{512\xi^3\delta^4\omega^7}.
\end{aligned}
$$

Note that considering a small damping hypothesis, these coefficients are equivalent to those given by Equation (36). However, in the general case, the expressions differ slightly and are longer in the present case.

H Nonlinear mapping, complex normal form for the forced cubic Duffing oscillator for the 3:1 superharmonic resonance

This appendix gives the nonlinear displacement mapping for the forced cubic Duffing oscillator when a 3:1 superharmonic forcing is considered. The parametrisation is pushed to order 3 only,
but higher order parametrisations can be easily obtained with MORFE_Symbolic. The velocity mapping can be obtained as the derivative of the displacement.

$$
\begin{aligned}
u = & \left( z_1 + z_2 - \frac{9i\kappa}{4(3\xi - 4i)\omega^2} z_3 + \frac{9i\kappa}{4(3\xi + 4i)\omega^2} z_4 \right) + \left( \frac{h}{4(3\xi(i\delta - \xi) + 2)\omega^2} z_1^3 \right. \\
& - \frac{3h}{4(\xi(i\delta - \xi) + 1)\omega^2} z_1^2 z_2 + \frac{243i\hbar\kappa}{8(3\xi - 4i)(3\xi(-6i\delta + 6\xi - i) - 6\delta - 14)\omega^4} z_1^2 z_3 \\
& - \frac{243i\hbar\kappa}{8(3\xi + 4i)(3\xi(-6i\delta + 6\xi + i) + 6\delta - 14)\omega^4} z_1^2 z_4 + \frac{3h(\xi - i\delta)}{4i\delta\omega^2} z_1 z_2^2 \\
& - \frac{243h\kappa}{4(9\xi^2 + 16)\omega^4} z_1 z_2 z_3 - \frac{243h\kappa}{4(9\xi^2 + 16)\omega^4} z_1 z_2 z_4 - \frac{2187i\hbar\kappa^2}{64(-3\xi + 4i)^2(3i\delta + i)\omega^6} z_1 z_3^2 \\
& + \frac{243h\kappa^2}{32(9\xi^4 + 7\xi^2 - 16)\omega^6} z_1 z_3 z_4 + \frac{2187i\hbar\kappa^2}{64(3\xi + 4i)^2(3i\delta - i)\omega^6} z_1 z_4^2 - \frac{h}{4(3\xi(i\delta + \xi) - 2)\omega^2} z_2^3 \\
& + \frac{243i\hbar\kappa}{8(3\xi - 4i)(3\xi(6i\delta + 6\xi - i) + 6\delta - 14)\omega^4} z_2^2 z_3 \\
& - \frac{243i\hbar\kappa}{8(3\xi + 4i)(3\xi(6i\delta + 6\xi + i) - 6\delta - 14)\omega^4} z_2^2 z_4 \\
& + \frac{2187i\hbar\kappa^2}{64(-3\xi + 4i)^2(3i\delta - i)\omega^6} z_2 z_3^2 + \frac{243h\kappa^2}{32(9\xi^4 + 7\xi^2 - 16)\omega^6} z_2 z_3 z_4 \\
& - \frac{2187i\hbar\kappa^2}{64(3\xi + 4i)^2(3i\delta + i)\omega^6} z_2 z_4^2 + \frac{729i\hbar\kappa^3}{128(3\xi - 4i)^3(\xi(i\delta + \xi) - \delta - 1)\omega^8} z_3^3 \\
& + \frac{19683h\kappa^3}{128(3\xi - 4i)^3(3\xi + 4i)\omega^8} z_3^2 z_4 + \frac{19683h\kappa^3}{128(3\xi - 4i)(3\xi + 4i)^3\omega^8} z_3 z_4^2 \\
& \left. - \frac{729h\kappa^3(i\delta + \xi - i)}{256\xi(3\xi + 4i)^3i\delta\omega^8} z_4^3 \right)
\end{aligned}
\qquad (H.1)
$$

# I Nonlinear mapping, complex normal form for the forced cubic Duffing oscillator for the 1:3 subharmonic resonance

This appendix gives the nonlinear displacement mapping for the forced cubic Duffing oscillator when a 1:3 subharmonic forcing is considered. The parametrisation is pushed to order 3 only, but higher order parametrisations can be easily obtained with MORFE_Symbolic. The velocity mapping can be obtained as the derivative of the displacement.

$$
\begin{aligned}
u = & \left( z_1 + z_2 - \frac{i\kappa}{4(3\xi + 4i)\omega^2} z_3 + \frac{i\kappa}{4(3\xi - 4i)\omega^2} z_4 \right) + \left( \frac{h}{4(3\xi(i\delta - \xi) + 2)\omega^2} z_1^3 \right. \\
& - \frac{3h}{4(\xi(i\delta - \xi) + 1)\omega^2} z_1^2 z_2 + \frac{3i\hbar\kappa}{8(3\xi + 4i)(\xi(-2i\delta + 2\xi - 3i) - 6\delta - 6)\omega^4} z_1^2 z_3 \\
& + \frac{3h\kappa(\xi(i\delta + \xi) - 3\delta - 1)}{16(9\xi^4 + 7\xi^2 - 16)\omega^4} z_1^2 z_4 + \frac{3h(\xi - i\delta)}{4i\delta\omega^2} z_1 z_2^2 \\
& - \frac{3h\kappa}{4(9\xi^2 + 16)\omega^4} z_1 z_2 z_3 - \frac{3h\kappa}{4(9\xi^2 + 16)\omega^4} z_1 z_2 z_4 - \frac{i\hbar\kappa^2}{64(3\xi + 4i)^2(i\delta + 3i)\omega^6} z_1 z_3^2 \\
& + \frac{3h\kappa^2}{32(9\xi^4 + 7\xi^2 - 16)\omega^6} z_1 z_3 z_4 + \frac{i\hbar\kappa^2}{64(-3\xi + 4i)^2(i\delta - 3i)\omega^6} z_1 z_4^2 \\
& - \frac{h}{4(3\xi(i\delta + \xi) - 2)\omega^2} z_2^3 + \frac{3i\hbar\kappa}{8(3\xi + 4i)(\xi(i\delta + \xi) + 3\delta - 1)\omega^4} z_2^2 z_3 \\
& - \frac{3i\hbar\kappa}{8(3\xi - 4i)(\xi(2i\delta + 2\xi + 3i) - 6\delta - 6)\omega^4} z_2^2 z_4 \\
& + \frac{i\hbar\kappa^2}{64(3\xi + 4i)^2(i\delta - 3i)\omega^6} z_2 z_3^2 + \frac{3h\kappa^2}{32(9\xi^4 + 7\xi^2 - 16)\omega^6} z_2 z_3 z_4 \\
& - \frac{i\hbar\kappa^2}{64(-3\xi + 4i)^2(i\delta + 3i)\omega^6} z_2 z_4^2 - \frac{h\kappa^3}{128(3\xi + 4i)^3(9\xi + 40i)\omega^8} z_3^3 \\
& \left. - \frac{h\kappa^3}{64(-3\xi + 4i)^3(i\delta + 3i)\omega^6} z_3 z_4^2 - \frac{h\kappa^3}{128(3\xi + 4i)^3(9\xi + 40i)\omega^8} z_3 z_4^3 \right)
\end{aligned}
$$
$$ + \frac{3h\kappa^3}{128(3\xi - 4i)(3\xi + 4i)^3\omega^8} z_3^2 z_4 + \frac{3h\kappa^3}{128(3\xi - 4i)^3(3\xi + 4i)\omega^8} z_3 z_4^2 - \frac{h\kappa^3}{128(3\xi - 4i)^3(9\xi - 40i)\omega^8} z_4^3 \quad (I.1) $$

## J Detailed calculation of the existence region for the 1:3 subharmonic resonance

This appendix is concerned with the derivation of the expression that gives the boundary of the region such that nontrivial solutions exist for the 1:3 subharmonic resonance of the cubic Duffing oscillator. The starting point is Equation (57). Upon substitution of coefficients $A_r, A_i, B_r$ and $B_i$, it is possible to find a biquadratic equation in $\rho$:

$$ \underbrace{\frac{9h^2}{64\delta^2\omega^2}}_a \rho^4 + \underbrace{\left[ \frac{3h}{4\delta\omega} \left( \delta\omega - \frac{\Omega}{3} \right) + \frac{27h^2\kappa^2}{256\delta^2\omega^6 (16 + 9\xi^2)} \right]}_b \rho^2 + \underbrace{\left( \delta\omega - \frac{\Omega}{3} \right)^2 + \frac{9h^2\kappa^4}{256\delta^2\omega^{10} (16 + 9\xi^2)^2} + \frac{3h\kappa^2}{8\delta\omega^5 (16 + 9\xi^2)} \left( \delta\omega - \frac{\Omega}{3} \right) + \xi^2\omega^2}_{c} = 0. \quad (J.1) $$

For the existence of real solutions, two conditions have to be satisfied, namely $b < 0$ and $\Delta = b^2 - 4ac \ge 0$. After algebraic manipulations, the first condition yields

$$ \Omega > 3\delta\omega + \frac{27h\kappa^2}{64\delta\omega^5 (16 + 9\xi^2)}, \quad (J.2) $$

while the second gives

$$ \frac{63h^2}{2048\delta^2\omega^{10}(16 + 9\xi^2)^2}\kappa^4 - \frac{h(\Omega - 3\delta\omega)}{16\delta\omega^5(16 + 9\xi^2)}\kappa^2 + 2\xi^2\omega^2 \le 0, \quad (J.3) $$

which is in itself a biquadratic inequality in $\kappa$. The above expression is negative if $\kappa^2$ is between its two roots, such that the boundaries of the existence region are given by

$$ \frac{63h\kappa^2}{64\delta\xi\omega^6(16 + 9\xi^2)} = \frac{\Omega - 3\delta\omega}{\xi\omega} \pm \sqrt{\frac{(\Omega - 3\delta\omega)^2}{\xi^2\omega^2} - 63}. \quad (J.4) $$

Explicitly isolating $\Omega$ gives the sought expression for existence condition of subharmonic resonance as:

$$ \Omega = 3\delta\omega + \frac{63h\kappa^2}{128\delta\omega^5(16 + 9\xi^2)} + \frac{32\delta\omega^7\xi^2(16 + 9\xi^2)}{h\kappa^2}. \quad (J.5) $$

Note that Equation (J.5) automatically guarantees that Equation (J.2) is verified.

## K Linear stability analysis for the parametrically excited oscillator

This Appendix is concerned with the linear stability analysis of the bifurcated branches for the case of parametric excitation, see Equation (66) in Section 3.3. The dynamical system from Equation (70) serves as the starting point. Assuming $\rho \ne 0$, substituting the coefficients from Equation (71), and considering the fixed points $\dot{\rho} = \dot{\psi} = 0$, it can be rewritten as

$$ \rho f_1^R + \frac{\rho^3}{4} f_2^R + \rho f_4^R + \rho f_3^R \cos \psi + \rho f_3^I \sin \psi = 0 \quad (K.1a) $$

$$ 2f_1^I - \Omega + \frac{\rho^2}{2} f_2^I + 2f_4^I - 2f_3^R \sin \psi + 2f_3^I \cos \psi = 0. \quad (K.1b) $$
Then, in order to perform the linear stability analysis, the Jacobian

$$ J = \begin{bmatrix} f_1^R + \frac{3\rho^2}{4} f_2^R + f_4^R + f_3^R \cos \psi + f_3^I \sin \psi & -\rho f_3^R \sin \psi + \rho f_3^I \cos \psi \\ \rho f_2^I & -2f_3^R \cos \psi - 2f_3^I \sin \psi \end{bmatrix} \quad (K.2a) $$

$$ = \begin{bmatrix} 0 & \rho \left( \frac{\Omega}{2} - \delta\omega - \rho^2 \frac{3h}{8\delta\omega} - \frac{\kappa^2(\delta-1)}{32\omega^3 \xi^2 \delta^2} \right) \\ \rho \frac{3h}{2\delta\omega} & -2\xi\omega \end{bmatrix}, \quad (K.2b) $$

of the vector field of the dynamical system from Equation (70) is calculated, where the coefficient expressions from Equation (69) have been substituted and Equation (K.1) has been used in order to simplify the Jacobian. Then, calculating its eigenvalues corresponds to solving the following two quadratic equations in $\lambda$:

$$ \lambda^2 + 2\xi\omega\lambda \pm \rho^2 \frac{3h}{2\delta\omega} \sqrt{\frac{\kappa^2}{16\omega^2\delta^2} - \xi^2\omega^2} = 0, \quad (K.3) $$

which has been obtained from det(J - λI) = 0 by substituting $\Omega$ from Equation (76). If the plus sign is considered first, one has

$$ \lambda = -\xi\omega \pm \sqrt{\xi^2\omega^2 + \rho^2 \frac{3h}{2\delta\omega} \sqrt{\frac{\kappa^2}{16\omega^2\delta^2} - \xi^2\omega^2}} \quad (K.4) $$

In this case, the real part of at least one eigenvalue is always positive, since the term under the square root is larger than $\xi\omega$. Thus, the right branch of the FRC, corresponding to the plus sign in Equation (76), is unstable. Analogously, for the minus sign case, the eigenvalues are

$$ \lambda = -\xi\omega \pm \sqrt{\xi^2\omega^2 - \rho^2 \frac{3h}{2\delta\omega} \sqrt{\frac{\kappa^2}{16\omega^2\delta^2} - \xi^2\omega^2}} \quad (K.5) $$

Therefore, using a similar argument, the left branch of the FRC is stable. From the stability of the bifurcated branches, the one from the main branch can be deduced due to the topology of the attractors in phase space: it is unstable between the bifurcation points and stable outside them. Note that the above analysis assumes a hardening Duffing oscillator with $h > 0$. In the softening case with $h < 0$, the FRC is bent to the left, and the analysis before has just to be changed: the unstable branch corresponds to the left bifurcation point, with smaller frequency, while the stable branch corresponds to the right point with higher frequency.

## L Nonlinear mappings for the two dofs system, free of internal resonance

The nonlinear mappings for the two dofs system with only cubic coefficients and no internal resonance are given below. Calculations were pursued up to order 5, and only displacement mappings are given since velocity ones can be obtained as their time derivatives.

$$ \begin{aligned} u_1 = (z_1 + z_2) &+ \left( \frac{h_{111}^1}{8\omega_1^2} z_1^3 - \frac{3h_{111}^1}{4\omega_1^2} z_1^2 z_2 - \frac{3h_{111}^1}{4\omega_1^2} z_1 z_2^2 + \frac{h_{111}^1}{8\omega_1^2} z_2^3 \right) \\ &+ \left( \frac{8h_{111}^1}{576\omega_1^6} \omega_1^2 + 9h_{111}^1}{64\omega_1^4\omega_2^2} \omega_1^2 - h_{111}^1}{390h_{111}^1}{69h_{111}^1} \omega_1^2 \omega_2^2 - 351h_{111}^1}{390h_{111}^1}{69h_{111}^1} \omega_2^4 \right) z_1^4 z_2 \\ &+ \frac{696h_{111}^1}{576\omega_1^8} \omega_1^4 - 120h_{111}^1}{64\omega_1^6\omega_2^2} \omega_1^2 \omega_2^2 - 351h_{111}^1}{64\omega_1^4\omega_2^4} \omega_1^2 + 390h_{111}^1}{69h_{111}^1}{69h_{111}^1} \omega_1^2 \omega_2^2 - 39h_{111}^1}{69h_{111}^1}{69h_{111}^1} \omega_2^4 \right) z_1^3 z_2 \\ &+ \frac{-1968h_{111}^1}{288\omega_1^8} \omega_1^4 + 240h_{111}^1}{320\omega_1^6\omega_2^2} \omega_1^2 \omega_2^2 + 621h_{111}^1}{320\omega_1^4\omega_2^4} \omega_1^2 - 690h_{111}^1}{69h_{111}^1}{69h_{111}^1} \omega_1^2 \omega_2^2 + 69h_{111}^1}{69h_{111}^1}{69h_{111}^1} \omega_2^4 \right) z_1^2 z_2^2 \\ &+ \frac{-1968h_{111}^1}{288\omega_1^8} \omega_1^4 + 240h_{111}^1}{320\omega_1^6\omega_2^2} \omega_1^2 \omega_2^2 + 621h_{111}^1}{320\omega_1^4\omega_2^4} \omega_1^2 - 690h_{111}^1}{69h_{111}^1}{69h_{111}^1} \omega_1^2 \omega_2^2 + 69h_{111}^1}{69h_{111}^1}{69h_{111}^1} \omega_2^4 \right) z_1^2 z_2^3 \end{aligned} $$
$$ 
\begin{aligned}
& + \frac{696h_{111}^2 \omega_1^4 - 120h_{111}^2 \omega_1^2 \omega_2^2 - 351h_{111}^1 \omega_1^4 + 390h_{111}^1 \omega_1^2 \omega_2^2 - 39h_{111}^1 \omega_2^4}{576\omega_1^8 - 640\omega_1^6\omega_2^2 + 64\omega_1^4\omega_2^4} z_1 z_2^4 \\
& + \frac{8h_{111}^2 \omega_1^2 + 9h_{111}^1 \omega_1^2 - h_{111}^1 \omega_2^2}{576\omega_1^6 - 64\omega_1^4\omega_2^2} z_2^5
\end{aligned}
\quad (L.1a) $$

$$ 
\begin{aligned}
u_2 = & \left( \frac{h_{111}^2}{9\omega_1^2 - \omega_2^2} z_1^3 + \frac{3h_{111}^2}{\omega_1^2 - \omega_2^2} z_1^2 z_2 + \frac{3h_{111}^2}{\omega_1^2 - \omega_2^2} z_1 z_2^2 + \frac{h_{111}^2}{9\omega_1^2 - \omega_2^2} z_2^3 \right) \\
& + \left( \frac{27h_{111}^2 h_{111}^1 \omega_1^2 - 3h_{111}^2 h_{111}^1 \omega_2^2 + 8h_{111}^2 h_{122}^1 \omega_1^2}{1800\omega_1^6 - 272\omega_1^4\omega_2^2 + 8\omega_1^2\omega_2^4} z_1^5 \right. \\
& + \frac{C_1}{162\omega_1^8 - 198\omega_1^6\omega_2^2 + 38\omega_1^4\omega_2^4 - 2\omega_1^2\omega_2^6} z_1^4 z_2 + \frac{C_2}{72\omega_1^8 - 152\omega_1^6\omega_2^2 + 88\omega_1^4\omega_2^4 - 8\omega_1^2\omega_2^6} z_1^3 z_2^2 \\
& + \frac{C_2}{72\omega_1^8 - 152\omega_1^6\omega_2^2 + 88\omega_1^4\omega_2^4 - 8\omega_1^2\omega_2^6} z_1^2 z_2^3 + \frac{C_1}{162\omega_1^8 - 198\omega_1^6\omega_2^2 + 38\omega_1^4\omega_2^4 - 2\omega_1^2\omega_2^6} z_1 z_2^4 \\
& \left. + \frac{27h_{111}^2 h_{111}^1 \omega_1^2 - 3h_{111}^2 h_{111}^1 \omega_2^2 + 8h_{111}^2 h_{122}^1 \omega_1^2}{1800\omega_1^6 - 272\omega_1^4\omega_2^2 + 8\omega_1^2\omega_2^4} z_2^5 \right)
\end{aligned}
\quad (L.1b) $$

With the coefficients on the numerators given by:

$$ C_1 = -81h_{111}^2 h_{111}^1 \omega_1^4 + 84h_{111}^2 h_{111}^1 \omega_1^2 \omega_2^2 - 3h_{111}^2 h_{111}^1 \omega_2^4 + 58h_{111}^2 h_{111}^1 h_{122}^1 \omega_1^4 \\ - 10h_{111}^2 h_{111}^1 h_{122}^1 \omega_1^2 \omega_2^2 \quad (L.2) $$

$$ C_2 = -1107h_{111}^2 h_{111}^1 \omega_1^4 + 582h_{111}^2 h_{111}^1 \omega_1^2 \omega_2^2 - 51h_{111}^2 h_{111}^1 \omega_2^4 + 656h_{111}^2 h_{111}^1 h_{122}^1 \omega_1^4 \\ - 80h_{111}^2 h_{111}^1 h_{122}^1 \omega_1^2 \omega_2^2 \quad (L.3) $$

## References

Arafat, H. N. and A. H. Nayfeh (2003). Non-linear responses of suspended cables to primary resonance excitation. *Journal of Sound and Vibration* 266:325–354. [DOI].

Benacchio, S., C. Giraud-Audine, and O. Thomas (2022). Effect of dry friction on a parametric nonlinear oscillator. *Nonlinear Dynamics* 108:1005–1026. [DOI], [HAL].

Breunung, T. and G. Haller (2018). Explicit backbone curves from spectral submanifolds of forced-damped nonlinear mechanical systems. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* 474(2213):20180083. [DOI], [OA], [ARXIV].

Cabré, X., E. Fontich, and R. de la Llave (2003a). The parameterization method for invariant manifolds. I. Manifolds associated to non-resonant subspaces. *Indiana University Mathematics Journal* 52(2):283–328. [URL], [OAI].

Cabré, X., E. Fontich, and R. de la Llave (2003b). The parameterization method for invariant manifolds. II. Regularity with respect to parameters. *Indiana University Mathematics Journal* 52(2):329–360. [URL], [OAI].

Cabré, X., E. Fontich, and R. de la Llave (2005). The parameterization method for invariant manifolds. III. Overview and applications. *Journal of Differential Equations* 218(2):444–515. [DOI], [OA], [OAI].

Dhooge, A., W. Govaerts, and Y. A. Kuznetsov (2004). MATCONT: a Matlab package for numerical bifurcation analysis of ODEs. *ACM SIGSAM Bulletin* 38(1):21–22. [DOI].

Elliott, A. J., A. Cammarano, S. A. Neild, T. L. Hill, and D. J. Wagg (2018). Comparing the direct normal form and multiple scales methods through frequency detuning. *Nonlinear Dynamics* 94:2919–2935. [DOI], [OA], [OAI].

Elphick, C., E. Tirapegui, M. Brachet, P. Coullet, and G. Iooss (1987). A simple global characterization for normal forms of singular vector fields. *Physica D: Nonlinear Phenomena* 29(1):95–127. [DOI].
Frangi, A., A. Colombo, A. Vizzaccaro, and C. Touzé (2025). Reduced order modelling of fully coupled electro-mechanical systems through invariant manifolds with applications to microstructures. *International Journal for Numerical Methods in Engineering* 126:e7641. [DOI], [OA].

Gabale, A. P. and S. Sinha (2009). A direct analysis of nonlinear systems with external periodic excitations via normal forms. *Nonlinear Dynamics* 55:79-93. [DOI], [OA].

Givois, A., A. Grolet, O. Thomas, and J.-F. Deü (2019). On the frequency response computation of geometrically nonlinear flat structures using reduced-order finite element models. *Nonlinear Dynamics* 97(2):1747-1781. [DOI], [HAL].

Gobat, G., L. Guillot, A. Frangi, B. Cochelin, and C. Touzé (2021). Backbone curves, Neimark-Sacker boundaries and appearance of quasi-periodicity in nonlinear oscillators: application to 1:2 internal resonance and frequency combs in MEMS. *Meccanica* 56(8):1937-1969. [DOI], [OA].

Grandi, A. A., S. Protière, and A. Lazarus (2021). Enhancing and controlling parametric instabilities in mechanical systems. *Extreme Mechanics Letters* 43:101195. [DOI], [HAL].

Guckenheimer, J. and P. Holmes (1983). *Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields*. Vol. 42. Applied Mathematical Sciences. Springer. [DOI].

Guillot, L., B. Cochelin, and C. Vergez (2019). A generic and efficient Taylor series-based continuation method using a quadratic recast of smooth nonlinear systems. *International Journal for Numerical Methods in Engineering* 119(4):261-280. [DOI], [HAL].

Haller, G. and S. Ponsioen (2016). Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction. *Nonlinear Dynamics* 86(3):1493-1534. [DOI], [OAI].

Haragus, M. and G. Iooss (2011). *Local Bifurcations, Center Manifolds, and Normal Forms in Infinite Dimensional Systems*. EDP Science & Springer.

Haro, À., M. Canadell, J.-L. Figueras, A. Luque, and J. M. Mondelo (2016). *The Parameterization Method for Invariant Manifolds. From Rigorous Results to Effective Computations*. Springer. [DOI].

Huseyin, K. and W. Zhang (2000). On the normal forms associated with high dimensional systems. *Journal of Vibration and Acoustics* 123(2):157-169. [DOI].

Iooss, G. (1988). Global characterization of the normal form for a vector field near a closed orbit. *Journal of Differential Equations* 76(1):47-76. [DOI].

Iooss, G. and M. Adelmeyer (1999). *Topics in Bifurcation Theory and Applications*. 2nd ed. World Scientific.

Jain, S. and G. Haller (2022). How to compute invariant manifolds and their reduced dynamics in high-dimensional finite-element models. *Nonlinear Dynamics* 107:1417-1450. [DOI], [OA], [ARXIV].

Jezequel, L. and C. Lamarque (1991). Analysis of non-linear dynamical systems by the normal form theory. *Journal of Sound and Vibration* 149(3):429-459. [DOI].

Jiang, D., C. Pierre, and S. Shaw (2005). Nonlinear normal modes for vibratory systems under harmonic excitation. *Journal of Sound and Vibration* 288(4-5):791-812. [DOI].

Lamarque, C.-H., C. Touzé, and O. Thomas (2012). An upper bound for validity limits of asymptotic analytical approaches based on normal form theory. *Nonlinear Dynamics* 70(3):1931-1949. [DOI], [HAL].

Leung, A. and Q. Zhang (1998). Higher order normal form and period averaging. *Journal of Sound and Vibration* 217(5):795-806. [DOI].

Leung, A. and Q. Zhang (2003). Normal form computation without central manifold reduction. *Journal of Sound and Vibration* 266(2):261-279. [DOI].

Li, M., S. Jain, and G. Haller (2022). Nonlinear analysis of forced mechanical systems with internal resonance using spectral submanifolds – Part I: Periodic response and forced response curve. *Nonlinear Dynamics* 110:1005-1043. [DOI], [OA].

Liu, X. and D. J. Wagg (2019). Simultaneous normal form transformation and model-order reduction for systems of coupled nonlinear oscillators. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* 475:20190042. [DOI], [OA], [OAI].

Llave, R. de la and F. Kogelbauer (2019). Global persistence of Lyapunov subcenter manifolds as spectral submanifolds under dissipative perturbations. *SIAM Journal on Applied Dynamical Systems* 18(4):2099-2142. [DOI], [ARXIV].
Manneville, P. (1990). *Dissipative Structures and Weak Turbulence*. Academic Press.

Martin, A., A. Opreni, A. Vizzaccaro, M. Debeurre, L. Salles, A. Frangi, O. Thomas, and C. Touzé (2023). Reduced order modeling of geometrically nonlinear rotating structures using the direct parametrisation of invariant manifolds. *Journal of Theoretical, Computational and Applied Mechanics*. [DOI], [OA].

Murdock, J. (2003). *Normal Forms and Unfoldings for Local Dynamical Systems*. Springer.

Nasir, A., N. Sims, and D. J. Wagg (2021). Direct normal form analysis of oscillators with different combinations of geometric nonlinear stiffness terms. *Journal of Applied and Computational Mechanics* 7:1167–1182. [DOI], [OA].

Nayfeh, A. H. and D. T. Mook (1995). *Nonlinear Oscillations*. John Wiley & Sons.

Neild, S. A. and D. J. Wagg (2011). Applying the method of normal forms to second-order nonlinear vibration problems. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* 467:1141–1163. [DOI], [OAI].

Neild, S. A., A. R. Champneys, D. J. Wagg, T. L. Hill, and A. Cammarano (2015). The use of normal forms for analysing nonlinear mechanical vibrations. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* 373(2051):20140404. [DOI], [OA].

Opreni, A., A. Vizzaccaro, A. Martin, A. Frangi, and C. Touzé (2022). *MORFEInvariantManifold.jl: Nonlinear model order reduction of large-dimensional finite element models using the direct parametrisation method for invariant manifolds*. GitHub repository github.com/MORFEproject/MORFEInvariantManifold.jl.

Opreni, A., A. Vizzaccaro, C. Touzé, and A. Frangi (2023). High order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to generic forcing terms and parametrically excited systems. *Nonlinear Dynamics* 111:5401–5447. [DOI], [OA].

Opreni, A., A. Vizzaccaro, A. Frangi, and C. Touzé (2021). Model order reduction based on direct normal form: application to large finite element MEMS structures featuring internal resonance. *Nonlinear Dynamics* 105:1237–1272. [DOI], [OA].

Poincaré, H. (1892). *Les méthodes nouvelles de la mécanique céleste*. Gauthiers-Villars.

Ponsioen, S., S. Jain, and G. Haller (2020). Model reduction to spectral submanifolds and forced-response calculation in high-dimensional mechanical systems. *Journal of Sound and Vibration* 488:115640. [DOI].

Ponsioen, S., T. Pedergnana, and G. Haller (2018). Automated computation of autonomous spectral submanifolds for nonlinear modal analysis. *Journal of Sound and Vibration* 420:269–295. [DOI], [ARXIV].

Rega, G., W. Lacarbonara, and A. H. Nayfeh (2000). Reduction methods for nonlinear vibrations of spatially continuous systems with initial curvature. *IUTAM Symposium on Recent Developments in Non-linear Oscillations of Mechanical Systems*. Ed. by N. Van Dao and E. J. Kreuzer. Springer, pp 235–246. [DOI].

Roberts, A. J. (2014). *Model Emergent Dynamics in Complex Systems*. SIAM series on Mathematical Modeling and Computation. [DOI].

Roberts, A. (1997). Low-dimensional modelling of dynamics via computer algebra. *Computer Physics Communications* 100(3):215–230. [DOI], [ARXIV].

Salas, A. H. and J. E. Castillo H. (2014). Exact solution to Duffing equation and the pendulum equation. *Applied Mathematical Sciences* 8(176):8781–8789. [DOI], [HAL].

Shami, Z. A., Y. Shen, C. Giraud-Audine, C. Touzé, and O. Thomas (2022). Nonlinear dynamics of coupled oscillators in 1:2 internal resonance: effects of the non-resonant quadratic terms and recovery of the saturation effect. *Meccanica* 57(11):2701–2731. [DOI], [HAL].

Shaw, S. W., C. Pierre, and E. Pescheck (1999). Modal analysis-based reduced-order models for nonlinear structures: an invariant manifold approach. *The Shock and Vibration Digest* 31(1):3–16. [OAI].

Shaw, S. and C. Pierre (1991). Non-linear normal modes and invariant manifolds. *Journal of Sound and Vibration* 150(1):170–173. [DOI], [HAL].

Shaw, S. and C. Pierre (1993). Normal modes for non-linear vibratory systems. *Journal of Sound and Vibration* 164(1):85–124. [DOI], [OAI].

Shen, Y., A. Vizzaccaro, N. Kesmia, T. Yu, L. Salles, O. Thomas, and C. Touzé (2021). Comparison
of reduction methods for finite element geometrically nonlinear beam structures. *Vibration* 4(1):175–204. [DOI], [OA].

Thomsen, J. J. (2003). *Vibrations and Stability: Advanced Theory, Analysis and Tools*. 2nd ed. Springer.

Touzé, C. (2003). *A Normal Form Approach for Non-Linear Normal Modes*. Tech. rep. Publications du LMA. [HAL].

Touzé, C. (2014). Normal form theory and nonlinear normal modes: theoretical settings and applications. *Modal Analysis of nonlinear Mechanical Systems*. Ed. by G. Kerschen. Springer, pp 75–160. [DOI], [HAL].

Touzé, C. and M. Amabili (2006). Nonlinear normal modes for damped geometrically non-linear systems: application to reduced-order modeling of harmonically forced structures. *Journal of Sound and Vibration* 298(4-5):958–981. [DOI], [HAL].

Touzé, C. and O. Thomas (2006). Non-linear behaviour of free-edge shallow spherical shells: effect of the geometry. *International Journal of Non-Linear Mechanics* 41(5):678–692. [DOI], [HAL].

Touzé, C., O. Thomas, and A. Chaigne (2004). Hardening/softening behaviour in non-linear oscillations of structural systems using non-linear normal modes. *Journal of Sound and Vibration* 273(1-2):77–101. [DOI], [HAL].

Touzé, C., C. Camier, G. Favraud, and O. Thomas (2008). Effect of imperfections and damping on the type of nonlinearity of circular plates and shallow spherical shells. *Mathematical Problems in Engineering* 2008(1):ID 678307. [DOI], [OA].

Touzé, C., A. Vizzaccaro, and O. Thomas (2021). Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques. *Nonlinear Dynamics* 105(2):1141–1190. [DOI], [ARXIV].

Vizzaccaro, A., G. Gobat, A. Frangi, and C. Touzé (2024). Direct parametrisation of invariant manifolds for non-autonomous forced systems including superharmonic resonances. *Nonlinear Dynamics* 112:6255–6290. [DOI].

Vizzaccaro, A., L. Salles, and C. Touzé (2021). Comparison of nonlinear mappings for reduced-order modeling of vibrating structures: normal form theory and quadratic manifold method with modal derivatives. *Nonlinear Dynamics* 103:3335–3370. [DOI], [HAL].

Vizzaccaro, A., A. Givois, P. Longobardi, Y. Shen, J.-F. Deü, L. Salles, C. Touzé, and O. Thomas (2020). Non-intrusive reduced order modelling for the dynamics of geometrically nonlinear flat structures using three-dimensional finite elements. *Computational Mechanics* 66(6):1293–1319. [DOI], [OA].

Vizzaccaro, A., A. Opreni, L. Salles, A. Frangi, and C. Touzé (2022). High order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to large amplitude vibrations and uncovering of a folding point. *Nonlinear Dynamics* 110(1):525–571. [DOI], [OA].

Vizzaccaro, A., Y. Shen, L. Salles, J. Blahoš, and C. Touzé (2021). Direct computation of nonlinear mapping via normal form for reduced-order models of finite element nonlinear structures. *Computer Methods in Applied Mechanics and Engineering* 384:113957. [DOI], [ARXIV].

Volvert, M. and G. Kerschen (2022). Resonant phase lags of a Duffing oscillator. *International Journal of Non-Linear Mechanics* 146:104150. [DOI], [ARXIV].

Wagg, D. J. (2022). Normal form transformations for structural dynamics: an introduction for linear and nonlinear systems. *Journal of Structural Dynamics* 1. [DOI], [OA].

Wang, W. and A. Roberts (2012). Average and deviation for slow-fast stochastic partial differential equations. *Journal of Differential Equations* 253(5):1265–1286. [DOI], [OA].

Waswa, P. M. B. and S. Redkar (2020). A direct approach for simplifying nonlinear systems with external periodic excitation using normal forms. *Nonlinear Dynamics* 99:1065–1088. [DOI].

Wiggins, S. (2003). *Introduction to Applied Nonlinear Dynamical Systems and Chaos*. 2nd ed. Springer.

Yagasaki, K. (1998). Higher-order averaging and ultra-subharmonics in forced oscillators. *Journal of Sound and Vibration* 210(4):529–553. [DOI].

Yagasaki, K. (1999). Detection of bifurcation structures by higher-order averaging for Duffing's equation. *Nonlinear Dynamics* 18:129–158. [DOI].

Yagasaki, K. and T. Ichikawa (1999). Higher-order averaging for periodically forced weakly