Nonlinear Model Reduction to Fractional and Mixed-Mode Spectral Submanifolds

A primary spectral submanifold (SSM) is the unique smoothest nonlinear continuation of a nonresonant spectral subspace $E$ of a dynamical system linearized at a fixed point. Passing from the full nonlinear dynamics to the flow on an attracting primary SSM provides a mathematically precise reduction of the full system dynamics to a very low-dimensional, smooth model in polynomial form. A limitation of this model reduction approach has been, however, that the spectral subspace yielding the SSM must be spanned by eigenvectors of the same stability type. A further limitation has been that in some problems, the nonlinear behavior of interest may be far away from the smoothest nonlinear continuation of the invariant subspace $E$. Here we remove both of these limitations by constructing a significantly extended class of SSMs that also contains invariant manifolds with mixed internal stability types and of lower smoothness class arising from fractional powers in their parametrization. We show on examples how fractional and mixed-mode SSMs extend the power of data-driven SSM reduction to transitions in shear flows, dynamic buckling of beams and periodically forced nonlinear oscillatory systems. More generally, our results reveal the general function library that should be used beyond integer-powered polynomials in fitting nonlinear reduced-order models to data.


I. INTRODUCTION
Reduced-order modeling seeks to identify a low-dimensional dynamical system that captures important features of a much higher (often infinite) dimensional physical system. This problem is too ambitious to have a general solution and, hence, a number of approaches with different assumptions have been under development. Recent reviews cover a number of these approaches (see, e.g., Amsallem et al. 1 and Leliévre et al. 2 ) while our focus here is specifically data-driven model reduction, which targets physical systems defined by data sets rather than equations (see, e.g., the focus issue by Ghadami and Epureanu 3 ).
Machine learning (ML) methods are broadly used for reducedorder modeling (Brunton et al., 4 Hartman and Mestha, 5 Mohamed, 6 Daniel et al., 7 and Calka et al. 8 ) but tend to produce overly complex models that lack physical interpretability and perform poorly outside their training range (Loiseau et al. 9 ). As an improvement, physics-inspired machine learning (PIML) seeks to encode physical features (such as governing equations, symmetries, and conservation laws) into the learning process (Karniadakis et al. 10 ). Remaining challenges for PIML are the handling of high-frequency features, application to large-scale examples, and the incorporation of multiphysics.
A broadly used alternative to ML is the dynamic mode decomposition (DMD) and its variants, reviewed recently by Schmid. 11 The original DMD algorithm of Schmid 12 seeks to fit a linear dynamical system to a set of observables, while the extended DMD (EDMD) of Williams et al. 13 performs such a fit to a user-specified set of functions defined on those observables. If the observables coincide with the phase space variables of the system, then DMD can be justified near fixed points as an approximation to the first-order term in the spatial expansion of the flow map at the fixed point. In the same setting, EDMD defined with respect to a set of polynomial functions serves as an approximation of a truncated linearizing transformation for a dynamical system in a neighborhood of the fixed points. For more general observables, Koopman operator theory is often invoked to view EDMD as an attempt to construct special functions of the observables that span an eigenspace of the Koopman operator (Rowley et al., 14 Budišić et al., 15 and Mauroy et al. 16 ). All this rationalizes the operational use of DMD or EDMD for continuous or discrete dynamical systems on domains where those systems exhibit linearizable dynamics. For such problems, DMD and EDMD are easily implementable, although the efficacy of the latter depends on the functions of observables used in the procedure.
A number of physical phenomena, however, are fundamentally non-linearizable over their region of interest given that they display coexisting isolated stationary states of various stability types, transitions among them, and bifurcations affecting them. No linearizing coordinate change may ever cover such a range of behaviors simultaneously, because no linear system can have more than one isolated stationary state. Most of the interest in data-driven reduced modeling stems originally from the need to describe, predict, and control coexisting states and the transitions connecting them, in problems ranging from transition to turbulence (Avila et al. 17 ) to tipping points in climate (Lenton et al. 18 ). DMD is often proposed for model reduction in such a non-linearizable setting as well, because a high-enough dimensional linear dynamical system can always be fitted with high accuracy to a small number of trajectories of a nonlinear system. Outside the linear range of the system, however, such a fit simply represents a black-box-type pattern matching with little predictive power (see, e.g., Cenedese et al. 19 for a demonstration of this on experimental data).
A logical next step toward the reduced modeling of nonlinearizable phenomena is to fit a simple, low-dimensional nonlinear differential equation to a small set of observables via regression. The resulting algorithm, the sparse identification of nonlinear dynamics (or SINDy; Brunton et al. 20 ), has been an influential tool for identifying low-dimensional systems whose nonlinearities are expected to fall in specific (typically polynomial) function classes. For unsteady dynamics with a priori unknown dimensionality and nonlinearities, however, the inherent sensitivities of this approach limit its applicability. Recently, model reduction approaches based on manifold learning have also been proposed (Floryan and Graham 21 and Farzamnik et al. 22 ). These methods can be very effective in locating candidate surfaces for model reduction operationally, without identifying the underlying dynamics of the phase space that creates those surfaces manifolds. It remains, therefore, unclear if these manifolds are robust under parameter changes or under the addition of external, time-dependent forcing.
A recent alternative to these model reduction methods is reduction to spectral submanifolds (SSMs), which were defined by Haller and Ponsioen 23 as the unique smoothest invariant manifolds of a continuous or discrete nonlinear system that are tangent and locally diffeomorphic to the spectral subspaces of the linearization of the system at a fixed point. Therefore, the internal dynamics of attracting SSMs constructed over the slowest set of eigenspaces provides a perfect, mathematically justifiable reduced-order model for a large set of trajectories. SSMs can cut through the boundaries of the domain of attraction of a fixed point and hence can carry nonlinearizable dynamics. These invariant manifolds, termed nonlinear normal modes (NNMs) at the time, were first envisioned and formally computed in seminal work by Shaw and Pierre 24 on nonlinear mechanical vibrations. The existence of such manifolds was later proved rigorously under nonresonance conditions even for infinitedimensional systems by Cabré et al., 25 who also showed the SSMs are unique in a high enough smoothness class. Automated computations of SSMs and their reduced dynamics for large finite-element problems are now available via the open source package SSMTool of Jain and Haller, 26 whereas an automated data-driven extraction of SSMs from arbitrary sets of observables can be carried out via the open source SSMLearn and fastSSM developed by Cenedese et al. 19 and Axås et al. 27 Applications of data-driven SSM-reduction in mechanical system identification and control of soft robots show strong performance of this approach in the domain of non-linearizable behavior (see Cenedese and Haller 28 , Alora et al. 29 and Kaszás et al. 30 ). Based on their mathematical construction, SSMs are also known to be robust with respect to parameter changes and the addition of external periodic or quasiperiodic forcing. The latter property enables them to predict experimentally observed forced response based solely on unforced training data (see Cenedese and Haller 28,31 ) and closed-loop response based on open-loop training data in applications to control (Alora et al. 29 ).
There have been, however, two main limitations in the underlying SSM theory that have been hindering the application of SSM reduction to some other important problems. First, the spectral subspace underlying the SSM had to be of like-mode type, i.e., spanned by eigenvectors of the same stability type, which does not hold for important transitions among exact coherent states in turbulence (see, e.g., Hof et al. 32 and Graham and Floryan 33 ). Similar, mixedmode connections are also known to exist in the dynamic buckling of structures (Abramovich 34 ) and dissipation induced instabilities of gyroscopic systems (Krechetnikov and Marsden 35 ). A further limitation has been that in problems with long transitions among different steady states, the actual transition manifold may divert substantially from the SSM defined as the smoothest nonlinear continuation of the invariant subspace E. In that case, an approximate polynomial reduced model will fail to signal the target of the transition. Beyond transitions in Couette flows (Kaszás et al. 30 ), this limitation also creates an issue with certain constructs in renormalization group theory (Li et al. 36 ), as noted by de la Llave. 37 Here, we remove both of these limitations for generic finitedimensional dynamical systems. Specifically, we construct explicit parametrizations for a significantly extended class of SSMs that also Chaos ARTICLE scitation.org/journal/cha contains invariant manifolds with mixed internal stability types and of lower smoothness class that arises from fractional powers in their parametrization. For both continuous and discrete dynamical systems, we derive these representations explicitly near hyperbolic, nonresonant, and semisimple fixed points. These three conditions may certainly be violated in small analytic examples with symmetries, but all of them will hold generically in datasets of physical systems. Importantly, resonances arising from repeated eigenvalues are not excluded by our results. Via the construct of Poincaré maps, the secondary SSMs obtained for maps imply the existence of such manifolds around periodic orbits of continuous dynamical systems.
Beyond the realm of SSM reduction, the specific fractionalpowered representations we derive for the full SSM family should be useful for other data-driven model reduction methods, such as EDMD. Indeed, all methods fitting dynamical systems to data have exclusively been using integer-powered polynomial functions. Here, we show what fractional powers will also appear generically in reduced-order models up to a given order of truncation. Using these additional fractional terms in the fit enables one to keep the order of the polynomial model relatively low and, hence, mitigate large derivatives in the model and spurious oscillations in the SSM in larger distances from the origin. We also discuss normal forms for the reduced dynamics in fractional SSMs, which have not been treated in the literature.
In Sec. II, we outline two simple, low-dimensional examples that illustrate the need for extending the prior theory of primary SSMs. In Sec. III, we state our main results for continuous dynamical systems along with several remarks on their interpretation and application. The proofs are lengthy and involve detailed calculations, which prompts us to relegate them to appendixes of a separate supplementary material. Section IV reformulates the same results for discrete dynamical systems. Section V develops a normal form theory for fractional reduced-order models on slow two-dimensional SSMs, which is the most frequent setting for dissipative mechanical systems exhibiting underdamped oscillatory behavior. Section VI discusses three examples of data-driven model reduction to fractional and mixed-mode SSMs arising in three physical problems: transition in Couette flows, dynamic buckling of a von Kármán beam, and transition in a forced nonlinear oscillator system. The first two examples involve autonomous systems, whereas the last example involves non-autonomous dynamics.

