Failing parametrizations: what can go wrong when approximating spectral submanifolds

Invariant manifolds provide useful insights into the behavior of nonlinear dynamical systems. For conservative vibration problems, Lyapunov subcenter manifolds constitute the nonlinear extension of spectral subspaces consisting of one or more modes of the linearized system. Conversely, spectral submanifolds represent the spectral dynamics of non-conservative, nonlinear problems. While finding global invariant manifolds remains a challenge, approximations thereof can be simple to acquire and still provide an effective framework for analyzing a wide variety of problems near equilibrium solutions. This approach has been successfully employed to study both the behavior of autonomous systems and the effects of non-autonomous forcing. The current computation strategies rely on a parametrization of the invariant manifold and the reduced dynamics thereon via truncated power series. While this leads to efficient recursive algorithms, the problem itself is ambiguous, since it permits the use of various approaches for constructing the reduced system to which the invariant manifold is conjugated. Although this ambiguity is well known, it is rarely discussed and usually resolved by an ad hoc choice of method, the effects of which are mostly neglected. In this contribution, we first analyze the performance of three popular approaches for constructing the conjugate system: the graph style parametrization, the normal form parametrization, and the normal form parametrization for “near resonances.” We then show that none of them is always superior to the others and discuss the potential benefits of tailoring the parametrization to the analyzed system. As a means for illustrating the latter, we introduce an alternative strategy for constructing the reduced dynamics and apply it to two examples from the literature, which results in a significantly improved approximation quality.


Introduction
Systems of ordinary differential equations (ODEs) arise in many theoretical and practical applications of dynamical systems, e.g., in finite element, finite volume and multibody models. High-fidelity representations usually lead to high-dimensional systems of ODEs whose numerical solution is computationally expensive. Therefore, reduced-order models that contain the distinctive characteristics of the original ODE system are desirable. For linear ODEs, this issue is well studied, and methods based on spectral theory [1] are available, while for nonlinear systems this is still an active research area.
Center manifold theory [2,3] provides an important result for nonlinear problems, namely, that the stable, unstable and center spectral subspaces of the linearized system persist as invariant manifolds for the original nonlinear one. Furthermore, a closely related and particularly relevant concept, the Lyapunov subcenter manifolds [4], explores the further partitioning of the center manifold, where any eigenspace corresponding to a non-resonant oscillatory mode of the linearized system persists as an invariant manifold for the original nonlinear one. In undamped vibration problems, this offers the possibility of reducing an arbitrarily large dynamical system to a two-dimensional manifold tangent to the eigenspace of a single oscillatory mode if it is not in resonance with any of the others. The extension of the eigenmodes of the linearized to the original nonlinear system is the subject of research on nonlinear normal modes. The Rosenberg definition of nonlinear normal modes [5] is based on synchronous periodic solutions of autonomous conservative systems which correspond to closed solution trajectories in the system's phase space, in analogy to the closed trajectories in the center subspace of the linearized system. Another notion of nonlinear normal modes was proposed by Shaw and Pierre [6] for non-conservative vibration problems, which they define as invariant manifolds in analogy to the invariance property of spectral subspaces of the linearized system. On this basis, Haller and Ponsioen [7] proposed the concept of spectral submanifolds (SSMs), as unique, smoothest, invariant manifolds tangent to spectral subspaces which satisfy suitable non-resonance conditions. Applying the results of Cabré et al. [8][9][10] for general Banach spaces, they derive existence conditions for spectral submanifolds [7] which form the theoretical basis for their numerical approximation. Series expansions of spectral submanifolds (and other invariant manifolds) can be calculated by means of the classical parametrization method as introduced in [8][9][10] and consequentially studied in [11]. Practical implementations usually employ later developments based thereon, such as [12][13][14][15]. While the parametrization method of Cabré et al. provides the most popular framework for addressing the problem of approximating invariant manifolds near equilibrium solutions, there are also prior developments discussing similar ideas, e.g., [16,17].
Motivated by the major advances in invariant manifold theory, a wide variety of methods for approximating low-dimensional invariant manifolds in non-linear systems for the purpose of model order reduction have emerged in recent decades, cf. [18] for a review. Although there are versatile and accurate global approaches such as shooting methods [19], harmonic balance [20], or discretization of the invariant manifold by a mesh extended via numerical continuation [21,22], these methods are computationally expensive. Therefore, there is also research interest in developing sufficiently accurate local techniques. Many of these techniques, such as the stiffness evaluation procedure (STEP) [23], implicit condensation [24,25] and modal derivatives (MD) [26], arose in the finite element community and focus on being non-intrusive to ensure straightforward integration with existing finite element solvers. These techniques typically require a slow/fast relationship between master and enslaved coordinates, i.e., a large frequency gap, in order to achieve an accurate approximation of the invariant manifold [18]. On the other hand, many of the most accurate local techniques for approximating invariant manifolds are based on the parametrization method of Cabré et al. [8][9][10][11]13,14] and normal form computations [27][28][29].
Different parametrizations are used in the literature to compute invariant manifolds tangent to spectral subspaces (cf. [18] for a survey). On the one hand, there is the classical approach where a master spectral subspace is chosen and the enslaved coordinates are represented as a function of the master coordinates [3,[8][9][10][11]. In this case, the parametrization of the invariant manifold can be interpreted as a graph over the (left) master spectral subspace. On the other hand, there is the normal form approach [13,27,28], which is based on the theorems of Poincaré and Poincaré-Dulac and strives to provide the simplest possible expressions for the reduced dynamics on the manifold. In this case, in the absence of resonances, the reduced dynamics is linear (Poincaré theorem), while in the presence of internal resonances, some resonant/essential nonlinear terms are required (Poincaré-Dulac theorem). However, these two approaches are not the only possibilities, since in general there are infinitely many parametrizations for an invariant manifold tangent to a spectral subspace. Yet, the influence of different parametrizations on the convergence and performance of approximations of this manifold is rarely studied. In the literature, only two aspects are usually considered: the limited convergence range of the graph style parametrization, which is not able to describe folded manifolds [14,18]; and the need to consider not only exact resonances but also "near resonances" in some cases. This intends to address the issue of poor conditioning in numerical computations [11][12][13] and in the case of mechanical structures to allow for correct representation of hardening/softening behavior [27,28].
The goal of the present work is to compare the properties and performance of various parametrizations of invariant manifolds tangent to spectral subspaces by means of suitable benchmark systems. For this purpose, we use a state-of-the-art algorithm for computing invariant manifolds [13] which is based on the parametrization method of Cabré et al. [8][9][10]. Additionally, a set of appropriate error analysis tools is used for analyzing the quality of the produced approximations as needed. The current state of the algorithm as summarized in Sect. 2 permits any parametrization which allows the coefficients of the invariant manifold to be determined from the invariance condition, resulting in a certain ambiguity. Of the infinite number of possible parametrizations, initially we focus on the graph style parametrization and the normal form parametrization mentioned above, which are special cases arising from additional requirements beyond mere invariance. In Sect. 3, both parametrizations for non-resonant systems and a variant of the normal form parametrization that considers resonances are applied to approximate invariant manifolds of seven autonomous and one non-autonomous dynamical systems. These systems are specifically constructed to investigate the properties and performance of different parametrizations and approximation orders and are intended as possible benchmarks for future development of the methodology. In regard of the latter, the last two examples are also used as a framework for introducing an alternative parametrization strategy. This ultimately aims at illustrating the untapped potential in developing specialized methods for constructing the reduced system. In particular, we show that a heuristically developed approach significantly improves the approximation quality for the investigated systems. The results for the parametrizations considered in this manuscript are discussed in Sect. 4, with particular emphasis on the instances of a failure and success among the different techniques. Recommendations for further development of the methodology are derived-based thereon in Sect. 5.

Background
Spectral submanifolds, as introduced by Haller and Ponsioen [7], are invariant manifolds asymptotic to an attractive or repelling set such as a fixed point, a periodic orbit or the closure of a quasi-periodic invariant torus. Algorithms for their approximation often deal with the case of a fixed point at the origin, to which the scope of this manuscript is limited. In Sect. 2.1, we briefly summarize the relevant definitions and existence and uniqueness conditions for those cases. Procedures for computing numerical approximations based on the parametrization method of Cabré et al. [8][9][10] are described for autonomous systems in Sect. 2.2, and for non-autonomous systems in Sect. 2.3. This procedure yields fewer equations than coefficients to be determined, allowing for different parametrizations. The resulting ambiguities and possible resolutions are discussed in Sects. 2.2.2 and 2.3.2. Lastly, a suitable set of error analysis tools is explored in Sect. 2.4.

Spectral submanifolds
Starting point is a nonlinear system given by the first order differential equations [13] Bż = Az + F(z) + F ext (z, t) (1) with z ∈ R N , A, B ∈ R N ×N , F : R N → R N , ∈ R K , F ext : R N × T K → R N and time derivatives( ·) = d(·) dt . Here, B is positive definite (B 0), F(z) are the purely nonlinear autonomous terms (∇ z F| z=0 = 0) and all non-autonomous terms F ext (z, t) with frequencies 1 , . . . , K are treated as a perturbation with parameter . In the limit = 0, (1) becomes the autonomous system Bż = Az + F(z) (2) with a fixed point at the origin z * = 0. The linearization of (2) at the origin yields The generalized eigenvalue problem of (3) gives the eigenvalues λ i and the corresponding eigenvectors v i , which are either real or complex-conjugate pairs, since A, B ∈ R N ×N . The eigenspaces E i ⊂ R N are spanned by the eigenvector v i ∈ R N in the first case, and by the real and imaginary parts of the complexconjugate eigenvector pair v i =v i+1 ∈ C N in the latter case (( ·) denotes complex conjugation). These eigenspaces are invariant for the linearized system (3) and can be combined into invariant spectral subspaces via direct summation [7,13]. Notable examples of spectral subspaces are the stable, unstable and center subspaces where ⊕ denotes the direct sum of vector spaces and the index sets I s , I u and I c refer to all eigenvalues with negative, positive and zero real parts, respectively. By the center manifold theorem, there exist stable, unstable and center invariant manifolds W s , W u and W c for the nonlinear system (2) which are tangent to the respective spectral subspaces E s , E u and E c at the origin. As introduced by Haller and Ponsioen [7] for stable fixed points at the origin (i.e., dim(E c ) = dim(E u ) = 0), a spectral submanifold is the smoothest invariant manifold tangent to a master spectral subspace E ⊂ E s with index set I E and dim(E) = M < N = dim(E s ). Define the relative spectral quotient as [7] where, by some abuse of notation, all eigenvalues of (4) whose associated eigenvector is (is not) in the span of E are considered in the denominator (numerator). Then, the SSM exists under the outer non-resonance condition as a unique, at least (σ (E)+1)-times differentiable invariant manifold tangent to E at the origin [7,Theorem 3].
In practice, we are interested in computing lowdimensional SSMs for high-dimensional systems of ODEs to obtain a reduced-order model. In this case, the determination of analytical solutions for the SSM is very costly or even impossible, and numerical approximations are desired. A common approach is to determine manifolds of (prescribed) finite order n which are tangent to E at the origin and satisfy the invariance condition in a neighborhood around the origin up to order n, e.g., by using the parametrization method as described in the next section. However, since usually n < σ(E), this approximation is not unique, but depends on details of the algorithmic procedure that are not obvious and often not deliberately chosen.

SSM approximations for autonomous systems
State-of-the-art approaches for the numerical approximation of SSMs [12,13] are based on the parametrization method as introduced by Cabré et al. [8][9][10]. First, the treatment of autonomous systems is discussed in this section, which is also the basis for treating non-autonomous systems in Sect. 2.3. The procedure for deriving recursively solvable systems of linear equations to determine the parametrization coefficients is described in Sect. 2

The parametrization method
Starting point is the autonomous system (2) and a selected master spectral subspace E ⊂ E s . A common first step is then to introduce modal coordinates which transform the linearized system (3) to diagonal or Jordan canonical form [5][6][7]12]. This, however, requires the computation of all eigenvalues and eigenvectors and results in unreasonably high computation times and memory requirements for large systems [13].
To avoid these drawbacks, we adopt the approach proposed by Jain and Haller [13], where only the eigenvalues and corresponding eigenvectors that span the master spectral subspace E are needed.
The approximate invariant manifold is defined by W : C M → R N with parametrization coordinates p ∈ C M parametrizing an M-dimensional manifold in the N -dimensional phase space of (2) [13]. The reduced dynamics on this manifold iṡ with the mapping R : C M → C M . In order to make the approximation suitable for numerical computation, both mappings are described by an ansatz in the form of a multivariate polynomial using the shorthand notation for multiple Kronecker products ⊗ proposed in [13]. The problem reduces to determining the unknown coefficient matrices W i ∈ C N ×M i and R i ∈ C M×M i , for which (8)-(11) are substituted into (2) and coefficients for powers of p are compared [7,8,13,30]. Substitution into the left hand side of (2) gives On the right hand side, the nonlinear terms are expanded into Taylor series at the origin which yields the invariance equation Comparing the coefficients for powers of p ⊗i , the linear terms return the generalized eigenvalue problem (4) but only for the M eigenvalue-eigenvector pairs that make up the master spectral subspace E [31], while the result for i ≥ 2 is with and the identity matrix 1 M of dimension M. Note that since F 1 = 0, (F • W) i and therefore also C i depend on all coefficients W j of orders j < i. This results in a recursive procedure where a linear system of equations determines the i-th order coefficients W i and R i , taking all previous orders into account. The equations for the first order can be solved by choosing the master modes and their eigenvalues as coefficients where V E contains the right master eigenvectors {v i | i ∈ I E } and E is the diagonal matrix of corresponding eigenvalues {λ i | i ∈ I E }. Using the calculation rules of the Kronecker product [32], (15) is reordered to the vectorized form [13] where L i is the co-homological operator of order i induced by the linear flow [11] which is determined by the linear parts of the original and reduced systems since it depends only on A, B and R 1 . Moreover, when R 1 is diagonal, as is the case with solution (17), the cohomological operator is block diagonal with M i blocks of size N × N , which allows for parallelization of the solution of (18) [13].