II. TWO MOTIVATING EXAMPLES
Consider the planar system of ODEṡ with the parameters a, b, c > 0. This system has two fixed points: a stable node x 0 = (0, 0) and a saddle-type fixed point at x 1 = a, b . Also note that both the x and the y axes are invariant lines. A representative phase portrait of (1) for a = b = 1 and c = 2.5 is shown in Fig. 1(a). This phase portrait mimics the geometry of a number of important transition problems in fluids, such as transitions between stationary states in Couette flows (see, e.g., Page and Kerswell 38 ) and transitions to turbulence in pipe flow (see, e.g., Skufca et al. 39 ). In the latter context, the saddle point mimics an exact coherent state (ECS) and its stable manifold models the edge of chaos, the border between laminar and turbulent behavior. An important objective of reduced-order modeling in these flows is to forecast whether specific initial conditions decay to the origin or move toward the interior of the turbulent region.

ARTICLE scitation.org/journal/cha
A Galerkin projection onto POD modes does not produce a feasible reduced-order model for this simple problem, as we show in Appendix A of the supplementary material. Linear data-driven techniques, such as DMD and its variants (Schmid 11 ), are clearly unable to capture simultaneously the two fixed points and the orbit connecting them, as was explicitly illustrated for Couette flows by Page and Kerswell. 38 The slow SSM of the origin, i.e., the unique smoothest nonlinear continuation of the slower eigenspace, is the x axis, on which the reduced dynamics isẋ = −bx. This gives a correct reduced-order model for the asymptotic behavior of all trajectories in the domain of attraction of the origin, yet fails to provide any information about the truly nonlinear dynamics. This is undesirable as both the coexisting saddle point and the separatrix M connecting it to the origin can be brought arbitrarily close to the origin by selecting small enough values for the parameters a and b.
A simple calculation reveals that any trajectory of the linear part of (1) can be written as a graph y = K |x| ac b for some choice of the parameter K ∈ R. This suggest that the manifold M containing the two fixed points may also just be ac b times differentiable at the origin, with the brackets referring to the integer part of a real number. Our first main objective in this paper is to establish the general existence and smoothness of secondary SSMs like M in nonlinear dynamical systems and obtain data-driven reduced models on such manifolds. These SSMs, however, turn out to lie outside the scope of usual equation-and data-driven approximation methods that use polynomials with integer powers. As we will see for this example in Appendix A of the supplementary material, use of fractional SSMs enables the data-driven construction for an accurate 1D reducedorder model for trajectories near the transition orbit. Additionally, we will show for the planar Couette flow in Sec. VI that the use of fractional SSMs increases the accuracy of the SSM-reduced model without the inclusion of higher-order terms in that model.
As a second example, consider the systeṁ for which M = x ∈ R 3 : x 3 = ax 2 1 is a unique, two-dimensional (2D), attracting C ∞ invariant manifold that is tangent to the x 3 = 0 spectral subspace. Within this manifold, the restricted dynamics is just that of a damped Duffing oscillator with two stable spirals and one saddle point at the origin, as shown in Fig. 1(b).
The x 3 = 0 spectral subspace contains both a stable and an unstable eigenvector of the origin, and hence the results of Cabré et al. 25 are not applicable to conclude the existence of M purely based on the spectrum of the linearized ODE at the origin. The more classic result of Irwin on pseudo-unstable manifolds (de la Llave and Wayne 40 ) does not apply for the parameter values used in Fig. 1 (see Appendix B of the supplementary material). Our second main objective in this paper is to derive reduced-order models on secondary SSMs such as M in Fig. 1(b). As we will see in Sec. VII, an important physical application of such mixed-mode SSMs is the reduced-order modeling of dynamic buckling of beams.

III. GENERALIZED SSMs NEAR FIXED POINTS OF AUTONOMOUS DYNAMICAL SYSTEMS
We consider a smooth system of n-dimensional ODEs of the formẋ and assume that its fixed point at x = 0 is hyperbolic, i.e., the spectrum spect with Re spect (A) denoting the real part of the spectrum of the linear operator A. We note that the C ∞ assumption on f can be relaxed to finite smoothness if Re spect (A) ⊂ (−∞, 0) or Re spect (A) ⊂ (0, ∞) holds [see statement (v) of Theorem 1]. We also note that our results in Sec. IV will also allow for the addition of small time-periodic forcing term to (3) (see Remark 10). We further assume that A is semisimple, i.e., the geometric and algebraic multiplicities of any potentially repeated eigenvalue of A are equal. This is generally true for physical oscillatory systems in which possible repeated eigenvalues arise from symmetry and hence each eigenvalue has its own associated oscillatory mode generated by an independent eigenvector. We finally assume that with the exception of possible 1 : 1 resonances created by repeated eigenvalues of A, there are no further resonances within spect (A), i.e., we have While these nonresonance conditions will typically hold for generic parameter configurations of typical dissipative systems, they can easily fail on simple toy examples with non-generic parameter choices. For instance, a simple 2D saddle-type fixed point with eigenvalues λ 1 = −1 and λ 2 = 1 will violate infinitely many of the conditions (5) for j = 1, m 1 = m 2 + 1 and j = 2, m 2 = m 1 + 1. The main motivation for the present paper is the data-driven modeling of dissipative physical systems in which such exceptional parameter configurations are unlikely. We stress, however, that repeated eigenvalues arising from a physical symmetry are not excluded by conditions (5).
We now select and fix an arbitrary spectral subspace E ⊂ R n of A and seek to find its nonlinear continuation in the full system (3). In practical terms, we are interested in nonlinear interactions along which all modes can be enslaved to a set of linear master modes spanning E. We assume that spect (A| E ) contains p purely real eigenvalues and q pairs of complex conjugate eigenvalues. Similarly, we assume that spect (A) − spect (A| E ) contains r purely real eigenvalues and s pairs of complex conjugate eigenvalues. In that case, after a linear change of coordinates, we can assume a partition of the x coordinate in the form with p + 2q + r + 2s = n, in which the coefficient matrix A of the linear part of system (3) takes its real Jordan canonical form with Some of the eigenvalues λ j , α k + iω k , κ and β m + iν m in the spectrum of A may be repeated, as 1 : 1 resonances are allowed by the nonresonance conditions (5). Our main result is the simplest to state in complexified coordinates. To this end, we introduce the variables and the piecewise constant functions with constant families K ± j ∈ R and O ± mj ∈ C. We also define the continuous function families where Q mk ∈ C and As seen in the proof of our upcoming main theorem, the function families (11) describe all continuous invariant graphs over E in the linear part of system (3). For any positive integer N, any vector e ∈ C N and any nonnegative integer vector k ∈ N N , we will use the notation Finally, we will use the notation for expansions involving monomials M k (u, z) of general positive (but potentially fractional) powers of u and z. These monomials will be indexed by a multi-index k ∈ N p+2q+r+2s but their order (referred to as monomial order) will be a homogeneous linear function of k, rather than k itself. We will explicitly write out the monomial order as functions in simpler settings (see, e.g., Propositions 1 and 2) in which the notation does not become too involved. With these definitions, our main result can be stated as follows. Theorem 1: (i) For any integer order K ≥ 2 of approximation, there exists a unique set of coefficients h k ∈ C r+s such that the full family W (E) of SSMs tangent to the spectral subspace E of the nonlinear system (3) can locally be written near the origin as where V(u, z) ∈ R r and E(u, z) ∈ C s are defined via formulas (11), with K j , O mj , L k , and Q mk satisfying the constraints (12) but otherwise chosen arbitrarily. (ii) The reduced dynamics on the individual members of the family W (E) are given by where f (u,z) denotes the (u, z) coordinate component of hold for all values of the indices j, k, , and m, then W (E) is unique and of class C ∞ . Otherwise, W (E) is a family of nonunique SSMs whose typical member is of smoothness class C η with Here the Int refers to the integer part of a positive number and min + refers to a minimum of the positive elements of a list of numbers that follow. Even in that case, however, there exists a unique member W ∞ (E) of the SSM family that is of class C ∞ and satisfies K ± j ≡ O ± mj = L k = Q mk = 0 for all j, k, , and m in Eq. (11).
is also C R in µ.

ARTICLE scitation.org/journal/cha
(v) If the function f only has finite smoothness, i.e., f = O |x| 2 ∈ C ρ and Re λ j has the same nonzero sign for all j = 1, . . . , N, then statements (i)-(iv) still hold for all K ≤ ρ.

(vi) If the function f is analytic in a neighborhood of the origin and
Re λ j has the same nonzero sign for all j = 1, . . . , N, then the series expansion (14) converges near the origin as K → ∞. Proof. See Appendix C of the supplementary material.
Remark 1: (Relation to primary SSMs) We note that to obtain classical stable and unstable manifolds as well as the primary (smoothest) SSMs defined by Haller and Ponsioen, 23 all constants K j , O mj , L k , and Q mk in the definitions (11) have to be set to zero and the signs of all λ j and α k involved must be the same. In contrast, the full family of invariant manifolds covered by Theorem 1 is significantly larger and includes secondary (fractional and mixed-mode) SSMs. The existing open-source SSMLearn package of Cenedese et al. 41 can be used to compute primary mixed-mode SSMs as well, even though the SSM theory originally supporting that package was only available for fixed points that are asymptotically stable in forward or backward time. This is because the integer-powered polynomial expansions for mixed-mode SSMs satisfy the same type of invariance equations as their like-mode counterparts and hence can be approximated from data in the same fashion.
Remark 2: (Relation to non-linearizable dynamics) Increasing K in the general formula (14) generally decreases the size of the open neighborhood of the origin in which this representation of the SSM family is valid. Specifically, in the limit of K → ∞, only the part of the SSM family W (E) is captured by formula (14) that carries linearizable dynamics, given that the main tool used in obtaining this asymptotic limit of this formula is the C ∞ linearization result of Sternberg. 42 In contrast, using a finite K in a data-driven identification of W (E) enables the accurate approximation of SSMs on larger domains containing non-linearizable dynamics, such as coexisting isolated compact invariant sets and transitions among them. We will illustrate this on specific examples in Secs. VI-VIII.
Remark 3: (Resonances) As noted after formula (5), 1 : 1 resonances are allowed among the eigenvalues of A by the assumptions of Theorem 1. For other resonances, the nonresonance condition (5) is substantially less restrictive than what one might be used to in the classical vibrations literature. For instance, 2λ 1 = 3λ 2 satisfies condition (5) and hence is not a resonance in our discussion unless further matches arise with the remaining eigenvalues, which has low probability in a physical system. Similarly, the case Imλ 1 = 3Imλ 2 is considered a 1 : 3 frequency-resonance in the nonlinear vibrations literature, but it is only a resonance in the sense of formula (5) if we also have the same resonance between the decay rates of the two corresponding modes, i.e., Reλ 1 = 3Reλ 2 holds simultaneously. This is again highly unlikely for a physical system, especially in a data-driven setting. If eigenvalue relations violating (5) do need to be considered for some reason, then primary, like-mode SSMs defined over spectral subspaces incorporating the resonant eigenvectors still exist by the theory of Cabré et al. 25 (see Li et al. 43,44 for specific examples). Finally, under resonances of order n k=1 m k > K, formula (14) still provides an expression for an approximately invariant manifold with invariance error of o |(u, z)| K . This follows from the fact that a finite polynomial (and hence analytic) transformation linearizing the system up to order K still exists. In a data-driven setting, such a small error in locating invariant manifolds is generally unavoidable for higher values of K even in the absence of resonances.

Remark 4: (Constraints on parameters)
The constraints on the parameters in (12) are relevant when the stability types of the eigenvectors aligned with v are the opposite of those aligned with u j or z k , in which case any non-flat invariant graph over E blows up at the origin. If contrast, if κ λ j > 0 and κ α k > 0 hold, then any nonzero K ± j and L k are allowed in expression (11). Remark 5: (Generic smoothness in the W (E) family) Statement (iii) of Theorem 1 specifies the order η of differentiability of a generic fractional SSM in the family of W (E). For such generic members of the family, all coefficients appearing in the definition (11) of the functions V and E are nonzero. Nongeneric subfamilies of the W (E) family will have higher smoothness classes equal to the integer part of one of the positive members in the list of quotients in formula (17). However, as we illustrate in Fig. 2 on a specific case, these nongeneric subfamilies form measure zero, nowhere dense subsets of the full W (E) family.
Finally, by their construction in the proof of Theorem 1, generic members of the fractional SSM family W (E) guaranteed by statement (i) of the theorem are actually always of class C ∞ away from the origin. At the origin, their Hölder smoothness class is C η,α for some 0 ≤ α < 1.