Methodological ambiguities
The presented procedure introduces some ambiguities, thus the resulting coefficients W i and R i are not unique. One source of ambiguity comes from the multiple Kronecker products introduced in (10). This is easily shown by expanding an example like for M = 2, which contains redundancies since the single and double underlined terms refer to the same monomial, respectively. The number of coefficients is larger than necessary, making the solution ambiguous from a certain point of view. On the other hand, the number of independent equations increases by the same amount, and the co-homological operator is regular, except in resonance cases which are treated in the next section. The derivation of the invariance equation introduces an implicit condition that artificially resolves this artificial ambiguity. While not pretty, this is nonetheless unproblematic and additionally allows for convenient exposition with clear and simple structure of the resulting equations. Furthermore, this can be resolved in the implementation, where for a firstorder solution of the form (17), the i-th co-homological equation can be solved independently for the coefficients of any i-th order monomial or permutation thereof (i.e., column in W i ) [13]. In this case, it is sufficient to solve the i-th co-homological equation only once for all coefficients of identical monomials (e.g., p 1 p 2 and p 2 p 1 ) as, e.g., in the implementation SSM-Tool2.1 [33]. This can be avoided by augmenting the system with basis vectors from the kernel of the cohomological operator at the expense of a larger dimension of the resulting equations, see, e.g., [14,17]. Since the present manuscript contains analytic investigations of low-dimensional dynamical systems, we choose the Kronecker notation as it yields compact and structured equations of the form A more relevant source of ambiguity is the underdetermination of (16) and (18). At the i-th order, the number of equations following from the invariance equation (14) is equal to the number of unknown coefficients W i ∈ C N ×M i ; there are no equations to determine the remaining coefficients R i ∈ C M×M i . For order 1, arbitrary vectors that span the master subspace can be chosen as columns of W 1 [13], the associated matrix R 1 is in general dense. However, this is only a matter of efficiency, since it destroys the advantageous block structure, but each of these solutions describes the same dynamics. As is shown in the remainder of this manuscript, the case is different for i ≥ 2. As long as the co-homological operators are not singular, any choice of R i yields a parametrization of the manifold for which (18) determines the corresponding coefficients W i [8][9][10][11]. Based on this, many sources [11,13,14,16,34] propose setting the coefficients R i to zero whenever possible, which results in a reduced order model with linear dynamics as long as there are no resonances. However, this choice impacts the convergence of the SSM approximation. Different parametrizations lead to quite different magnitudes of the oscillations up to which the results are useful. Moreover, resonance cases where the co-homological operators are singular need to be addressed, which is done in the next section.

Treatment of resonances
In the case of resonances, the co-homological operator becomes singular [7] and additional steps must be taken to determine a parametrization of the approximate manifold and the reduced dynamics thereon. If L i is singular, the method fails unless (18) is still solvable, which is the case if and only if the right hand side is in the image of the operator. Since the coefficients R i are part of the right hand side and have yet to be determined, this boils down to the question of whether there is any choice of coefficients that ensures that this condition is satisfied. Many authors [7,11,14,34] distinguish between inner resonances, in which only master modes are involved and where this is the case, and outer resonances, where no such M-dimensional manifold exists.
The condition that the right hand side of (18) lies in the image of L i is equivalent to the projection of the right side onto its kernel vanishing. To perform this projection, a basis for the left kernel of L i is required, which can be constructed from the left eigenvectors of the two eigenvalue problems where (·) * denotes complex-conjugate transposition. Due to the definition (19a), is in the left kernel of L i , since Assuming that R 1 is semisimple and chosen to be diagonal, as in (17), the left eigenvectors and corresponding eigenvalues of R i,i are where λ E are master eigenvalues from the diagonal of R 1 and e j = 0 · · · 0 1 0 · · · 0 ∈ C M is a unit vector that is zero everywhere except for its j-th entry. There are in total M i such eigenpairs, since this is the number of possible combinations .
Any combination ( , j) with ∈ {1, . . . , M} i and j ∈ {1, . . . , N } for which is called a resonance, and I R i is the index set of all resonances of the i-th order co-homological operator. If I R i = ∅, L i is singular and its left kernel is spanned by all Furthermore, any combination ( , j) ∈ I R i with j ∈ I E is called an inner resonance where all involved eigenvalues belong to master modes, while the other case is called a (low order) outer resonance [7]. The distinction between both cases becomes clear when n * ( , j) is used to project both sides of (18) [13] In the case of an outer resonance, u * j is B-orthogonal to the master spectral subspace E spanned by the columns of W 1 and the term u * j BW 1 = 0 [31]. Since the only remaining term (e ⊗ u * j )c i is in general not zero, the i-th order invariance equation is not solvable and there is no M-dimensional SSM tangent to that master subspace. In this case, the outer resonant mode must be added to the base of E and the method restarted, which then leads to an inner resonance.
In the case of an inner resonance, the j-th resonant mode is equal to the k-th master mode, meaning λ j = λ E k , and the choice of W 1 = V E via (17) in combination with B-orthonormality of the left and right eigenvectors yields Since the Kronecker product of two unit vectors with only one nonzero entry is again a unit vector with one nonzero entry, albeit of different dimension, just returns the ( , k)-th entry of r i . Thus, every inner resonance determines one coefficient which equals one of the undetermined coefficients R i via (19e). Each inner resonance determines one coefficient of R i , and by resolving them all, (18) becomes solvable. Any coefficient in R i that is not fixed by an inner resonance can still be chosen arbitrarily, which means the number of possible parametrizations is still infinite. Since all n * ( , j) , ( , j) ∈ I R i are linearly independent, row-by-row stacking in a matrix N * i provides a basis for the left kernel of L i . This can be used to project (18) onto the entire kernel to determine the coefficients all at once. The Boolean matrix E i = N * i D i is obtained by a similar derivation as in (29) and E i in (32) assigns the same nonzero values to the same coefficients of r i as (30), and zero otherwise. This solution, where all inner resonances are resolved via (32) and any remaining coefficients R i are set to zero, is called the normal form parametrization [13].
Note that this solution of (31) is not unique, since all coefficients of R i that are not needed to resolve inner resonances can be chosen arbitrarily. Contrary to the statement in [13], the product Rather, E i is the pseudoinverse of E i ; hence, the normal form parametrization results from the least-squares solution (32) of (31).
Of the infinite number of alternative parametrizations, another distinguished choice of coefficients which ensures solvability of (18) is where the matrix U * is a row-wise stacking of all (complex conjugate transpose) left master eigenvectors via (21b) and C i is defined in (15). This parametrization results from expressing the enslaved coordinates as a function of the master ones; hence, it is termed the graph style parametrization [13]. While so far exact expressions have been assumed, the treatment of resonances requires further considerations in the case of finite precision numerical calculations. Many authors [11,13,27,28] suggest the consideration of "near resonances," where the resonance condition (25) is approximately satisfied which causes numerical inaccuracies due to poor conditioning of the co-homological operator. As a remedy, it is proposed to treat these "near resonances" like real resonances and to determine the coefficients of the reduced dynamics according to (32). However, this proposal does not include a formal definition of "near resonances" and the case in which (25) should be considered to be approximately satisfied, making the application of this modification to the SSM parametrization ambiguous.
While the description so far has been limited to the autonomous system (2), the treatment of nonautonomous terms is covered in the next section.

SSM approximations for non-autonomous systems
The notion of spectral submanifolds and their approximation based on the parametrization method can be extended to the treatment of systems with nonautonomous forcing of order O( ) as introduced in (1) [7,13,30]. For small enough > 0, the assumed hyperbolic fixed point at the origin becomes a periodic orbit or a quasi periodic invariant torus, while the invariant manifold persists tangent to this new attractor (or repellor) under the appropriate non-resonance conditions between eigenfrequencies and forcing frequencies [7]. The treatment of the general case is beyond the scope of the investigations intended in this manuscript, so in what follows we restrict ourselves to the special case of a single forcing frequency = 1 = K that turns the hyperbolic fixed point into a periodic orbit.
To further focus the investigations in this manuscript, we restrict ourselves to the simplest case of a twodimensional SSM which is tangent to the eigenspace of a complex conjugate pair of eigenvalues. Already here, further ambiguities arise which have a strong influence on the quality of the SSM approximation and whose future treatment is the prerequisite for the application of the method in more general cases. To this end, an extension of the parametrization method for systems with mono-frequent non-autonomous forcing is introduced in Sect. 2.3.1 and the resulting ambiguities are discussed in Sect. 2.3.2. As proposed in [7,13], this is the basis for computing forced response curves and backbone curves that are used later in our analysis, as described in Sects. 2.3.3 and 2.3.4, respectively.

Extended parametrization
The treatment of non-autonomous forcing as a perturbation of the autonomous system gives the ansatz [7,13] for the invariant manifold W (p, t) and its corresponding reduced dynamics R (p, t). Substitution into F ext (z, t) and Taylor expansion around = 0 gives the external forcing on the invariant manifold as Further substituting (34) and (35) into (1) and comparing powers of yields the autonomous invariance equation (14) for O( 0 ) and for O( 1 ). To proceed from here, Haller and Jain [7,13] propose expanding into Taylor series in p and considering only the zeroth order, yielding Next, the remaining non-autonomous terms are expanded into Fourier series and substituted into (38). Comparing Fourier coefficients at order k yields the determining equations with Note that there are also more sophisticated developments [15] that consider higher order expansions in p to the effect of better reproducing non-autonomous systems and successfully treating parametric resonances. Those are not discussed further, since the focus of this manuscript is on autonomous systems; the treatment of non-autonomous systems is used as a unified framework for introducing backbone curves in Sect. 2.3.4.

Methodological ambiguities
The structure of this system of equations is similar to the autonomous case in (18) and (19a): the number of equations equals the number of coefficients x 0,k ∈ C N , but also s 0,k ∈ C M must be determined. Moreover, the operator L 0,k becomes singular if there are resonances which satisfy λ j = ik for any generalized eigenvalue λ j of (A, B).
Since the fixed point of the autonomous system is assumed to be hyperbolic, there are no generalized eigenvalues λ j of (A, B) with zero real part and the resonance condition (41) never holds. An analogous approach as for the non-resonant autonomous case would be to set all s 0,k to zero, thereby eliminating any influence of the forcing on the reduced dynamics and only considering it in the manifold coefficients x 0,k . However, this makes the analysis of quantities of interest, such as forced response curves and backbone curves, based only on the reduced dynamics (34b) impossible; hence, the "near resonance" assumption is used [13].

Forced response curves
Recall that the scope of this manuscript is restricted to the simplest case of approximating a two-dimensional SSM which is tangent to the eigenspace E of a complex conjugate pair of eigenvalues λ E 1 = λ j = λ and λ E 2 = λ j+1 =λ. To account for some effect of the non-autonomous forcing in the parametrization of the reduced dynamics (34b), Jain and Haller [13] propose considering the orders ± with i ≈ Im {λ} as "near resonances" and taking them into account via the normal form parametrization. The further treatment focuses on this case, since the graph style parametrization also proposed in [13] yields identical coefficients and results for the special case of a two-dimensional SSM and mono-frequent non-autonomous forcing. An analogous procedure as in the derivation of (30) yields two nonzero coefficients Noticing p 1 =p 2 = p, the reduced dynamics can be expressed in the forṁ where due to the special structure of the reduced dynamics, resulting from the normal form parametrization procedure in Sect. 2.2, all nonzero coefficients at the n-th order can be represented by γ n andγ n , respectively. Further simplification is achieved by transformation to polar coordinates via p = ρ e iθ and introduction of the phase shift ψ = θ − t which yields with Periodic orbits with frequency correspond to fixed points of (44) that are given by the zero level set of [13,34] The zero level set of (46) in the (ρ, )-plane is called the frequency response curve. Any periodic orbit that corresponds to a point on the frequency response curve is stable, if the real parts of both eigenvalues of the Jacobian that follows from (44) are negative [13]. For constant forcing magnitude | f |, the frequency response curve provides the amplitude-frequency relationship ρ( ), which is a useful tool for illustrating and analyzing the response of system (1) to mono-frequent external forcing.

Backbone curves
Backbone curves are another useful tool for this kind of analysis. In the context of this manuscript, we use the definition from [34], where they are introduced as "the curve of maximal amplitude of the periodic response on the SSM (...) as a function of the frequency of the external forcing". Points on the backbone curve satisfy the necessary condition [34] which gives the frequency-amplitude relationship The comparison with (43) shows that the backbone curve depends only on the SSM coefficients of the autonomous system. Note that the simple expressions (46) and (49) for the forced response curve and the backbone curve are based on the assumptions that the non-autonomous forcing is mono-frequent and that a two-dimensional SSM tangent to the eigenspace of a pair of complex conjugate eigenvalues is used to reduce the system. An extension of the method for computing forced response curves based on higher dimensional SSMs is presented in [35]. However, these forced response curves have more than one peak and a generalization of the corresponding backbone curve expressions has not been provided yet.