Remark 6: [Use of piecewise constant coefficients]
The use of the piecewise constant functions (10) enables us to pair up arbitrary positive and negative branches of fractional SSMs to obtain a single, continuous invariant manifold. This is illustrated in Fig. 3 for the linear part of system (3), which already exhibits the same geometry. Replacing the function K j (u j ) simply with a single constant K j is possible but limits the general family of graphs to symmetric ones, as seen in Fig. 3. The same observation is valid for the function O mj (u j ).
Remark 7: (Data-driven construction of fractional and mixed-mode SSMs) Primary mixed-mode SSMs require no detailed information on the spectrum of A outside their underlying spectral eigenspaces. For this reason, mixed-mode SSMs can be directly constructed from data featuring decaying nonlinear oscillations near E via the open source Matlab codes SSMLearn and fastSSM developed by Cenedese et al. 19 and Axås et al. 27 (see https://github.com/hallergroup). In contrast, the construction of fractional SSMs requires decay-rate information both inside and outside E, the latter of which is not immediately recognizable from decaying oscillation data near E. In Sec. IX, we mention a few data-driven strategies to obtain this spectral information external to E. Once this information is available, the construction of mixed-mode SSMs follows the same algorithm used in SSMLearn and fastSSM but with the fractionalpowered terms added to the integer-powered representation of the invariant manifolds. If, however, some of the fractional powers appearing in the formulas (11) are close to integers, then those fractional terms are advisable to omit in a data-driven construction of the SSM to avoid an overfit. All these steps in a data-driven implementation of Theorem 1 are currently under development and will appear elsewhere. In all our upcoming examples in Secs. VI-VIII, we obtain linear spectral information outside E from the governing equations but construct the fractional and mixed-mode manifolds in a data-driven fashion, using only numerically generated solutions of those governing equations.

Remark 8: (Uniqueness of invariant manifolds)
The C ∞ linearizing conjugacies of Sternberg 42 used in proving Theorem 1 are only known to be unique under the assumptions of statement (v) of the theorem (Kvalheim and Revzen 45 ). This, however, does not affect the uniqueness of W ∞ (E) or the uniqueness of the members of the full SSM family W(E). Indeed, as members of W(E) provide a well-defined foliation by invariant surfaces of a full neighborhood of the origin under a given C ∞ linearizing conjugacy, another C ∞ linearization conjugacy cannot change the elements and their smoothness in this foliation without violating their invariance. Only the specific parametrization of the fractional and primary SSMs in W(E) can change when one picks another linearizing conjugacy. This is also the reason behind the uniqueness of W (E) in statement (iii), given that conditions (16) cause the linearized system to have a unique continuous invariant manifold (E itself) that is tangent to E at the origin and has the same dimension as E. This manifold then has the same unique C ∞ preimage under any C ∞ linearizing transformation, and hence the corresponding W (E) is a unique manifold in the nonlinear system (2).
Remark 9: (Choice of the dimension of the SSM) Theorem 1 provides a hierarchical chain of SSM families, , with each family acting as a nonlinear, invariant continuation of a corresponding spectral subspace in the hierarchical chain E 1 ⊂ E 2 ⊂ . . . ⊂ E n−1 of spectral subspaces. In applications, the asymptotic dynamics are generally captured by the slowest family W (E 1 ) of these SSM families. If predictions are required for shorter time scales as well, one can gradually increase the size of the family by considering W (E k ) under increasing k. This procedure includes more and more of the transient dynamics while keeping the core asymptotic dynamics unchanged. One reason for considering W (E k ) with k > 1 in a data-driven setting is if the data show significant interaction of the first mode with higher modes up to order k. In most practical applications of SSMs with larger spectral gaps, however, including the lowest-order primary SSM suffices (see Cenedese et al. 19 and Cenedese et al. 46 ). A similar experience is emerging from applications of SSMs in model-predictive control of soft robots, where a two-dimensional SSM constructed along each spatial degree of freedom already provides a highly accurate reduced-order model (see Alora et al. 29 ).
The simplest case in applications is wherein E is 1D (p = 1, q = 0), and there are also r further non-oscillatory modes of the same stability type (s = 0). For that case, Theorem 1 takes the following more specific form.
Proposition 1: Assume that p = 1, q = 0, and s = 0 hold. Then, for any integer order K ≥ 2 of approximation, the full family W (E) of SSMs tangent to the spectral subspace E of the nonlinear system (3) can locally be written as with the arbitrary, piecewise constant functions C k (u 1 ) ∈ R r satisfying

ARTICLE scitation.org/journal/cha
Proof. See Appendix D of the supplementary material. Note that for the primary (smoothest) SSM in the W (E) family covered by Proposition 1, there exists a unique parameter vector sequence C k with k 1 ≥ 2 and k 4 = 0.
Another frequent case in practice is wherein E is 2D, carries oscillations (p = 0, q = 1), and there are also s further oscillatory modes (r = 0). For that case, Theorem 1 takes the following more specific form.
Proposition 2: Assume that p = 0, q = 1, and r = 0 hold. Then, for any integer order K ≥ 2 of approximation, the full family W (E) of SSMs tangent to the spectral subspace E of the nonlinear system (3) can locally be written as with the arbitrary constants D k ∈ C s satisfying Proof. See Appendix D of the supplementary material. Note that for the primary SSM W ∞ (E) in the W (E) family covered by Proposition 2, there exists a unique parameter vector sequence D k with |k| ≥ 2, k 5 = k 6 = 0.

IV. GENERALIZED SSMs NEAR FIXED POINTS OF DISCRETE DYNAMICAL SYSTEMS
The techniques and prior results used in the proof of Theorem 1 are also applicable with appropriate modifications to discrete dynamical systems defined by iterated mappings. This enables us to derive generalized SSM results for sampled datasets of autonomous dynamical systems and for the Poincaré maps of time-periodically forced dynamical systems. The latter extension, in turn, implies the existence of time-periodic generalized SSMs along periodic orbits of such systems. We only list the related results for maps in this section with their proofs relegated to appendixes in the supplementary material.
We consider an iterated nonlinear mapping of the form and assume that its fixed point at x = 0 is hyperbolic, i.e., the spectrum spect As in the case of continuous dynamical systems, we can also allow f to be of finite smoothness if spect (A) is fully inside or fully outside the complex unit circle (see Theorem 2). We again assume that A is semisimple and the following nonresonance conditions hold among its eigenvalues: As in the case of continuous dynamical systems that we have already treated, repeated eigenvalues arising from a physical symmetry are not excluded by conditions (23). We select and fix an arbitrary spectral subspace E ⊂ R n and assume that the spectrum of A| E contains p purely real eigenvalues and q pairs of complex conjugate eigenvalues. Similarly, we assume that the remainder of the spectrum of A contains r purely real eigenvalues and s pairs of complex conjugate eigenvalues. In that case, possibly after a linear change of coordinates, we can assume a partition of the x coordinate in the form with p + 2q + r + 2s = n, in which the coefficient matrix A of the linear part of system (21) takes its real Jordan canonical form where To simplify our upcoming formulas, we also assume that which means that the linearized mapping is orientation preserving along all of the v directions.
Using again the complexified variables z = a + ib ∈ C q and w = c + id ∈ C s , as well as the piecewise constant functions (10), we define the discrete analogues of formulas (11) as As in the case of Theorem 1, the function families (28) turn out to describe all continuous invariant graphs over E in the linear part of the mapping (21). Our main result on generalized SSMs for the full nonlinear mapping can then be stated as follows. Theorem 2: (i) For any integer order K ≥ 2 of approximation, there exists a unique set of coefficients h k ∈ C r+s such that the full family W (E) of SSMs tangent to the spectral subspace E of the nonlinear mapping (21) can locally be written as where V(u, z) ∈ R r and E(u, z) ∈ C s are defined via formulas (28), with K j , O mj , L k and Q mk satisfying the constraints (29) but otherwise chosen arbitrarily. (ii) The reduced dynamics on the individual members of the family W (E) are given by where f (u,z) denotes the (u, z) coordinate component of hold for all values of the indices j, k, and m, then W (E) is unique and of class C ∞ . Otherwise, W (E) is a family of nonunique SSMs whose typical member is of smoothness class C η with where the Int refers to the integer part of a positive number and min + refers to a minimum of the positive elements of a list of numbers that follow. Even in that case, however, there exists a unique member W ∞ (E) of the SSM family that is of class C ∞ and satisfies K ± j ≡ O ± mj = L k = Q mk = 0 for all j, k, and m in Eq. (28). (iv) If system (21) is C R in a parameter vector µ ∈ R P , then W (E) is also C R in µ. Appropriately rephrased versions of Remarks 1-8 continue to apply in the present, discrete setting. We also add the following remark that shows how Theorem 2 can be used to extend the results of Theorem 1 to small, time-periodic perturbations of the autonomous system (3).
Remark 10 (Extension to time-periodically forced continuous dynamical systems): Consider a small, temporally T-periodic For = 0, if the conditions of Theorem 1 are satisfied, then the resulting SSM family W (E) also acts as an SSM family for the period-T Poincaré map defined for the = 0 limit of system (33). By statement (iv) of Theorem 2, the SSMs of this Poincaré map are smooth functions of and hence will smoothly persist for the Poincaré map of system (33) for small > 0 parameter values. This, in turn, allows us to conclude the existence of T-periodic SSMs for the full dynamical system (33) for small enough > 0. These SSMs and their reduced dynamics can the sought in the form of a Taylor expansion in combined with a Fourier expansion in time, as explained and demonstrated for primary, non-mixed-mode SSMs by Haller and Ponsioen, 23 Breunung and Haller, 47 and Ponsioen et al. 48 Next, we also restate the appropriately modified versions of Propositions 1 and 2 for maps, whose proofs we write out for completeness in Appendix F of the supplementary material.
Proposition 3: Assume that p = 1, q = 0, and s = 0 hold. Then, for any integer order K ≥ 2 of approximation, the full family W (E) of SSMs tangent to the spectral subspace E of the nonlinear mapping (21) can locally be written as with the arbitrary, piecewise constant functions C k (u 1 ) ∈ R r satisfying Proposition 4: Assume that p = 0, q = 1, and r = 0 hold. Then, for any integer order K ≥ 2 of approximation, the full family W (E) of SSMs tangent to the spectral subspace E of the nonlinear system (3) can locally be written as

V. NORMAL FORMS ON 2D FRACTIONAL SSMs
It is often advantageous to bring the dynamics on a generalized SSM to a normal form. Normal form transformations simplify the reduced SSM dynamics (15) and (31) by removing as many terms from these expansions as possible via near-identity changes of coordinates. These transformations have been well-understood for primary SSMs on which the reduced dynamics are expressible via classical Taylor expansions (see Ponsioen et al., 49 Jain and Haller, 26 and Cenedese et al. 19,46 ). This classical normal form procedure (see, e.g., Arnold 50 and Guckenheimer and Holmes 51 ) is also applicable to mixed-mode primary SSMs, as those have no fractional-powered terms in their expansion either. For that reason, available equation-and datadriven model packages, such as SSMTool and SSMLearn, can directly be used to compute mixed-mode SSMs and normal forms of their reduced dynamics, even though the mixed-mode SSM theory described here in Secs. III and IV was not available at the time of the development of these packages. Normal forms on fractional SSMs, however, have not yet been treated in the literature and require a separate discussion here.
Each normal form transformation is a near-identity change of coordinates that simplifies a dynamical system near a fixed point up to a given polynomial order. The simplified dynamical system, however, is only topologically conjugate to the original system on a domain on which the transformation is at least continuously invertible. As a consequence, each normal form transformation further limits the ability of the normalized system to capture nonlinear dynamical features away from the origin. A clear manifestation of this limitation is the smooth linearization of Sternberg 42 used in proving Theorems 1 and 2 of this paper. Indeed, this linearization result succeeds in constructing a C ∞ normal form transformation that removes all nonlinear terms from the polynomial expansion of the right-hand side of the dynamical system. At the same time, the transformation only converges on a domain around the fixed point in which the original system exhibits linearizable behavior, completely missing all interesting dynamical phenomena that arise from the coexistence of several isolated invariant sets. For this reason, we do not advocate linearization for model reduction and only use linearization to identify the local signatures of SSMs which we then study globally without linearization.

ARTICLE scitation.org/journal/cha
More generally, we advocate a sparing use of normal form transformations for SSM-reduced models in order to preserve the predictive power of SSM-reduced dynamics on larger domains around the origin. To simplify our discussion, we only cover here the setting of 2D underdamped, oscillatory SSMs on which loworder normal forms have proven to be very effective in predicting coexisting, nontrivial steady states (see Ponsioen et al., 48,52 Jain and Haller, 26 and Cenedese et al. 19,46 ). These extended normal forms only seek a partial (as opposed to maximal possible) simplification of the dynamics. Namely, they do not remove near-resonant terms that would be in principle removable but would result in small denominators that seriously limit the domain of invertibility of the normal form transformation (see Cenedese et al. 19 for a detailed discussion for data-driven SSMs). Specifically, in the setting of Proposition 2, we obtain the following result.
Theorem 3: Assume that p = 0, q = 1 and r = 0 hold. Then: (i) For any integer order K ≥ 2 of approximation, there exists a near-identity change of coordinates that transforms the 2D reduced dynamics on each member of the SSM family W (E) to the forṁ The coordinate change is of class C ∞ if min m (β m /α 1 ) < 0 and of class C Int[minm(βm/α 1 )] if min m (β m /α 1 ) > 0. (iii) Assume further that s = 1 and the relation holds between the decay exponent β 1 outside E and the decay exponent α 1 inside E. Then in polar coordinates (r, φ) defined via the relation ξ = re iφ , a truncation of the normal form (37) that incorporates all possible terms below K = 3 (and possibly more, depending on the exact value of β 1 /α 1 ) can be written aṡ for appropriate constants A, B, P j , Q, and R j that depend on both the original full system and the specific SSM chosen.
Proof. See Appendix G of the supplementary material.

Remark 11:
In nonlinear vibration studies, the backbone curve for a given mode is the instantaneous frequency viewed as a function of the instantaneous amplitude. From the polar normal form (39), this backbone curve can be approximated up to cubic order by the formula Note that on the unique smoothest (or primary) SSM, only the first two terms of this expression can have nonzero coefficients. Remark 12: Another customary notion in nonlinear vibration studies is the instantaneous damping curve, defined as the quotient of the rate of change of the instantaneous amplitude and the instantaneous amplitude itself. From formula (39), we obtain this curve up to cubic order in the following form: Again, on the primary SSM, only the first two terms in κ(r) are nonzero; all other terms appear only on secondary SSMs. Remark 13: Formulas (40) and (41) remain valid for the slowest, 2D SSM family of a general, multi-degree-of-freedom dynamical system (n ≥ 2) as long as the relation 1 2 < β 1 α 1 < 2 holds between the decay rates of the first two modes, and all other modes decay at rates satisfying β m ≥ β 1 for m = 2, . . . , s. This is because under this relation, further factional terms are too high in order to enter the normal form truncated at cubic order.

VI. EXAMPLE 1: DATA-DRIVEN FRACTIONAL SSM FOR TRANSITIONS IN COUETTE FLOW
As a first example, we consider a plane Couette flow between two plates that move in opposite directions with velocity U, as shown in Fig. 4(a).
As already mentioned in Sec. II, Page and Kerswell 38 showed in detail the inevitable failure of linear model reduction methods to capture transitions among coexisting steady states in this flow. In contrast, Kaszás et al. 30 used data-driven SSM reduction to construct 1D and 2D models for such transitions. These models were obtained from an integer-powered polynomial approximation of the SSM. Here, we will show that a lower-order polynomial approximation to these manifolds is more accurate if we account for fractional-powered terms in the data-driven construction of the SSM.
The Navier-Stokes equations governing the evolution of the velocity field u(x, y, z, t) defined over the domain (x, y, z) ∈ = [0, where Re = Uh/ν denotes the Reynolds number with viscosity ν and p is the pressure. To eliminate the pressure from (42), we recall the Helmholtz-decomposition to write u as where is a scalar field and q is divergence free. The Leray projection of u, denoted by P[u], extracts the divergence-free component of u, i.e., returns P[u] = q. With this notation, the Leray projection of (42) becomes Denoting perturbations from the laminar solution U(y) = (y, 0, 0) by v, we can rewrite the Leray-projected Navier-Stokes equation (43) as whose linear part, where the fundamental wave numbers are defined as α x = 2π Lx and α z = 2π Lz . This turns the linearized Eq. (45) into an eigenvalue problem, with the eigenvalue denoted by λ. Using the explicit expressions of the differential operators involved (see, e.g., Schmid and Henningson 53 ), one finds that the least stable mode is independent of x and has eigenvalues and eigenfunctions of the form For this mode, a direct calculation shows that is the primary, slow SSM of the laminar state. Indeed, substituting an integer-powered polynomial expansion for an invariant manifold tangent to E 1 into the nonlinear Eq. (44) yields vanishing Taylor coefficients at all orders. In other words, such an expansion trivially converges and yields the single flat, analytic SSM shown as a dashed line in Fig. 4(b). The dynamics on this primary SSM is governed by the linear system (45), and hence no connection from the laminar base state U(y) to any nontrivial steady state can be contained in W ∞ (E 1 ). As a result, the transition to the base state to the other nontrivial steady state in system (44) must be a fractional SSM, W * (E 1 ), anchored at the base state, just as in our simple example (1). We confirm this numerically in Fig. 4(b), whose geometry bears a close similarity to that of our simple example in Fig. 1(a).
The fractional SSM W * (E 1 ) in Fig. 4(b) has limited differentiability. Yet, a data-driven global approximation for this SSM still gives satisfactory results for high enough orders of regression, as shown by Kaszás et al. 30 This is because every continuous curve can be arbitrarily closely approximated by smooth polynomials (Rudin 54 ). The required order for such an approximation, however, is generally higher than necessary and leads to larger errors further away from the origin.
Here, we will seek the globally defined manifold W * (E 1 ) in the asymptotic form (34), which has the same approximation power as polynomial expansions but also has the correct local form and correct degree of smoothness near the origin. We show that a datadriven fit of this exact fractional form for W * (E 1 ) is indeed more accurate at a given order than the one with purely integer-powered polynomials.
We use Channelflow of Gibson et al. 55 with a spectral discretization of n = O 10 5 for the nonlinear system (44). For the parameter values chosen (L x = 5π/2, L z = 4π/3), we have p = 1, q = 0, and r + s = n − 1 in terms of the notation used in Sec. III. We use the discrete formulation of our results from Sec. IV with an integer time sampling, which Kaszás et al. 30 found to give more robust reduced models than the continuous time formulation. The eigenvalues of the linearized map λ map and the linearized flow λ flow at the base state are related to each other via λ map = e λ flow .
For these parameters, formulas (34) give the fractional SSM parametrization In this expression, we have omitted the modulus of u 1 in the fractional powered terms as we are only interested in the positive branch of W * (E 1 ) along which u 1 ≥ 0 can be chosen. A linear term with k = (1, 0) also appears in expression (48) because the linear part of system (44) has not been diagonalized in the current set of coordinates.
Furthermore, Kaszás et al. 30 show that this heteroclinic connection between the base state is a graph over the square root of the energy input rate, J = √ |I|, with the energy input rate defined as Based on this, we will construct the fractional SSM W * (E 1 ) as a graph over J. Formally, this can be achieved by a change of coordinates. The observable J is a function of the state, F(u 1 , . . . , u p , a, b, v, c, d), where we recall the decomposition of (24). By the inverse function theorem, on a neighborhood of x = 0, there exists a function G such that u 1 = G(J, u 2 , . . . , u p , a, b, v, c, d), as long as Let us now introduce the new state vectorx = J,û 2 , . . . ,û p ,â,b,v, Under this change of coordinates, (3) must be transformed aṡ Since the linear part of (50) is related to the linear part of (3) through a similarity transformation, they have the same spectrum. This means, that we may view the manifolds as graphs over the variable J, and the parametrization contains the same fractional powers as the parametrization in the original phase space.
To compute W * (E 1 ) up to order 5 , we need to find all terms in (48) with integer and fractional powers up to K = 5. These powers Since log κ 5 log |λ 1 | > 5, we only have to use fractional terms arising from the eigenvalues with κ 1 , . . . , κ 4 in (48). We also see from Table I that the fractional powers log κ j log|λ 1 | and log κ j log|λ 1 | are very close to the integer powers 2 and 4, which are also present in expansion (48). While some of the near-integer spectral ratios arise due to resonances in the infinite-dimensional system, they are never exact resonances in the finite truncation of the system that we work with here. Nevertheless, to avoid sensitivity and overfit (see Remark 7), we set k 41 = k 43 = 0. Based on these considerations, a quintic-order approximation for the SSM-reduced map J(ı) → J(ı + 1) can be written as We determine the coefficients R k 1 ,k 42 ,k 44 from linear regression by minimizing the squared error  Table II for the reduced mapping (51).
In Fig. 5, we compare the predictions of the model (51) with those of the model which we obtain from a classical, integer-powered polynomial regression up to order K = 5. Note that this model is purely a numerical fit, given that the only SSM with an integer-powered expansion in the problem is W (E 1 ) = E 1 , whose expansion coefficients are all zero up to any order. As seen from Fig. 5, there is a noticeable improvement in the correct, fractional expansion of the model relative to the formal, polynomial fit to the same data setup to the same order. The polynomial fit ignores five terms in the regression that do appear in the actual reduced dynamics below quintic order. The accuracy of the polynomial fit can only be increased by adding higher-order, integer-powered polynomials, which in turn add higher derivatives and hence more sensitivity to the model.

ARTICLE
scitation.org/journal/cha  In our second example, we discuss the data-driven construction of a mixed-mode, SSM-reduced model for buckling using a discretized nonlinear von Kármán beam model. This model is a second-order approximation of geometrically exact beam theory that reproduces the bending-stretching coupling arising from bending displacements in the order of the beam thickness. We use the finite-element code of Jain et al., 56 who applied a combination of slow-fast reduction and SSM reduction to study oscillations of this beam model. In contrast, here we apply a sufficiently large longitudinal compression force and seek to derive an SSM-reduced nonlinear model for the ensuing buckling dynamics.
Schematically shown in Fig. 6, the finite element model of the von Kármán beam consists of 13 nodes, each with three degrees of freedom (DOF): axial, transversal, and rotational. By setting boundary conditions as pinned-pinned, we constrain two DOFs of the leftmost node and one DOF of the rightmost node. Therefore, the full system before model reduction has 12 elements and 72 DOFs. The beam is made of steel, with a Young's modulus E = 1.9 × 10 11 Pa, density ρ = 7850 kg/m 3 , viscous damping rate κ = 7 × 10 6 kg m/s, length L = 2 m, height h = 0, 01 m, and width w = 0.05 m. Euler's critical horizontal load for buckling is given by P N = N 2 π 2 EI/L 2 , with I denoting the area moment of inertia of the cross section of the beam. Accordingly, to generate the first buckling mode (N = 1), we apply the force A. 2D primary mixed-mode SSM Under the external force (53), the beam has three equilibria: one unstable equilibrium corresponding to purely axial displacement and two stable equilibria corresponding to the upper and lower buckled states. The first 10 eigenvalues of the linearized finite element model at the unstable equilibrium are listed in Table III, where we have already labeled the eigenvalues based on the role they will play in the construction of a mixed-mode SSM tangent to the 2D eigenspace with the 1D eigenspaces, E 1 and E 2 corresponding to the eigenvalues λ 1 = 11.06 and λ 2 = −11.10, respectively. The remaining eigenvalues are all of the form β m ± iν m with m > 4.
To obtain a reduced model that captures the dynamics of buckling, we seek reduction to a member of the SSM family W (E). Therefore, we have p = 2, q = r = 0, and s = n − 2 with n = 144 in terms of the notation used in Theorem 1. By statement (ii) Theorem 1, a typical member of this mixed-mode family is only continuous, given that All these fractional SSMs and the primary SSM, W ∞ (E), contain all three fixed points and are all tangent to the 2D stable subspaces of the two buckled fixed points, as illustrated in Fig. 6(b). For this reason, we select W ∞ (E) for constructing a reduced-order model because we can then use a classical polynomial expansion without fractional powers. By Remark 1, the SSMLearn package of Cenedese et al. 41 can be used directly to approximate the primary mixed-mode SSM in this example. We use two training trajectories whose initial conditions TABLE III. The first 10 eigenvalues of the linearized finite-element model for a buckling beam at its unstable fixed point. Our notation is relative to the 2D eigenspace E defined in (54). are obtained via a positive and a negative vertical force applied to the midpoint of the beam, respectively. We then turn off the vertical forces and record the approach of the resulting two trajectories to the upper and lower buckled state of the beam. The two trajectories obtained in this fashion are shown in Fig. 7(a), projected to a 3D coordinate space comprising the transverse displacement q 18 and velocityq 18 of the midpoint, as well as the axial displacement q 35 of right endpoint. Figure 7(b) shows the primary mixed-mode SSM, W ∞ (E), projected onto the same coordinate system, along with truncated versions of the two training trajectories. Figure 7(c) shows truncated training trajectories in model coordinates from the spectral subspace E = E 1 ⊕ E 2 . Finally, Fig. 7(d) shows two truncated test trajectories of the full mechanical model system and predictions for them by the W ∞ (E)-reduced and normalized 2D model. Following the procedure outlined by Cenedese et al. 41 for datadriven primary SSMs, we find that both the invariance error of the manifold and the accuracy of the reduced-order model for the dynamics on the manifold are minimal at a polynomial order of 11. This minimum, however, is quite weak for the invariance error and allows us to choose an order K = 7 approximation that reaches almost the same accuracy, resulting in a 0.1% normalized error in fitting W ∞ (E) to the trajectories (see Cenedese et al. 41 on computing this error). Figure 7 Fig. 3, but now labeled with respect to the modified definition (55) of the spectral subspace E. two training trajectories used in the procedure. Each training trajectory is truncated to a segment with no initial transients, which ensures that the trajectory is close to the SSM to be inferred from the data. The training trajectories are also shown in Fig. 7(c) in modal coordinates that parametrize the spectral subspace E. Finally, Fig. 7(d) shows a highly accurate prediction by the 2D seventhorder reduced model on W ∞ (E) for two test trajectories of the full system that were not used in the construction of the SSM-reduced model. We applied no normal form transformation to this model, because the two stable buckled states are far from the origin and hence are fully expected to lie outside the validity of a truncated normal form.

B. 4D primary mixed-mode SSM
The reduced-order model derived in Sec. VII A already captures the core dynamics (i.e., all stable and unstable steady states) of the bucking beam. To illustrate how enlarging the dimension of the SSM reveals more transient behavior for larger initial conditions (see Remark 9), we repeat the above procedure for a 4D primary SSM constructed over the enlarged spectral subspace Here, the newly added, 2D real eigenspace E 3 corresponds to the eigenvalue pair −0.36 ± i119.36 in Fig. 3. This means that we now have p = 2, q = 1, r = 0, and s = n − 4 with n = 144 in the application of Theorem 1. Accordingly, our notation for the eigenvalues featured in the theorem changes from that in Fig. 3 to that in Table IV. By Theorem 1, a typical fractional SSM in the W (E) family is still only continuous, given that The geometry depicted in Fig. 6(b) remains similar for the choice of E in (55) which again prompts us to choose the primary, mixedmode SSM W ∞ (E) for reduced-order modeling. This SSM, however, is now 4D and hence cannot be easily visualized.
We initialize four new training trajectories in a way that they involve oscillations from the second buckling mode as well. To generate each training trajectory, we apply two vertical forces of opposing orientations at 1/4 and 3/4 of the total length of the beam to create a static, horizontal S-shaped equilibrium state for the beam. We then turn off these forces and use the resulting trajectories shown in Fig. 8(a) for extracting the 4D mixed-mode SSM W ∞ (E) and its reduced dynamics. Figures 8(b) and 8(c) show extended model prediction accuracy that now even capture transients of the test trajectories that were not used in constructing the model.
We close by recalling that within mixed-mode SSM families, just as within like-mode SSMs, there are primary and secondary (fractional) SSMs. For instance, the SSMs computed in Secs. VII A and VII B are primary mixed-mode SSMs. They also have fractional mixed-mode counterparts, which we did not have to compute in this case because of the manifold geometry shown in Fig. 6(b).

VIII. EXAMPLE 3: FRACTIONAL AND MIXED-MODE SSMs IN A FORCED OSCILLATOR SYSTEM
Our third example involves a forced version of the two-DOF nonlinear oscillator system originally proposed by Shaw and Pierre, 24 which has been widely used to illustrate the concept of nonlinear normal modes in mechanical systems. A forced and slightly modified version of this example (shown in Fig. 9) was considered by Ponsioen et al. 52 with small forcing amplitudes for SSM-based model reduction. That system comprises two identical rigid bodies of mass m > 0 that move horizontally without friction, with their positions marked by the coordinates q 1 and q 2 . The bodies are connected to each other and two walls via identical springs of linear stiffness k > 0. The first spring also has a softening cubic nonlinearity with coefficient γ > 0, while the other two springs are subject to linear damping with coefficient c > 0. The first body is subjected to external forcing of amplitude A and frequency .
With the notation p 1 =q 1 , p 2 =q 2 , the first-order form of the equations of motion for the system is We first identify fractional SSMs near the trivial equilibrium of this system in the absence of external forcing (A = 0). As a next step, under the addition of external time-periodic forcing (A > 0), we construct a data-driven reduced model on a primary mixed-mode SSM of an unstable periodic orbit to capture transitions between coexisting periodic orbits. This second part of the example, therefore, serves as an illustration of our mixed-mode SSM results for discrete dynamical systems. where we have already applied the notation of Theorem 1 to the 2D spectral subspace E corresponding to the slower decaying pair of complex eigenvalues. Then, after a linear change of coordinates, system (57) with A = 0 will take the form This system is of the general form (3) with j = 0, k = 1, = 0, and m = 1. Therefore, formula (11) simplifies to the single complex function where Q ∈ C is a complex parameter.
In terms of the real (a, b, c, d) coordinates, the full family of 2D invariant graphs of the linear part over the E spectral subspace takes the form which defines a two-parameter-family of invariant graphs in the linearized system, parametrized by P 1 , P 2 ∈ R.
As guaranteed by statement (i) of Theorem 1, the corresponding 2D, slow SSM family W (E) in the full nonlinear system (59) can be written near the origin as By the complexification formula (9), when truncated at order K and re-indexed, this SSM family can be written in real coordinates as where C 0010 = D 0001 = 1, and all other C k and D k coefficients with |k| = 1 are zero. In all nonlinear normal mode and SSM calculations available in the literature for this model, only the k 3 = k 4 = 0 terms of the above expansion have been used and hence only the unique primary SSM, W ∞ (E), has been exploited for reduced-order modeling. This restriction to W ∞ (E) is fully justified when β/α is relatively large (β/α > K), in which case the fractional terms do not yet appear in the truncated expansion (62). For the parameter configuration used by Shaw and Pierre 24 and Ponsioen et al., 49 the eigenvalues in (58) give β α = 5.0729, which means that fractional powers only appear in the expansion for any member of the W (E) family for orders K > 5. The sixth order Taylor coefficients happen to vanish in this case, but the fractional coefficient √ a 2 + b 2 5.0729(k 3 +k 4) does appear and starts shaping the secondary SSMs and their reduced dynamics already below order 6.
To show the difference between the primary SSM W ∞ (E) and the secondary (fractional) SSMs in the full SSM family (62), we explicitly calculate the Taylor expansion of the C ∞ linearizing transformation underlying the results of Theorem 1 for the present model up to and including order 7. In Fig. 10, we show the transformed images of some of the SSMs of the linearized system (60) under the inverse of that linearizing transformation. Note the difference between the smooth W ∞ (E) [P = 0, Q = 0 in (60)] and the general fractional SSMs, which arises because fractional terms do start appearing beyond order K = 6. The figure also shows a general, decaying trajectory that follows the fractional SSM and only approaches W ∞ (E) closer to the origin. from SSMTool (see Jain and Haller 26 ), which are based on a leading-order approximation to the time-periodic perturbation of W ∞ (E). Note that for larger forcing amplitudes with A > 0.085, even the O (15) theoretical predictions based on W ∞ (E) divert from the actual periodic response.
To improve on these predictions, we focus on the forcing frequency = 1.07 and forcing amplitude A = 0.11 for which three periodic orbits coexist: a high-amplitude stable periodic orbit, a saddle-type periodic orbit, and a low-amplitude stable periodic orbit. Their stability types can be more specifically determined by computing their Floquet multipliers numerically, which we show in Fig. 12. For our analysis, the most important is the spectrum of the saddle-type periodic orbit. With the notation used in Theorem 2, the Floquet multipliers of that periodic orbits are λ 1 , λ 2 , and β ± iν, where λ 1 = 1.0835, λ 2 = 0.7726, β = −0.4132, ν = 0.6474.
To explore transitions among the coexisting periodic orbits, we consider the time-T Poincaré map associated with the forced system (57). All three periodic orbits have period T and hence all are hyperbolic fixed points of this map, as seen from their Floquet multipliers in Fig. 12. To capture connections among the fixed points, we seek to construct the 2D mixed-mode SSM family tangent to the 2D spectral subspace E that is spanned by the stable and unstable real eigenvectors of the saddle points. Of these mixedmode SSMs, we only consider the primary one, W ∞ (E), for the same reason as in the beam buckling problem [see Fig. 6(b)]. The existence and uniqueness of this W ∞ (E) follows from statement (iii) of Theorem 2, given that the Floquet multipliers of the saddle fixed point are nonresonant.
Since the Poincaré map is not known explicitly, we have no direct access to the quantities featured in Sec. III for the map formulation of our main results. For a data-driven construction, we launch four trajectories along the tangent space with initial conditions where x 0 is the location of the saddle-type fixed point, e 1 is its unstable eigenvector, e 2 is its weakest stable eigenvector, and δ = 0.01. We use these trajectories for extracting the primary mixed-mode SSM W ∞ (E) via polynomial regression. In Fig. 13, we show the iterates of the four training trajectories under the time-T map of system (57) via a projection onto the (q 1 , p 1 ) coordinate plane.
To illustrate the accuracy of this reduced-order model on W ∞ (E), we generate 20 test trajectories that are initialized near the saddle. Each initial condition has a distance of 2 × 10 −2 from the saddle point. Figure 14 shows that these test trajectories indeed closely track the extracted W ∞ (E) in a 3D projection from the 4D phase space of the Poincaré map.
We can also make predictions for the time evolution of these test trajectories using solely the reduced discrete model. The resulting trajectories obtained from the iterations of the full Poincaré map and the predictions for them using the 2D reduced dynamics match closely, as seen in Fig. 15.

IX. CONCLUSIONS
We have extended the existing theory of (primary, or smoothest) spectral submanifolds to the full family of continuous invariant manifolds that emanate from a nonresonant fixed point of a generic finite-dimensional, class C ∞ continuous or discrete dynamical system. This extended family now includes nonlinear continuations of spectral subspaces of mixed stability type and/or of finite smoothness. For the latter type of manifolds (fractional SSMs), we have derived a formal expansion with specific fractional powers that can be inferred from the spectrum of the linearized system. We have obtained similar results for discrete dynamical systems defined by iterated mappings. We have illustrated our results for Chaos ARTICLE scitation.org/journal/cha both discrete and continuous dynamical systems on the data-driven reduced-order modeling of non-linearizable behavior in three physical systems: a canonical shear flow, a buckling beam, and a forced nonlinear rigid body system. The main technical tool we rely on to prove these results is a classic linearization result of Sternberg, 42 coupled with the explicit construction of all continuous invariant manifolds through the origin of a general, multi-dimensional linear system with a hyperbolic fixed point. Neither this linearization result nor its refinements (see, e.g., Sell 57 ) are specific enough about the smoothness of the linearizing transformation, except for the class of C ∞ dynamical systems considered here, for which the linearization is also of class C ∞ . We note that the linearizing transformation only converges on a domain around the fixed point in which the original system exhibits linearizable behavior, completely missing all interesting dynamical phenomena that arise from the coexistence of several isolates invariant sets. For this reason, we do not advocate linearization for model reduction and only use linearization here to identify the function basis over which an SSMs should be more globally approximated from data. Specifically, we use truncated, finite fractional expansions for the SSMs and its reduced dynamics in our three data-driven examples. These truncated expansions continue to be well-defined outside the domain of linearizable behavior as well and hence have the potential to capture nontrivial dynamics, as we illustrate on examples.
A conceptual limitation of the results in this paper is that it assumes the underlying dynamical system to be finite dimensional. This is not an issue for numerical data sets of solid or fluid mechanics, as those are always generated from finite-dimensional discretization of their originally infinite-dimensional governing equations. Nevertheless, for experimental data obtained from a physical system, it would be desirable to have a general, infinite-dimensional extension of Sternberg's linearization results that we use. Such an extension is available for the most frequent case in practice, wherein the spectrum of the linearized system is either fully in the left-hand or the right-hand side of the complex plane, i.e., satisfies the conditions of statement (v) in Theorem 1 or 2 (see Elbialy 58 ). This extension implies that the finite-dimensional like-mode, primary and fractional SSMs we have identified infinite-dimensional limits under appropriate spectral gap conditions. As for primary mixedmode SSMs, their existence and uniqueness in infinite dimensions has been known under appropriate spectral gap conditions as long as they are also pseudo-unstable or pseudo-stable manifolds, i.e., normally hyperbolic and purely repelling or purely attracting, respectively (Elbialy 59 ). Primary pseudo-stable manifold SSMs, which also include primary like-mode SSMs, have also been proven to exist by Buza 60 specifically for the Navier-Stokes equations.
A more practical challenge for the fully data-driven application of our results is that the fractional powers in the parametrization of secondary SSMs depend on the real part of the spectrum of the linearization outside the underlying spectral subspace E. While traces of frequencies outside of E can be extracted from the frequency analysis of the available data, the real parts of the eigenvalues are more difficult to identify. An option is to treat the leading-order fractional powers as unknowns and determine them from a maximum likelihood regression. This approach shows promise in some of our preliminary numerical tests, but may also identify physically incorrect values for the fractional powers in some other examples we have considered. Another option is to estimate the real parts of the dominant eigenvalues outside E from a linear data-driven method, such as DMD or EDMD. A third option is to excite individual modes at small amplitudes and infer their decay rates by fitting one-degreeof-freedom linear oscillators to their forced response curves. The detailed implementation of these strategies will be pursued in other publications.
The backbone curves and damping curves derived in Sec. IV from a fractional normal form should be helpful in explaining some unexpected backbone curves and instantaneous damping curves observed in detailed experimental data from gravity measurements and hydrogel oscillations (Cenedese 61 and Eriten 62 ). In those experiments, the observed shapes of backbone curves and instantaneous damping curves cannot be explained without including the fractional-powered terms in formulas (40) and (41), which provided the original motivation of our present study. These experimental datasets along with others will be analyzed elsewhere with the methods developed here.
We close by noting that the fractional polynomial powers identified in Theorems 1 and 2 generically arise in the reduced dynamics of any non-unique invariant manifold emanating from hyperbolic fixed points. This implies that even other data-driven modeling approaches ignoring the existence of SSMs (e.g., EDMD and SINDy) would benefit from including such fractional powers in the dictionary of functions that they seek to fit to data.

SUPPLEMENTARY MATERIAL
See supplementary material for technical details for Example 1 (Appendix A) and Example 2 (Appendix B). We also prove Theorem 1 (Appendix C), Propositions 1 and 2 (Appendix D), Theorem 2 (Appendix E), Propositions 3 and 4 (Appendix F), and Theorem 3 (Appendix G).