Error analysis
Different approaches for determining the error of series approximations or constraining the domain of the approximation in order to keep the error beneath a certain threshold have been successfully applied as a-posteriori error analysis tools in the context of the parametrization method, c.f. [36][37][38][39][40][41]. Some of the simpler techniques use defects like some norm of the residual of the invariance equation evaluated at the approximate solution to estimate the truncation error; others produce validated error bounds by employing computer assisted Newton-Kantorovich-type arguments. While the former is usually less computationally expensive and more flexible with regard to different parametrizations or dynamical systems to which the invariant manifold is conjugated, the latter yields rigorous error bounds.
The literature on the computation of validated error bounds for the parametrization method mainly focuses on the normal form parametrization for non-resonant [37,38,40] or finitely resonant [39] systems. Since the focus of this manuscript falls on investigating the approximation quality of different parametrizations, the general framework presented in the aforementioned works is briefly introduced below and applied in Sect. 3.
Define w = w n n≥1 as the infinite vector of all (vectorized) invariant manifold coefficients stacked one below the other in increasing order. For a given approximation order N , this infinite vector can be split into the sum of two components w = w + w where w represents the approximate solution that consists of the first N orders of the exact solution w * complemented with an infinite string of zeros, and w represents the (unknown) tail terms. Further define the map T(w) with where in the case of (inner) resonances, an appropriate choice of the reduced dynamics r n is assumed for all orders, L † n denotes the Moore-Penrose pseudo-inverse of L n and b is a freely chosen vector in the right nullspace of L n . In the non-resonant case, L † n ≡ L −1 n denotes the exact inverse of L n and b ≡ 0 is the zero vector.
Since the (infinite) operator T(w) is defined by arranging the invariance equations for determining w in increasing order, w = w * is a fixed point. This makes it a good candidate for bounding the size of the undetermined tail terms by applying the Banach fixed-point theorem [37,40]. In particular, the infinite vector w can be interpreted as an infinite (multi-)sequence. If only cases with ∞ n=1 |w n | < ∞ are considered, these form a 1 space, which in turn endowed with the · 1 norm, yields a Banach space. From here, the aim is to show that T(w) is a contraction on some closed ball B r (w) around the approximate solution [37,38,40], which together with the norm induced metric forms a complete metric space. Based on this, the conclusion that T(w) admits a fixed point w * in B r (w) follows from the Banach fixed-point theorem. Furthermore, the radius of the ball B r (w) can be used as an upper bound for the norm of the tail terms w * −w 1 = w * 1 ≤ r , since w contains the terms of the exact solution w * up to order N and zeros thereafter. In order to show that T(w) is a contraction, it is necessary to find a closed ball B r (w) with radius r , which T(w) maps back onto itself and on which the supremum of its Fréchet derivative is strictly less than one. Assuming w ∈ B r (w), the aforementioned radius can be determined from [37,40] which shows that T maps any w ∈ B r (w) back into B r (w). Furthermore, Y > 0 and Z (r ) is a polynomial in r > 0 with only non-negative coefficients, hence 1 − Z (r ) > 0 and Z (r ) = sup DT(w) < 1, respectively. The main task in the application of this error analysis tool, i.e., the radii polynomial method [38][39][40], is to find an upper bound for Y and Z (r ) for a given approximation w and to show that an r > 0 that satisfies (52) exists. It is straightforward to calculate the bound Y in our application, since the (infinite) vector T(w) − w has only finitely many nonzero elements. The first N orders are identically zero because the map T(w) consists of the invariance equations which produce the nonzero elements for the first N orders in the approximate solution w in the first place. Furthermore, T(w) can have at most mN nonzero elements, where m is the highest order product between elements of the argument within the invariance equations.
Note that the presented reasoning presupposes that an upper bound for the tail terms is shown to be below some desired (finite) value on a closed unit (poly-)disk, i.e., p n ∈ [−1, 1], ∀n ∈ [1, M]. Therefore, the scaling of the master eigenvectors used as a solution to the firstorder invariance equation should be chosen such that the aforementioned conditions are met [38,40]. However, this is not the only perspective of this problem. An alternative that leads to the same conclusions is to take an arbitrary scaling of the master eigenvalues in W 1 , e.g., unit length, and to define a weighted 1 ν space of infinite (multi-)sequences were ∞ n=1 |w n ||ν| n < ∞ and to endow it with the corresponding weighted · ν 1 norm. The equivalence of those two perspectives and the benefits in numerical stability of the former are discussed, e.g., in [38,40].
As mentioned above, the application of the radii polynomial method is a much studied problem in the case of the normal form parametrization for non-resonant or finitely resonant systems [37][38][39][40]. However, its application to different parametrizations requires a specifically adapted algorithm in each case, since the fixed point problem T must satisfy the invariance equations that change with the method of constructing the reduced dynamics. Moreover, for systems with infinitely many resonances, such as conservative nonlinear vibration problems, the normal form parametrization can lead to terms whose order is proportional to the order of the invariance equation to which they belong. The reason for this can be outlined as follows: in the normal form parametrization for resonant systems, the nonzero coefficients in the reduced dynamics at order n have the form (30) (i.e., (e ⊗ u * j )c n ), which among other, involves a term that results from the (standard) inner product of a known constant vector with some column of These terms could propagate through the invariance equations indefinitely or up to unreasonably high orders, thus resulting in Z (r ) being a polynomial of (arbitrarily) high order. This in turn renders the radii polynomial approach impractical for such problems. In extreme cases, the additional problem of finding an upper bound for the first positive root of the infinite polynomial Z (r ) needs to be solved, which might be difficult. Furthermore, the bound Y becomes increasingly difficult to calculate, since a higher order of nonlinearity within the invariance equation means more nonzero terms in T(w), infinitely many in the extreme case. The development of a general method for constructing the required bounds for all considered parametrizations is beyond the scope of this manuscript; hence, we develop the necessary results on an individual basis for the examples in the next section.

Parametrization study for benchmark systems
In the previous section, a state-of-the-art method for approximating spectral submanifolds in the neighborhood of fixed points and periodic orbits is described, based mainly on [13], where remaining ambiguities are highlighted and discussed. Additionally, some suitable a posteriori error analysis tools are examined there. In the following subsections, several benchmark systems without exact resonances and with hyperbolic fixed points at the origin are proposed to study the performance of the method. In Sects. 3.1 through 3.6, specifically constructed autonomous systems, for which an analytic expression of the invariant manifold and the reduced dynamics thereon is known, are studied. Thereafter, two vibration problems adopted from the literature are examined. Section 3.7 deals with a conservative oscillator originally studied by Shaw and Pierre [6] with a slight modification that allows for a closed form expression of the exact solution. A non-autonomous system adopted from [13] is considered in Sect. 3.8.
While there is an infinite number of possible parametrizations for these systems, we focus on the three variants NFP-L: the normal form parametrization for systems without resonances and linear reduced dynamics, NFP-NR: the normal form parametrization for systems with "near resonances" and nonlinear reduced dynamics, GSP: the graph style parametrization with nonlinear reduced dynamics.
The consideration of "near resonances" in the case of autonomous systems is based on the observation that pairs of complex conjugate eigenvalues with "small" real parts approximately satisfy (25), where the condition "small" is not well-defined, cf. [13, Remark 3] and Sect. 2.2.3. Conversely, in the limit of vanishing real parts the resonances become exact.
In Sects. 3.1 and 3.2, two-dimensional autonomous systems with a one-dimensional slow invariant manifold are considered. The performance of NFP-L and GSP is fairly compared, where the low dimensionality of these benchmark problems allows the computation of all coefficients in terms of infinite power series as well as their domain of convergence in physical coordinates. This, combined with the inherently non-resonant nature of the systems, also allows all parametrizations to be carefully studied in terms of the truncation error.
In Sects. 3.3 through 3.6, three-dimensional autonomous systems are considered which possess a two-dimensional slow invariant manifold that is tangent to the eigenspace of a complex conjugate pair of eigenvalues at the origin. Each system is specifically chosen to favor one parametrization style over the others and to highlight unexpected behavior.
In particular, Sect. 3.3 deals with a system with planar slow SSM; hence, GSP shows good performance. In fact, for systems with a hyperbolic fixed point at the origin and no exact resonances, linear invariant manifolds can always be calculated exactly using GSP, which can be shown as follows.
A (smooth) invariant manifold W of a dynamical systemż =F(z) is characterized by the defining vector fieldF(z) being tangential to the manifold W at all points z ∈ W whereF(z = 0) [2,3]. In the case of a dynamical system of the form (2) with a linear invariant where G(p) is an M-dimensional vector valued function that depends on p. The substitution of a generic multivariate power series for both the manifold and the reduced dynamics in (2) yields Using GSP as parametrization gives Finally, substituting W(p) = W 1 p solves the invariance equation and since the solution at all orders is unique (hyperbolic fixed point at the origin and no exact resonances), it is also guaranteed to be recovered by solving recursively, starting from W 1 as defined in (17). The reduced dynamics on the exact SSM for the system in Sect. 3.4 is linear; therefore, NFP-L yields the best approximation. Furthermore, under the assumption of a system with hyperbolic fixed point at the origin and no exact resonances, i.e., unique solution of the invariance equation at all orders, NFP-L is guaranteed to recover the correct solution for any system with linear reduced dynamics by Poincaré's theorem [28]. Section 3.5 highlights an unexpected case where NFP-NR recovers the exact solution for a system with a stable limit cycle on the unstable manifold. Moreover, it is shown that the good performance of NFP-NR is not due to "small" real parts of the master eigenvalues or the conditioning of the co-homological operators, thus challenging the common practice [12,13,34] of choosing the parametrization based on linear algebra arguments.
In Sect. 3.7, a conservative oscillator originally studied by Shaw and Pierre [6] and slightly modified to permit a closed form solution is used to illustrate the significant affect that the method of constructing the reduced dynamics has on the approximation quality. In addition, the example is used as a basis for introducing an alternative parameterization strategy tailored to systems of this type that is demonstrated to improve the approximation accuracy.
The treatment of a non-autonomous system from [13] is discussed in Sect. 3.8. The focus is on the approximation of forced response curves and backbone curves for a moderate forcing amplitude. It is shown that the poor performance of the approximation of the autonomous SSM directly affects the non-autonomous case. Finally, the heuristic approach from Sect. 3.7 is adopted to this case further illustrating the potentials of developing specialized parametrizations.
In conclusion, the first two problems focus on a fair comparison between normal form and graph style parametrization. Their analysis includes information that ranges from the exact solution, trough explicit expressions for all parametrization coefficients, to a comparison of different error analysis tools with the exact truncation error. For all remaining systems, the performance of NFP-L, NFP-NR and GSP is compared. The examples are chosen to favor one parametrization to illustrate that any of the others could fail, partly in unexpected ways. This is intended to promote the idea that more thoughts should be put in the parametrization choice, which should ideally be tailored to the system at hand. A detailed discussion of the presented results is performed in Sect. 4.

2D system with quadratic SSM
Consider the autonomous systeṁ with two hyperbolic fixed points: the stable fixed point z * 1 = 0 0 at the origin and the unstable fixed point There is a one-dimensional invariant manifold that contains both fixed points and is tangent to their respective slow eigenspaces. A parametrization for this manifold is with the corresponding reduced dynamics beinġ as can be verified via direct substitution into (57). However, this parametrization is not unique and other polynomial parametrizations of finite degree can be obtained in the following way: first, the transformation ζ = P M (ξ ) with a polynomial P Q (ξ ) = Q i=0 a i ξ i of degree Q is substituted into (59); the coefficients a i are then chosen so that the right hand side of the transformed reduced dynamicsξ = −P Q (ξ ) 1 + P Q (ξ ) /P Q (ξ ) is a polynomial of degree Q + 1. This procedure yields a nonlinear system of equations for which we have no general solution. Nevertheless, we conjecture that this construction generates countably infinitely many polynomial parametrizations of finite degree, all of which are analytic at the origin, and where the corresponding reduced dynamics has Q + 1 fixed points and the immersion of the manifold covers the segment between z * 1 and z * 2 Q times. Three examples of such transformations are ζ = 2ξ + ξ 2 , ζ = −1 + ξ 2 and ζ = 1 4 (ξ + 6ξ 2 + 9ξ 3 ), where the first one gives The flow of (57), both fixed points z * 1 , z * 2 and the parametrizations (58) and (60) are depicted in Fig. 1.
The parametrization (59) describes an embedded manifold (red line) whose corresponding reduced dynamics (59) has two fixed points ζ * 1 = 0 and ζ * 2 = −1 that correspond to the two fixed points z * 1 and z * 2 , respectively. In contrast, (60) describes an immersed manifold that covers the blue line twice since z(−1 + η) = z(−1 − η). The corresponding reduced dynamics (61) has three fixed points ξ 1 = 0 and ξ 2 = −1 and ξ 3 = −2, where ξ 1 and ξ 3 correspond to z * 1 , and ξ 2 to z * 2 . The spectral submanifolds literature ( [7,12,13], etc.) does not impose any restrictions on the parametrization, where the caused indeterminacy is deemed to be irrelevant as long as the invariance equations are solvable at all orders. However, the manifold (60) forms a set of points in phase space {z(ξ ) : ξ ∈ R} which is a proper subset of the one formed by the manifold (59) {z(ζ ) : ζ ∈ R}, hence the two sets are not identical. Therefore, there are at least two distinct invariant manifolds which are differentiable infinitely many times at the origin, which technically violates the uniqueness claim in [7,Theorem 3]. Nevertheless, (58) is the only one of these finite degree polynomial parametrizations that describes an embedded manifold, so even in this case some uniqueness is preserved and a simple addition to the assumptions of the SSM existence and uniqueness theorem might suffice to resolve this issue.
The fact that the parametrization of the manifold is not unique raises the question of what approximation the procedure described in Sect. 2.2 produces starting from either one of the fixed points. Starting with the origin z * 1 , the linearization of (57) gives the eigenvalues λ 1 = −1 and λ 2 = − 7 2 with corresponding eigenvectors v 1 = 1 0 and v 2 = 0 1 . We are interested in approximating the slow spectral submanifold tangent to the master spectral subspace E = span{v 1 } belonging to the eigenvalue λ 1 with the smallest absolute value of the real part. The relative spectral quotient (6) is and since the non-resonance condition (7)
For any finite approximation order N , an error analysis tool could be employed in order to determine an estimate or a validated upper bound for the approximation error. In the following, the radii polynomial method [38][39][40] as introduced in Sect. 2.4 is employed to compute a validated upper bound for the error. In this case, the fixed-point problem for T(w) is defined as T 1 (w) = w 1 for n = 1 and as for n ≥ 2. The corresponding norm for the residual of the invariance equation evaluated at the approximate The where the · 1 norm of that infinite matrix is equivalent to the supremum of the · 1 norms of its columns Note that w k = T k (w) = const. ∀n ∈ [1, N ] by design, therefore, n > m > N ∈ N holds for all nonzero terms in DT(w) 1 . Furthermore, all pairs of columns for m > N + 1 have a smaller · 1 norm than their m = N + 1 counterparts, since the nonzero terms start at a greater n, hence, the multipliers 1 1−n and 1 7/2−n are smaller while running through the same linear combinations of w n−m,1|2 . Therefore, the only remaining task is to construct an upper bound for the · 1 norm of the "larger" column, e.g., Recall that w ∈ B r (w) and therefore w − w 1 = w 1 = ∞ N +1 |w k,1 | + |w k,2 | ≤ r . Hence, the supremum of the upper bound above can be expressed as where In this case, applying the results discussed in Sect. 2.4 to show that T(w) is a contraction on B r (w) and that a unique fixed-point w * of T(w) exists in B r (w), respectively, reduces to finding an r > 0 that satisfies Since the least (available) upper bound for the error is of interest, the radii polynomial approach [38][39][40] yields As discussed at the beginning of Sect. 2.4, simpler methods for estimating the error based on evaluating the residual of the invariance equation for a given approximate solution are also commonly used in the literature [36,38]. This approach has the benefit of relatively low computational effort and straightforward implementation, especially since in this case the residual is identical to the bound Y , i.e., it has already been calculated.
On top of that, the simplicity of this example also allows the norm of the tail terms (i.e., the approximation error) to be determined exactly, since the solution (66) for the normal form parametrization is known explicitly. The · 1 norm of the tail terms is A comparison between the validated upper bound for the approximation error by the radii polynomial method (RPM), the residual-based error estimation (REE) and the exact error (EXE), all measured w.r.t. the · 1 norm for ν = 0.1 and N = 5, is provided in Table 1.

Graph style parametrization
Since the procedure described in Sect. 2.2 does not result in a unique parametrization for the SSM approximation, a possible alternative is to use the graph style parametrization to determine the coefficients. This approach yields the coefficients for the reduced dynamics of order n ≥ 2 where the corresponding invariance equation is with the Kronecker delta δ j i . Its solution yields the invariant manifold and the corresponding reduced dynamicṡ where (·) n denotes the Pochhammer symbol. For νp ∈ − as shown in Fig. 2; hence, any finite-dimensional truncation of (81) yields a reasonable approximation of the SSM (at most) in this range (cf. Fig. 3). Nevertheless, the nonlinear reduced dynamics allows a better approximation of the stability behavior on the manifold, since sections in which the system does not converge to the evolution point z * 1 are also possible, cf. Fig. 3. As for the NFP-L approximation, an error analysis can be performed in the case of the graph style parametrization, however, there are some differences. In particular, verifying that the approximation error is below a certain threshold may not be sufficient if the reduced dynamics becomes unstable within its domain of admissible values, i.e., the unit (poly-)disk; therefore, it should also be verified that this is not the case. Moreover, since GSP produces an additional relationship between the coefficients of the invariant manifold for n ≥ 2 (must be B-orthogonal to the left master space), the corresponding fixed-point problem can be simplified. In this case, GSP results in a linear relationship between the reduced and the master coordinates and there are no tail terms in the equations corresponding to the master subspace, as it can be observed in the invariance equation (80). Therefore, the fixed-point problem used for constructing an upper bound for the approximation error/tail terms can be simplified to The index indicating the row in w k (i.e., w k,2 ) has been dropped for the sake of brevity, since only the second coordinate depends on the reduced system in a nonlinear manner. Next, the residual of the invariance equation is evaluated The Fréchet derivative can be expressed as an infinite matrix of the form with the · 1 norm and its supremum sup DT(w) 1 Once more, this result is a quadratic equation for determining an upper bound for the approximation error The series solution (81) for GSP is also explicitly known, as well as the residual for the invariance equation (86). Based thereon, a comparison between the radii polynomial method (RPM), the residual error estimate (REE) and the exact error (EXE) for ν = 0.1 and N = 5 is given in Table 2.

2D system with cubic SSM
Consider the autonomous systeṁ with three hyperbolic fixed points: the stable fixed point z * 1 = 0 0 at the origin and the unstable fixed points z * 2 = 0 3 2 and z * There is a one-dimensional invariant manifold that contains all fixed points and is tangent to their respective slow eigenspaces. A parametrization for this manifold is with the corresponding reduced dynamics beinġ as can be verified via direct substitution into (98). Linearization of (98) at the origin z * 1 gives the eigenvalues λ 1 = −1 and λ 2 = − 9 2 with corresponding eigenvectors v 1 = 1 −2 and v 2 = 0 1 . The relative spectral quotient (6) for the slow SSM tangent to E = span{v 1 } is and since the non-resonance condition (7)

Normal form parametrization
Following the procedure described in Sect. 2.2.1, the first-order invariance equation Since there are no resonances, the reduced dynamicṡ following from the normal form parametrization is linear and the manifold is given by the power series that converges for νp ∈ (−1, 1) to The transformation ζ = νp √ 1+(νp) 2 recovers (99) and, when substituted into (103), (100). The domain of convergence is depicted in Fig. 5, and the O(5) approximation is compared to the exact invariant manifold in Fig. 6.
Analogous to the previous example, the radii polynomial method [38][39][40] can be used as a tool for error analysis. The corresponding fixed point-problem is with the residual (108) This results in a cubic equation for determining an upper bound for the error, where the smallest positive root is of interest, since the least (available) upper bound is sought. In this case, the series solution (104) is explicitly known as well and a comparison between the radii polynomial method (RPM), the residual error estimate (REE) and the exact error (EXE) for ν = 0.1 and N = 5 is summarized in Table 3.

Graph style parametrization
The graph style parametrization, cf. Sect. 2.2.3, yields the manifold approximation with the reduced dynamicṡ and the radius of convergence r = 2 and (112)  approximation is compared to the exact invariant manifold in Fig. 6.
The fixed-point problem for the application of the radii polynomial method [38][39][40] is where the row-index (i.e., w k,2 ) has been dropped for the sake of brevity as in Sect. 3.1.2. The residual of the invariance equation is The Fréchet derivative can be expressed as the infinite matrix for which the · 1 norm satisfies where Again, the series solution (111) and the residual (116) are known explicitly and a comparison between the radii polynomial method (RPM), the residual error estimate (REE) and the exact error (EXE) for ν = 0.1 and N = 5 as summarized in Table 4 can be performed.

3D system with planar SSM and cubic reduced dynamics
Consider the autonomous systeṁ with the stable fixed point z * = 0 0 0 at the origin.
There is a two-dimensional invariant manifold that contains the origin and is given by with the reduced dynamics beinġ as can be verified by substitution into (121). Note that the manifold described by (122) is a two-dimensional plane embedded in the three-dimensional phase space. Linearization of (121) at the origin gives the eigenvalues λ 1,2 = − 1 2 ± i which means that the non-resonance condition (7)   the normal form parametrizations is less predictable, as discussed next.

Normal form parametrizations
Both NFP-L and NFP-NR result in presumably infinite power series expressions for which we do not have closed-form solutions. In both variants of the normal form parametrization, the embedding of the invariant manifold approximation is coincident with the eigenspace E, which is the exact solution. However, the approximations z = W(p) contain nonzero coefficients for orders greater than one because the parametrization within this plane is distorted. The deviations between the normal form parametrizations and the original system lead to significant approximation errors in the dynamics. This is illustrated in Fig. 7, where, starting from the same initial condition, the solution trajectories for the O(5) approximations of all three parametrizations are calculated by a Runge-Kutta method and compared with that of the full system (121).

3D system with cubic SSM and linear reduced dynamics
Consider the autonomous systeṁ . For 0 < a, b ∈ R, the origin z * = 0 0 0 is a stable fixed point, which is contained in the twodimensional invariant manifold with the reduced dynamicṡ The relative spectral quotient (6) for the slow SSM tangent to E = span{v 1 , v 2 } is and since the non-resonance condition (7) is satisfied for all c 1 , c 2 ∈ N, [7, Theorem 3] provides existence for a unique, at least (σ (E)+1)-times continuously differentiable SSM tangent to E at the origin. The approximation of the slow SSM by the normal form parametrizations NFP-L and NFP-NR is discussed in Sect. 3.4.1, the graph style parametrization GSP in Sect. 3.4.2.

Normal form parametrizations
For small values of a, the complex conjugate eigenvalues λ 1 and λ 2 lead to the "near resonances" described in Sect. 2.2.3. Therefore, with the justification discussed there, the NFP-NR approximation should be used to obtain better results and convergence regions than via NFP-L. However, this is a system with a hyperbolic fixed point at the origin and no exact resonances; hence, the invariance equation has a unique solution at all orders and since the reduced dynamics (127) corresponding to the exact manifold (126) is linear, NFP-L is guaranteed to recover the correct solution by Poincaré's theorem. Indeed, approximations of order 3 and higher provide this exact solution for any parameter set 0 < a < b. Yet, for this system, despite the ill-conditioning of the co-homological operators used as justification for considering "near resonances" in [13], the NFP-NR approximation is significantly worse than NFP-L. In Fig. 8, the O(5) approximation is shown. Clearly, this is a good approximation only in a small neighborhood around the origin z * , and large deviations occur quickly. As expected, the numerical simulation of (125) with an initial condition on the NFP-NR-O(5)-manifold converges rapidly toward the actual slow SSM and slowly toward the origin. Figure 9 reveals another discrepancy manifested in the backbone curves for the NFP-NR approximation computed according to Sect. 2.3.4. While the exact reduced dynamics (127) is linear, meaning that there is no amplitude dependence of its free oscillation frequency, this is not true for the nonlinear reduced dynamics of the NFP-NR approximation. Rather, with increasing approximation order, the backbone curves indicate a slight amplitude-frequency dependence that appears to converge for amplitudes corresponding to roughly z 1 ∞ < 0.7.

Graph style parametrization
The graph style parametrization yields the infinite power series which converges to , p 2 ∈ R. This is equivalent to (126) with the transformation ζ 1 = 2 √ 3 sinh 1 3 sinh −1 3 √ 3 2 p 1 , ζ 2 = p 2 + ζ 3 1 . Although the spectral submanifold is a graph over the master subspace to which the graph style parametrization eventually converges, the range of convergence is limited. As a result, finite-order GSP approximations perform similarly poorly here as they do for NFP-NR, as exemplified in Fig. 10 for the GSP approximation of order 5.
3.5 3D system with cubic SSM and stable limit cycle I Consider the autonomous systeṁ with 0 < a < b and the unstable fixed point z * = 0 0 0 at the origin. The two-dimensional manifold with the reduced dynamicṡ is invariant with respect to (133) and contains both z * and a stable limit cycle at ζ 2 1 +ζ 2 2 = 1. The linearization of (133) at the origin gives the eigenvalues λ 1,2 = a ±i and λ 3 = −b and eigenvectors v 1,2 = 1 ±i 0 and v 3 = 0 0 1 .
Since Re {λ 1,2 } > 0 and Re {λ 3 } < 0, the relative spectral quotient σ (E) according to (6) is negative and [7, Theorem 3] cannot be used to prove existence and uniqueness for the slow SSM tangent to E = span{v 1 , v 2 }. However, since E is the entire unstable spectral subspace of (133), the existence and uniqueness of the unstable invariant manifold W u , which is equivalent to the slow SSM we are looking for, follows from the center manifold theorem, cf. Sect. 2.1. Under this premise, we discuss how this slow SSM is approximated by the normal form parametrizations NFP-L and NFP-NR in Sect. 3.5.1 and by the graph style parametrization in Sect. 3.5.2.

Normal form parametrizations
Depending on the value of a, "near resonances" can be more or less well justified, with which the discussion in Sect. 2.2.3 suggests the use of either the NFP-L or the NFP-NR approximation. However, for this system, the NFP-NR approximation of order 3 and higher always yields the exact SSM and the associated cubic reduced dynamics, independent of a. Nevertheless, by increasing a, the condition number of the co-homological operator can be arbitrarily improved, thus questioning the formal justification for using the NFP-NR approximation. The NFP-NR approximation of order 5 with  Fig. 11.
The NFP-L approximation, on the other hand, is only applicable in a small neighborhood around the origin. However, since the origin is unstable and the system dynamics of interest is represented by the stationary stable limit cycle, this approximation is of limited use.

Graph style parametrization
Analogous to the above case of NFP-NR, GSP approximations of order 3 and higher yield the exact SSM and the corresponding reduced dynamics. In contrast to the behavior of GSP in the previous Sect. 3.4, the convergence region is not bounded and the cubic manifold is found globally.
Note that this behavior is not limited to the case of both stable and unstable spectral subspaces. If we choose b < 0, [7, Theorem 3] is applicable and NFP-NR, GSP still yield the correct SSM and reduced dynamics thereon, regardless of the (absolute) parameter values. The only difference is that the stable limit cycle does not attract the vast majority of trajectories any more.
3.6 3D system with cubic SSM and stable limit cycle II Consider the autonomous systeṁ . For 0 < a < b, there is an unstable fixed point z * = 0 0 0 at the origin. The two-dimensional manifold with the reduced dynamicṡ is invariant with respect to (136) and contains both z * and a stable limit cycle. The linearization of (136) at the origin gives the eigenvalues λ 1,2 = a ± i and λ 3 = −b and eigenvectors v 1,2 = 1 ±i 0 and v 3 = 0 0 1 . Analogous to the argument in Sect. 3.5, the existence and uniqueness of the slow SSM tangential to E = span{v 1 , v 2 } follows from the center manifold theorem. To calculate numerical SSM approximations, the parameters are set to a = 1 2 and b = 10 3 givinġ z = ⎡ ⎢ ⎣ The approximation of this slow SSM by the normal form parametrizations NFP-L and NFP-NR is discussed in Sect. 3.6.1 and by the graph style parametrization in Sect. 3.6.2.

Normal form parametrizations
The NFP-L approximation with linear reduced dynamics is not able to represent the limit cycle and gives poor results, analogous to the previous section. However, in contrast to (133), the NFP-NR approximation for (136) gives useful results only in a small neighborhood of the origin and is not able to describe the manifold and the reduced dynamics in the vicinity of the limit cycle.
The NFP-NR approximation of the invariant manifold degenerates rapidly and the immersion of the manifold in phase space interpenetrates itself, as shown in Fig. 12 for the approximation of order 5.

Graph style parametrization
Analogous to the results discussed in Sect. 3.4.2, the graph style parametrization yields an infinite power series that converges to the correct solution. However, the range of convergence is again limited, meaning the limit cycle cannot be represented and the GSP approximation is also of limited use.

Conservative oscillator with linear SSM
Consider the conservative oscillator given bẏ with γ > 0, that is a slightly modified version of the system studied by Shaw and Pierre [6] to simplify the exact solution. The origin is an elliptic fixed point z * = 0 0 0 that is contained in the two-dimensional invariant manifold with the corresponding reduced dynamicṡ as can be verified by substitution into (140). The linearization of (140) at the origin gives the eigenvalues λ 1,2 = ±i and λ 3,4 = ±i √ 3 with eigenvectors v 1,2 = ±i ±i 1 1 and v 3,4 = ±i √ 3 ∓i √ 3 −1 1 . Since the real parts of all eigenvalues are zero, the relative spectral quotient σ (E) according to (6) is undefined and [7, Theorem 3] cannot be used to prove the existence and uniqueness of an invariant manifold tangent to E = span{v 1 , v 2 }. However, since E is the spectral subspace corresponding to a single, non-resonant mode of a purely imaginary pair of complex conju-gate eigenvalues, the existence and uniqueness of the desired invariant manifold follows from the Lyapunov subcenter manifold theorem [4].
Note that the manifold described by (141) is a twodimensional plane embedded in the four-dimensional phase space. Therefore, GSP recovers the exact expressions (141) and (141), as discussed at the beginning of Sect. 3. Moreover, the system has (infinitely many) exact resonances; thus, the approximation via NFP-L is not possible. This specific example choice renders further investigation of those parametrizations unnecessary, thus, allowing us to focus on the interesting behavior of the NFP-NR, which ultimately leads to one of the main observations of this paper. Indeed, the following discussion of the Lyapunov subcenter manifold tangent to E = span{v 1 , v 2 } reveals the potential for an alternative normal-form-based parametrization technique for invariant manifold approximations.

Normal form parametrization
The invariance equation for this system is resonant at all odd orders greater than one, since any linear combination of the form (n + 1)λ 1 + nλ 2 = λ 1 ∀n ∈ N (143) represents a resonance between the master eigenvalues. Therefore, the reduced dynamics needs to be chosen such that it eliminates the components of the right hand side that do not lay in the range of the co-homological operator, cf. Sect. 2.2.3. Even if the parametrization is specified a priori, i.e., normal form style (32), there still remains some indeterminacy since the invariance equation at all odd orders has infinitely many solutions. This can be artificially resolved by augmenting the equations and bordering the co-homological operators with basis vectors from their left and right kernel, as described in [14,15]. However, we want to highlight how the solution choice affects the approximation quality, as this has a large impact and can be the starting point for developing better parametrizations. There are several obvious possibilities, such as the minimum · 1 and minimum · 2 norm solutions, which yield the conservative backbone curves depicted in Figs. 13 and 14, respectively.
Note that choosing two of the infinitely many solutions leads to two significantly different approximations in terms of accuracy. This is relatively surprising considering that the algorithm (32) for constructing the coefficients of the reduced dynamics from the coefficients of the manifold does not change. This directly raises the question of how to choose the "best" normal form parametrization. While we do not have a definitive answer, there is a natural approach. Since the normal form parametrization aims to produce the simplest possible reduced dynamics, the natural approach is to choose the solution at each resonance order that minimizes the coefficients for the reduced dynamics at the next resonance order. While this can be a complicated task in general, in the special case of a reduction to a single, undamped oscillatory mode, there is only one pair of purely imaginary, complex conjugate coefficients at any resonant order [7,11,14], which greatly simplifies the matter. Namely, at any resonant order the problem reduces to choosing the solution that minimizes the absolute value of the coefficient for the reduced dynamics at the next resonant order. This could also be useful approach for trying to "force" the coefficients of the reduced dynamics to converge to zero, thus, improving numerical stability at higher orders. For the sake of full disclosure, we must admit that there is no guarantee this would result in parametrization of a better quality. Yet, in this particular case, it does, as can be deduced from the conservative backbone curves produced by this parametrization, cf. Fig. 15. To illustrate the effects of parametrization on the treatment of damped, non-autonomous, mechanical systems and to continue the ideas in the last Sect. 3.7, we study the von Kármán beam introduced in [30] and treated in [13] by the method presented in Sect. 2.
The system equations result from the discretization of the one-dimensional beam continuum by ten finite elements of equal length. The beam with density 2700 kg/m 3 , length 1 m, width 100 mm and height 1 mm is fixed at one end and harmonically forced by a distributed load f (t) = A cos( t) with constant amplitude A, as depicted in Fig. 16. We use the Matlab toolbox SSMTool 2.1 by Jain et al. [33], which already provides an implementation of the discretized von Kármán beam, to study this system. To produce the results discussed therein, the forcing amplitude is set to A = 0.002 N/m, which is larger than the one used in [13], yet still small compared to the total weight force of the beam (≈ 2.65 N). Frequency response curves and backbone curves are calculated as described in Sects. 2.3.3 and 2.3.4. From these, the corresponding peak vertical deflection of the beam tip z 29 ∞ is calculated for the visualization in Fig. 17. The result of a numerical continuation method by means of the toolbox coco [22] is used as a reference approximation for the true solution. Figure 17 suggests that the backbone curves approach the true solution roughly for z 29 ∞ < 0.005, however, not for larger amplitudes. Since the backbone curves are the curves of maximal amplitude of the forced response curves, the latter show the same qualitative behavior. However, the autonomous system without forcing is sufficient for the calculation of the backbone curves, i.e., the limitation of the convergence range already results from the approximation of the autonomous SSM. Therefore, the domain of sufficient approximation quality of the backbone curves and thus also of the forced response curves can be significantly increased by choosing a better parameterization for the approximation of the autonomous SSM.
The results of the previous section indicate a heuristic approach for an improved parametrization in this case: if the damping is ignored, the von Kármán beam system has the same structure as the conservative oscillator from above, therefore, the same modifications could be applied to the normal form parametrization. Moreover, the introduction of small damping results in relatively small changes in the relationship between frequency and amplitude in a nonlinear oscillator. The idea of adopting the approach from Sect. 3.7 to the damped von Kármán beam arises naturally when the reduction to a SSM tangent to a single oscillatory mode is considered. In this case, however, there are no exact resonances and the modifications of the solution for the invariance equation, i.e., the coefficients for the invariant manifold are not directly possible, as they would brake the equality between right and left hand side. Nevertheless, there is a direct connection between adding a multiple of a basis vector for the (right) quasikernel of the co-homological operator to the solution for the invariant manifold and changing the real part of the coefficient for the reduced dynamics, since Furthermore, note that the real part of the coefficients for the reduced dynamics originates from the real parts of the master eigenvalues. Those are precisely the terms that brake the exact resonances on which the introduction of nonlinear terms under the normal form parametrization is premised. Therefore, one can argue that the real part of the coefficients in the reduced system is a hallmark of the fact that the resonances are not exact. Based on this, we (heuristically) reason that changes of the aforementioned real parts could be used to construct a parametrization with more desirable properties. Specifically, we tune the real part of those coefficients at all "near resonant" orders which results in changes to the solution of the corresponding invariance equations. This manipulation is used to minimize the imaginary part of the coefficients in the next "near resonant" order. Figure 18 shows the backbone curves that result from such a parametrization, whose coefficients are also available in [43].
The backbone curves of this alternative parametrization give a much better approximation in the same plot area compared to Fig. 17.

Discussion
The goal of the study in Sect. 3 is twofold. First, to compare the properties and performance of three parametrizations commonly used to approximate spectral submanifolds (SSMs): the normal form parametrization for systems without resonances and linear reduced dynamics (NFP-L), the normal form parametrization for systems with "near resonances" and nonlinear reduced dynamics (NFP-NR), and the graph style parametrization with nonlinear reduced dynamics (GSP). Second, to illustrate the fact that the parametrization and by extension the method for constructing the conjugate (reduced) system has a strong influence on the approximation quality and should be tailored to the analyzed system. The latter is also supported by the introduction of an alternative, heuristic procedure for constructing the reduced dynamics that is tailored to a specific system type. The exemplary systems in Sects. 3.7 and 3.8 show that the resulting parametrization yields better approximation results. The common algorithmic approach for all considered approximations is to recursively determine coefficients of a power series representation of the SSM. For the low-dimensional systems in Sects. 3.1 and 3.2, this approach yields a power series representation of the exact SSM when developed ad infinitum. However, the range of convergence of this series representation is limited, meaning that also any finite dimensional truncation of this power series can approximate the exact SSM only in a limited range, which is at best as large as this range of convergence. Although this is an expected result, there is no way to determine a priori the range of convergence for a particular parametrization. This means that a posteriori error analysis tools like [36][37][38][39][40] should be employed to determine the quality of the approximation. Therefore, we also developed the necessary results by adopting the framework of the aforementioned contributions to our specific systems and parametrizations as needed.
Moreover, the study in Sect. 3.1 shows that there can be infinitely many exact polynomial representations of an SSM. However, of all the exact parametrizations considered in the study of this system, only the lowest-order one can be embedded in phase space. The immersions of the higher-order parametrizations represent points which are all elements of the embedded manifold. However, there are multiply covered sections and inversion points at which the immersions of these higher-order parametrizations are not differentiable. Although neither the NFP-L nor the GSP parametrization leads to a finite order parametrization of the exact SSM, this observation is important for possible future developments of the methodology. Even if the algorithmic approach is improved to produce a finite order representation of the exact SSM (assuming its existence), this criterion alone does not guarantee uniqueness.
Furthermore, the investigations in Sects. 3.1 and 3.2 show limitations for the convergence range of the considered parametrizations. NFP-L can be applied to systems without resonances, in this case yielding linear reduced dynamics-the simplest form of the reduced dynamics in the sense of Poincaré's normal forms. However, this means that the range of convergence of NFP-L cannot contain other stationary solutions such as fixed points, limit cycles, or quasiperiodic tori, since these cannot be represented by the linear reduced dynamics. This limitation is evident in Sect. 3.1.1, since the range of convergence is restricted only in the direction of the second fixed point, cf. Fig. 2. Although limited by the further fixed points on the invariant manifold, the extension of the convergence region of NFP-L is even smaller and does not reach these fixed points. Thus, for strongly nonlinear systems with further stationary solutions near the origin, the possible convergence range of NFP-L is quite small, and the advantage of approximating the nonlinear invariant manifold with linear reduced dynamics is negligible compared to using the spectral subspace of the linearized system. As described above, different error analysis tools can be employed to either estimate or find a rigorous upper bound for the approximation error. However, in order to fully resolve the problem, it would be best to find a more suitable parametrization strategy.
In the graph style parametrization (GSP), the variables of the reduced dynamics p directly describe the position in the master spectral subspace E (cf., e.g., (81), (131)), which is spanned by the right master eigenvectors of the linearized system. Conversely, the non-linear terms of the approximation W(p) of the invariant manifold describe the location orthogonal to E w.r.t. the inner product u i∈ [1,M] B , · , where u i∈ [1,M] denotes all left master eigenvectors. The description as a graph maps each point on E onto a corresponding point on W(p). Therefore, this parametrization is not able to describe manifolds beyond folds, since the mapping would then no longer be one-to-one. Note that the fold comprises all points on the manifold whose tangent hyperplane is spanned by vectors v for which u i∈ [1,M] B , v = 0 holds. This can be seen in Figs. 2, 4 and 5, where in each case B ≡ 1 2 and there is only one left master eigenvector in the horizontal direction. This ultimately results in the convergence range of GSP being directly restricted by folds with tangent hyperplane in the vertical direction. Therefore, if the subspace spanned by the left eigenvectors is nearly orthogonal to E, folds can occur much sooner than in the case of coinciding left and right master subspaces. An instance of the restrictive effects this has on the approximation can be observed in Fig. 5. Moreover, the convergence range of GSP for system (57) is symmetric about the origin, although only one side is constrained by a fold. Hence, folds of the invariant manifold could affect the entire convergence range of GSP and potentially further restrict the useful range of finite dimensional approximations. In contrast, the normal form parametrization is not limited by folds, which is considered to be a major advantage of NFP-L and NFP-NR [13,18]. Figure 5 shows a case where the convergence region of NFP-L extends beyond the fold that limits GSP. However, it is limited by the other fixed points; thus, the difference in terms of the convergence range between the two parametrizations is rather small. Also note that for neither one of the systems in Sects. 3.1 and 3.2 does NFP-L manage to go beyond the fold w.r.t. the orthogonal complement of E, which is the fold one would usually be interested in. Hence, in a certain sense the NFP-L fails to deliver on the promise of correctly representing folds, even though its limiting factors, i.e., other fixed points, lay well beyond the aforementioned folds.
Further comparison between NFP-L and GSP for Sects. 3.1 and 3.2 shows that NFP-L does lead to larger convergence ranges when determining the coefficients of the power series ad infinitum. Furthermore, comparison of the fifth-order approximations also shows significantly smaller errors for NFP-L within a given domain size for the reduced dynamics. Nevertheless, the same type of comparison in Figs. 3 and 6 shows that GSP gives a better description of the qualitative behavior in a larger range. In particular, this includes the sections of the approximation for the invariant manifold that do not converge to the origin. As discussed above, NFP-L is not able to represent this behavior at all.
Another interesting observation is that, while the radii polynomial approach is expected to yield tight upper bounds for the approximation error, given tight bounds on all necessary quantities related to the corresponding fixed point problem, the residual based error estimation also performs surprisingly well. Furthermore, as discussed in Sect. 2.4 there are certain drawbacks of the otherwise rigorous and accurate radii polynomial method. Namely the computational burden, the handling of systems with infinitely many resonances and the necessity of tailoring it to the method for constructing the conjugate (reduced) system, which might be very challenging for more sophisticated parametrization schemes. Based on this and the quite reasonable accuracy of the residual based error estimation, we believe that for many systems the latter should be the error analysis tool to be tried first. Of course, in the case of insufficient accuracy more advanced techniques like the radii polynomial method could be used to refine those results.
In Sects. 3.1 and 3.2, only systems with twodimensional phase space and one-dimensional invariant manifold are investigated, since that arrangement allows to conveniently determine closed form solutions for the coefficients of the power series. This in turn facilitates calculating the radius of convergence and simplifies the illustration of the system's behavior. Since the master spectral subspace is one-dimensional, however, no internal resonances can occur. Therefore, the investigations in Sects. 3.3 to 3.6 consider systems with three-dimensional phase space and twodimensional invariant manifold whose master spectral subspace is spanned in each case by two complex conjugate eigenvectors. The differential equations of all four systems, the exact parametrizations of the invariant manifolds and the respective exact reduced dynamics thereon contain only terms of up to third order. There are no exact inner resonances in those examples, as is often the case with non-conservative systems.
The study shows that none of the three considered parametrizations is fundamentally superior to the others. Depending on the system, one or more parametriza-tions fail without a priori expectation. For instance, in Sect. 3.3, GSP provides the exact solution for the manifold and the reduced dynamics at order three, while both variants of the normal form parametrization are unable to do so within any finite approximation order. However, the limited range of validity of NFP-L and NFP-NR is unexpected, since all points of the invariant manifold lie in a plane in phase space, and the origin is the only fixed point on it. Given the exact manifold is planar and the origin is the only stationary solution on it, a large range of convergence of NFP-L might be expected. However, as shown in Fig. 7, NFP-L and NFP-NR produce significantly different approximations, which, however, do not even manage to qualitatively describe the correct solution trajectories. Nevertheless, the trajectory of the NFP-L approximation is close to the actual one near the origin, while NFP-NR seems to be more problematic. This is surprising given the high damping ratio (D ≈ 0.92) of the system. Section 3.4 shows a dynamical system with linear reduced dynamics and a cubic parametrization of the exact invariant manifold. As expected, NFP-L is able to recover this exact solution for all values 0 < a < b, although the conditioning of the co-homological operator becomes arbitrarily bad for a → 0, in which case Jain and Haller [13] recommend the use of NFP-NR. However, an application of NFP-NR leads to a poor approximation which cannot represent the exact solution with a finite order approximation. The convergence range of NFP-NR is limited, and for large amplitudes of the reduced dynamics, the approximation deviates strongly from the exact solution, as shown in Fig. 8. Moreover, the nonlinear terms in the reduced dynamics result in a slight dependence of the free oscillation frequency on the amplitude, as manifested by the backbone curves in Fig. 9. Although the invariant manifold and reduced dynamics for this system can be expressed very simply, the NFP-NR approximation of order 5 suggests much more complicated relationships and is only valid in a relatively small range. Similar to Sect. 3.2, within the (presumed) convergence range, the difference between the reduction to the eigenmode of the linearized model, and the nonlinear approximation by NFP-NR is small, making the benefit questionable compared to the effort. Section 3.5 shows the opposite behavior for NFP-L and NFP-NR. In this case, for arbitrary values of 0 < a < b the invariant manifold and the reduced dynamics thereon can be reproduced exactly by NFP-NR approximations of order three and higher, while NFP-L does not yield satisfactory results. Yet, the conditioning of the co-homological operator can be set arbitrarily close to one by changing a and b, so that the approach of Jain and Haller [13] suggests the use of NFP-L rather than NFP-NR. Nevertheless, NFP-L is guaranteed to perform poorly, since there is a stable limit cycle on the invariant manifold and that bounds the convergence range of the NFP-L approximation around the fixed point at the origin. However, the existence of this further stationary solution on the invariant manifold is not trivially evident from the definition of the dynamical system in (133). On the other hand, it is noteworthy that NFP-NR is able to accurately represent this stable limit cycle on the manifold, although only a local evolution around the fixed point at the origin is performed. This shows the great influence of the parametrization on the validity range of the approximations and the possible potential for further development of the whole methodology by improvements to the algorithm for constructing the conjugate (reduced) dynamical system. Furthermore, note that the system in Sect. 3.5 is the result of adjoining the extra pair of ODEs dz 1 dt = az 1 + z 2 − az 1 (z 2 1 + z 2 2 ) and dz 2 dt = −z 1 + az 2 − az 2 (z 2 1 + z 2 2 ), (i.e., a Stuart-Landau equation with a stable limit cycle solution z 1 = sin(t) and z 2 = cos(t)) to the linear differential equation dz 3 dt = −(3a +b)z 3 +(3a +b) sin 3 (t)+3 sin 2 (t) cos(t). This is a common approach to turn the non-autonomous system dz 3 dt = . . . with mono-frequent forcing into an autonomous system of larger dimension. Firstly, these results suggest that the commonly used linear algebra argument that "near resonant" terms should only be considered in the case of ill-conditioned cohomological operators and otherwise NFP-L should be used is not an optimal choice and further research into this decision process is needed. Secondly, it suggests that a more convenient way to treat non-autonomous forcing might be possible by appropriate construction of the conjugate system.
Like NFP-NR, GSP provides the exact manifold and reduced dynamics for approximations of order three and higher for the system in Sect. 3.5, but not for the system in Sect. 3.4. Even though for this system the exact reduced manifold can be expressed globally as a graph over the master spectral subspace, i.e., there are no folds, the convergence range of GSP is limited. The same result can be seen in Sect. 3.6. For both sys-tems, the convergence range of GSP is so small that ple is adopted to this damped system. Under the rough guideline of using the real part of the reduced dynamics at any "near resonant" order to minimize the imaginary part of the coefficients in the next "near resonant" order, as those are, in a sense, artifacts from the lack of an actual resonance, we developed the second reported parametrization. Somewhat surprisingly, this produces a better approximation whose backbone curves stay close to the correct one for significantly larger amplitudes. This further highlights the potential for further development of the methodology by using better criteria for determining the coefficients of the parametrization. However, no systematic methodology is yet available to implement this algorithmically since we only presented a heuristic that works well for two examples.

Conclusion
The goal of the present work is to compare the properties and performance of various commonly used parametrizations of invariant manifolds tangent to spectral subspaces by means of suitable benchmark systems. We present a study of three parametrizationsthe normal form parametrization without and with the consideration of "near resonances", and the graph style parametrization-which are applied to approximate invariant manifolds of seven benchmark systems. These different parametrizations are possible since the state-of-the-art methodology for parametrizing approximations of invariant manifolds permits a potentially infinite number of methods for constructing the reduced dynamics thereon. The study shows that all considered parametrizations may produce unusably bad approximations in more or less unexpected cases. There is no a priori criterion to determine which, if any, of the considered parametrizations should be used to produce a good approximation with large convergence range. While a posteriori error analysis tools can be employed to evaluate the produced approximations, as illustrated on the benchmark systems, there are still certain drawbacks in terms of computational and analytical effort. Furthermore, since the convergence domain of the parametrization, even if well estimated, still directly limits the usable range of any finite dimensional approximation, this presents a major limitation for the whole methodology.
The presented study shows two key issues that should be addressed by future developments and improvements of the methodology: (i) the lack of uniqueness of the parametrization that may lead to small convergence ranges if resolved sub-optimally, and (ii) the potential for tailoring the reduced dynamics to the analyzed system. Both issues and development needs are closely connected, yet distinct. The lack of uniqueness should be resolved by additional criteria to produce a parametrization with favorable convergence properties, since the usable range of any finite dimensional approximation cannot exceed the domain of convergence of the power series. While expanding the convergence range of the parametrization directly increases the quality of any finite dimensional approximation, it is difficult to develop a generally applicable strategy for doing so. This brings us to the idea of tailoring the reduced dynamics to the analyzed system. The last two examples in this contribution focus on outlining a simple strategy, attempting to partially achieve that goal. In this regard, they are intended to serve as a means to the end of incentivizing further developments in that direction by illustrating the possibility of significant gains in the approximation quality.
Finally, since the proposed benchmark systems are specifically designed to yield good results only for some or none of the considered parametrizations, they might also be used to evaluate the performance of any improved method for constructing the reduced dynamics.
Author contributions Both authors contributed to the study conception, design and analysis. Simulations were implemented and performed by A. K. Stoychev. The first draft of the manuscript was written by A. K. Stoychev, U. J. Römer reworked the discussion and revised the first draft. All authors read and approved the final manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL. The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Data Availability
Code to simulate all presented results is available in the KITOpen repository, https://doi.org/10.5445/IR/ 1000146477.

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