We extend the theory of spectral submanifolds (SSMs) to general non-autonomous dynamical systems that are either weakly forced or slowly varying. Examples of such systems arise in structural dynamics, fluid–structure interactions, and control problems. The time-dependent SSMs we construct under these assumptions are normally hyperbolic and hence will persist for larger forcing and faster time dependence that are beyond the reach of our precise existence theory. For this reason, we also derive formal asymptotic expansions that, under explicitly verifiable nonresonance conditions, approximate SSMs and their aperiodic anchor trajectories accurately for stronger, faster, or even temporally discontinuous forcing. Reducing the dynamical system to these persisting SSMs provides a mathematically justified model- reduction technique for non-autonomous physical systems whose time dependence is moderate either in magnitude or speed. We illustrate the existence, persistence, and computation of temporally aperiodic SSMs in mechanical examples under chaotic forcing.

Reduced models for complex physical systems are of growing interest in various areas of applied science and engineering. Mathematically justifiable reduction approaches, such as spectral submanifold (or SSM) reduction, seek to identify the internal dynamics on lower-dimensional, attracting invariant sets in the phase space of the system. These reduced dynamics then become viable reduced models for general trajectories that approach the invariant set and synchronize with its inner motions. SSM reduction also allows for periodic or quasiperiodic time dependence in the full system but has been inapplicable to systems with more general time dependence, such as impulsive, chaotic, or discontinuous forcing. This has hindered applications of SSM reduction to a number of problems in structural dynamics. Here, we remove this limitation by extending SSM theory of temporally aperiodic dynamical systems. We obtain exact results for cases of smooth small or smooth slow forcing, but find that our formulas for SSM-reduced dynamics extend to larger and faster forcing in physical examples, including even discontinuous chaotic forcing.

## I. INTRODUCTION

In its simplest form, a spectral submanifold (SSM) of an autonomous dynamical system is an invariant manifold $W(E)$ that is tangent to a spectral subspace $E$ of the linearized system at a fixed point (Haller and Ponsioen^{1}). Classic examples of spectral submanifolds are the stable, unstable, and center manifolds tangent to spectral subspaces in which the linearized spectrum has eigenvalues with purely negative, positive, and zero real parts, respectively. The stable and unstable manifolds are well known to be unique and as smooth as the dynamical system, while center manifolds are non-unique and not all of them are guaranteed to be as smooth as the dynamical system (Guckenheimer and Holmes^{2} and Hirsch *et al.*^{3}). Near fixed points that only have eigenvalues with negative and zero real parts, center manifolds attract all trajectories, and hence the dynamics restricted to it provides a mathematically exact reduced-order model for the full system (Carr^{4} and Roberts^{5}).

Dissipative physical systems, however, generically have hyperbolic equilibria and hence admit no center manifolds. Instead, near their stable equilibria, such systems tend to have a set of fastest-decaying modes that die out quickly, leaving a set of slower decaying modes to govern the longer-term dynamics. Slow manifolds (i.e., SSMs constructed over such slower decaying modes) then replace center manifolds as targets for mathematically justified model reduction. Such slow SSMs were first targeted via Taylor expansions as nonlinear normal modes (NNMs) by Shaw and Pierre.^{6} These insightful calculations were then extended by the same authors to periodically and quasi-periodically forced mechanical systems to approximate forced mechanical response in a number of settings (see, e.g., the reviews by Kerschen *et al.*,^{7} Mikhlin and Avramov,^{8} Touzé *et al.*,^{9} and Mikhlin and Avramov^{10}).

Later mathematical analysis of general SSMs yielded precise existence, uniqueness, and smoothness results for these manifolds. Specifically, if the spectral subspace $E$ comprises either only decaying modes or only growing modes (i.e., $E$ is a like-mode spectral subspace) with no integer resonance relationships to the modes outside of $E$, then the slow SSM family $W(E)$ has a unique, primary member, $ W \u221e(E)$, that is as smooth as the full dynamical system (Cabré *et al.*^{11} and Haller and Ponsioen^{1}). The remaining secondary (or fractional) members of the SSM family have reduced but precisely known order of differentiability that depends on the ratio of linearized decay rates outside $E$ to those inside $E$ (Haller *et al.*^{12}). If we further disallow any integer resonance in the full linearized spectrum, then primary and fractional SSMs also exist when $E$ is of mixed-mode type, i.e., spanned by a combination of stable and unstable linear modes (Haller *et al.*^{12}).

All these results also hold for discrete autonomous dynamical systems, and hence SSM results also cover time-periodic continuous dynamical systems when applied to their Poincaré maps. Based on this fact, numerical implementations of time-periodic SSM calculations for periodically forced finite-element structures have appeared in Ponsioen *et al.*,^{13–15} Jain and Haller,^{16} and Vizzaccaro.^{17} An open source, equation-driven MATLAB toolbox (SSMTool) with a growing collection of worked problems is also available for mechanical systems with general nonlinearities (Jain *et al.*^{18}). Data-driven construction of time-periodic SSMs have been developed and applied to numerical and experimental data by Cenedese *et al.*,^{19,20} Kaszás *et al.*^{21} and Axås *et al.*,^{22} with open-source MATLAB implementations (SSMLearn and fastSSM) available from Cenedese *et al.*^{23}

Outside autonomous and time-periodic non-autonomous systems, the existence of like-mode SSMs has only been treated under small, time-quasiperiodic perturbations of autonomous dynamical systems (Haro and de la Llave,^{24} Haller and Ponsioen,^{1} Opreni *et al.*,^{25} and Thurner *et al.*,^{26}). In such systems, the role of hyperbolic fixed points as anchor points for SSMs is taken over by invariant tori whose dimension is equal to the number of rationally independent frequencies present in the forcing. An SSM in such a case perturbs from the direct product of an unperturbed invariant torus with an underlying spectral subspace $E$. Such SSMs have been shown to exist for like-mode spectral subspaces $E$ under small quasi-periodic perturbations, provided that the real part of the linearized spectrum within $E$ has no integer relationships with the real part of the spectrum outside $E$ (Haro and de la Llave,^{24} Haro *et al.*,^{27} and Haller and Ponsioen^{1}).

Related work by Fontich *et al.*^{28} covers the persistence and smoothness of primary SSMs emanating from an arbitrary attracting orbit of an autonomous dynamical system. While all non-autonomous systems become autonomous when their phase space is extended with the time variable, this theory does not apply under such an extension. The reason is that all attracting orbits lose hyperbolicity in the extended phase space due to the presence of the neutrally stable time direction.

In summary, while available SSM results have proven highly effective in equation-driven and data-driven reduced-order modeling of autonomous, time-periodic, and time-quasiperiodic systems, they offer no theoretical basis or computational scheme for physical systems with general aperiodic time-dependence. Yet the aperiodic setting is clearly of importance in a number of problems, including turbulent fluid–structure interactions, civil engineering structures subject to benchmark aperiodic forcing (such as those mimicking earthquakes) and control of robot motion.

In this paper, we extend available SSMs results to systems with general time dependence. While several powerful linearization results imply the existence of invariant manifolds for such non-autonomous dynamical systems, these results guarantee either no smoothness for SSM-type invariant manifolds (see, e.g., Palmer^{29}) or rely on conditions involving the Lyapunov spectra, Sacker–Sell spectra, or dichotomy spectra of some associated non-autonomous linear systems of ODEs (see, e.g., Yomdin,^{30} Pötzsche,^{31} and Cuong *et al.*^{32}). The latter types of conditions are intuitively clear but not readily verifiable, especially not in a data-driven setting. Here, our objective is to conclude the existence of time-aperiodic SSMs under directly computable conditions that also lead to explicitly computable SSM-reduced models in equation-driven and data-driven applications.

To this end, we consider two settings that arise frequently in practice: weak and additive non-autonomous time dependence and slowly varying (or adiabatic) time dependence. The first setting of weak non-autonomous external forcing is common in structural vibrations, wherein a structure’s steady-state response is of interest under various moderate loading scenarios. So far, related studies have been restricted to temporally periodic or quasiperiodic forcing (see, e.g., Ponsioen *et al.*,^{13–15} Li *et al.*,^{33,34} Jain and Haller,^{16} Vizzaccaro *et al.*,^{35} and Opreni *et al.*,^{25}), because the existence and exact form of a steady state and its associated SSM have been unknown for more general forcing profiles. The second setting of slow time dependence arises, for instance, in control applications wherein the intended motion of a structure is generally much slower than the characteristic time scales of its internal vibrations. In those applications, the lack of an adiabatic SSM theory has so far confined model-reduction studies to small-amplitude trajectories along which a single, autonomous SSM computed at a nearby fixed point was used for modeling purposes (Alora *et al.*^{36,37}).

In both of these non-autonomous settings, we use, modify, or extend prior invariant manifold results and techniques to conclude the existence of weakly aperiodic or adiabatic SSMs in the limit of small enough or slow enough time dependence, respectively. We then derive explicit recursive formulas for the arising non-autonomous SSMs and the aperiodic anchor trajectories to which they are attached. These formulas also cover and extend temporally periodic and aperiodic SSM computations to arbitrarily high order of accuracy. Using simple mechanical examples subjected to chaotic excitation, we illustrate that the new asymptotic formulas yield accurate reduced-order models even for larger and faster forcing.

## II. NON-AUTONOMOUS SSMs UNDER WEAK FORCING

### A. Setup

*eigenspace*$ E k$ of $A$ is then the linear span of the real and imaginary parts of the eigenvectors corresponding to the eigenvalue $ \lambda k$, i.e.,

*spectral subspace*$E$ is the direct sum of a selected group of $\u2113$ eigenspaces, i.e.,

*stable subspace*$ E s$ and the

*unstable subspace*$ E u$, defined as

### B. Existence and computation of a non-autonomous anchor trajectory

We first state general results on the fate of the $x=0$ fixed point under the non-autonomous forcing term $ f 1(x,t)$. Specifically, we give conditions under which this fixed point perturbs for moderately large $ | f 1 ( x , t ) |$ into a unique nearby hyperbolic trajectory $ x \u2217(t)$ of (1). This distinguished trajectory remains uniformly bounded for all times and has the same stability type as the $x=0$ fixed point of system (3).

#### (Existence of anchor trajectory for non-autonomous SSMs)

**(Existence of anchor trajectory for non-autonomous SSMs)**

See Appendix A 1.

#### (Computation of anchor trajectory for non-autonomous SSMs)

**(Computation of anchor trajectory for non-autonomous SSMs)**

See Appendix A 2.

#### (Applicability to temporally discontinuous forcing)

**(Applicability to temporally discontinuous forcing)**

The asymptotic approximation (15) only requires $ f 1$ and its $x$-derivatives at $x=0$ to be uniformly bounded in $t$, as we noted earlier. No derivatives of $ f 1(x,t)$ are required to exist with respect to $t$, and hence temporally discontinuous forcing is also covered by these formal expansions, as we will also see on a specific example with discontinuous chaotic forcing in Sec. IV A 1.

#### (Simplification for state-independent forcing)

**(Simplification for state-independent forcing)**

In applications to structural vibrations, the non-autonomous forcing term in system (1) often arises from external forcing that does not depend on $x$, i.e., $ f 1(x,t)\u2261 f 1(t)$. In that case, the second term in formula (15) is identically zero and the summands in the first term are well defined as long as $ f 0(x)$ and its derivatives are bounded at $x=0$.

### C. Existence and computation of non-autonomous SSMs

We first recall available SSM results for the autonomous system (3). Let $E$ be a spectral subspace, as defined in Eq. (8). Following the terminology of Haller *et al.*,^{12} we call $E$ a *like-mode spectral subspace*if the sign of $ Re \lambda j$ is the same for all $ \lambda j\u2208 spect ( A | E )$. Otherwise, we call $E$ a *mixed-mode spectral subspace.*

*externally nonresonant*if

^{1}and then the same condition will also be satisfied for

*all*choices of $ m k$ with $ \u2211 k = 1 n m k\u22652$ (Cabré

*et al.*

^{11}).

If $E$ is an externally nonresonant, like-mode subspace, then $E$ has a unique, smoothest nonlinear continuation in the form of a primary spectral submanifold (or *primary SSM*), denoted $ W \u221e ( E )\u2208 C \u221e$. This was deduced by Haller and Ponsioen^{1} in this context from the more abstract, general results of Cabré *et al.*^{11} The primary SSM $ W \u221e ( E )$ is an invariant manifold of system (3) that is tangent to $E$ at the origin and has the same dimension as $E$, as shown in Fig. 1 for $ f 1(x,t)\u22610$ in the extended phase space.

In general, an infinite family, $ W ( E )$, of SSMs exists for system (3) with the same properties, but $ W \u221e ( E )$ is the unique smoothest among them: the other manifolds in $ W ( E )$ are no more than $\sigma ( E )$-times differentiable (see Haller *et al.*^{12}). A computational algorithm for $ W \u221e ( E )$ was developed by Ponsioen *et al.*,^{13} then extended to general, finite-element-grade problems in second-order ODE form by Jain and Haller.^{16} The latest version of the latter algorithm with further extensions is available in the open-source MATLAB live script package SSMTool (see Jain *et al.*^{18}).

*fully nonresonant*, i.e.,

*secondary*(or

*fractional*)

*SSMs*, which have a finite order of differentiability consistent with their fractional-powered polynomial representations derived explicitly by Haller

*et al.*

^{12}

Note that the nonresonance conditions (17)–(18) are less restrictive than what is customary in the nonlinear vibration literature. Indeed, conditions (17)–(18) are only violated if both the real and the complex parts of the eigenvalues satisfy simultaneously exactly the same resonance relationship. Also note that 1:1 resonances involving eigenvalues with nonzero real parts do *not* violate (17)–(18). We will, however, only allow a 1:1 resonance within $E$ to guarantee the normal hyperbolicity of $E$ for the existence of $ W ( E )$. (More specifically, a 1:1 resonance between the spectrum of $A$ within $E$ and outside $E$ would render our upcoming assumption (19) to fail for any $\rho >1$.) If a 1:1 resonance arises in a given application, one can simply enlarge $E$ to include all resonant modes. This, in turn, removes any issue with the 1:1 resonance.

*slow SSMs*(also called

*pseudo-unstable manifolds*), which are continuations of $d$-dimensional, $\rho $-

*normally attracting*(like-mode or mixed-mode) spectral subspaces $E$ of linearized system (7), i.e., can be written as

Note that the third condition in (19) always holds for arbitrary $\rho >1$ if $ Re \lambda k$ and $ Re \lambda k + 1$ have different signs. If $ Re \lambda k$ and $ Re \lambda k + 1$ have the same sign, then the third condition in (19) holds only if the decay exponents of solutions outside $E$ are at least $\rho $-times stronger than those of solutions inside $E$. This $\rho $ will then determine the maximal smoothness that we can *a priori* guarantee for the SSM emanating from $E$ without further assumptions. Formal expansions for such SSMs will, however, indicate higher degrees of smoothness under appropriate nonresonance conditions.

The following theorem gives our main result on the existence of non-autonomous SSMs associated with a spectral subspace $E$ satisfying (19). We will use the term *locally invariant manifold* when referring to a manifold carrying trajectories that can only leave the manifold through its boundary (see, e.g., Fenichel^{38}).

#### (Existence of non-autonomous SSMs)

**(Existence of non-autonomous SSMs)**

*1*are satisfied in a ball $ B \delta \u2282 R n$ of radius $\delta >0$ around $x=0$. Assume further that $ f 0, f 1\u2208 C r( B \delta )$ for some $r>1$ and the uniform bound,

#### (Uniqueness and smoothness of non-autonomous SSMs)

**(Uniqueness and smoothness of non-autonomous SSMs)**

The manifolds discussed in Theorem 3 are generally non-unique. In principle, an argument involving the smoothness of the Lyapunov and dichotomy spectra under generic perturbation (see Son^{39}) could be invoked to deduce sharper smoothness results for these manifolds from linearization (see Yomdin^{30}). These results, however, require the identification of the full Lyapunov spectrum of the linear part of the non-autonomous linear differential equation (B1) in our Appendix B 1 in order to exclude resonances. This is generally challenging for equations and unrealistic for datasets. As an alternative, Theorem 4 below will provide unique Taylor expansions for $W(E,t)$ under explicitly verifiable nonresonance conditions. These expansions are varied on each persisting non-autonomous SSM up to the degree of smoothness of the SSM. We conjecture that, just as in the time-periodic and time-quasiperiodic case treated by Cabré *et al.*^{11} and Haro and de la Llave,^{24} a unique member of the $W(E,t)$ family of manifolds will be smoother than all the others. Our Taylor expansion will then approximate those unique, smoothest manifolds at orders higher than the spectral quotient $\Sigma (E)$ defined in Haller and Ponsioen.^{1} The remaining manifolds can be constructed via time-dependent versions of the fractional expansions identified in Haller *et al.*^{12}

#### (Non-autonomous SSMs without anchor trajectories)

**(Non-autonomous SSMs without anchor trajectories)**

For stable hyperbolic fixed points in the $ f 1(x,t)\u22610$ limit, an alternative to the proof in Appendix B 1 for Theorem 3 can also be given. This alternative uses the theory of non-compact normally hyperbolic invariant manifolds (see Eldering^{40}) coupled with the “wormhole” construct of Eldering *et al.*^{41} that enables the handling of inflowing-invariant normally attracting invariant manifolds. Following the steps of our proof in Appendix C 2 for the adiabatic case, this alternative proof yields results similar to those in Theorem 3 but does not rely on a persisting anchor trajectory $ x \u2217(t)$ near the origin. As a consequence, it can also capture non-autonomous SSMs in weakly damped physical systems for higher forcing levels at which $ x \u2217(t)$ is already destroyed. This is because the strength of hyperbolicity of the unforced fixed point at $x=0$ (measured by $ | Re \lambda 1 |\u226a1$) is generally much weaker than the strength of hyperbolicity of the unforced SSM, $ W \u221e ( E )$ $ ( measured by | Re \lambda k + 1 | | Re \lambda k | > 1 )$ in weakly damped systems. Such persisting SSMs without a stable hyperbolic anchor trajectory have been well-documented in equation- and data-driven studies of time-periodically forced systems. Those SSMs are signaled by overhangs near resonances in the forced response curves at higher forcing levels (see, e.g., Jain and Haller^{16} and Cenedese *et al.*^{19}).

As our examples will show, $W(E,t)$ will generally persist even under $ f 1(x,t)$ perturbations that are significantly larger than those allowed by the rather conservative assumptions of Theorem 1. In addition, $W(E,t)$ will also be smoother than $ C \rho $ under additional nonresonance conditions. To this end, we will next derive numerically implementable, recursive approximation formulas that are valid for $W(E,t)$ as long as it persists.

#### (Computation of non-autonomous SSMs and their reduced dynamics)

**(Computation of non-autonomous SSMs and their reduced dynamics)**

*3*exists in the coordinates $(u,v)$ for system (26) for all $\u03f5\u2208 [ 0 , \u03f5 \u2217 ]$. Assume further that $ W \u03f5(E,t)$ is $N$-times continuously differentiable for any fixed $t$ and the nonresonance conditions,

*,*:

- The SSM $ W \u03f5(E,t)$ admits a formal asymptotic expansion,The uniformly bounded $ h k p(t)$ coefficients in this expansion can be computed recursively from their initial conditions,$ W \u03f5 ( E , t ) = { ( u , v ) \u2208 U \u2282 R n : v = h \u03f5 ( u , t ) \u2211 | ( k , p ) | \u2265 1 N = \u2211 | ( k , p ) | \u2265 1 N h k p ( t ) u k \u03f5 p + o ( | u | q \u03f5 N \u2212 q )} .$via the formula,$ h 0 p ( t ) \u2261 0 , t \u2208 R , p \u2208 N ; h k 0 ( t ) \u2261 \u2212 A k \u2212 1 M k 0 ( h j 0 ) , | j | < | k | ; h 0 0 = 0 ,$where the functions $ M k p$ are defined as$ h k p ( t ) = \u222b \u2212 \u221e \u221e G k ( t \u2212 s ) M k p ( s , h j m ( s ) ) d s , | ( j , m ) | < | ( k , p ) | , t \u2208 R ,$$ M k p ( t , h j m ) = \u2202 | ( k , p ) | \u2202 u 1 k 1 \u2026 \u2202 u d k d \u2202 \u03f5 p [ f ^ v ( u , \u2211 | ( j , m ) | \u2265 1 | ( k , p ) | \u2212 1 h j m ( t ) u j \u03f5 m , \u03f5 ; t ; k , p ) \u2212 \u2211 | ( j , m ) | \u2265 1 | ( k , p ) | \u2212 1 \u03f5 p ( h 1 j m ( t ) j 1 u j u 1 \cdots h 1 j m ( t ) j d u j u d \vdots \ddots \vdots h n \u2212 d j m ( t ) j 1 u j u 1 \cdots h n \u2212 d j m ( t ) j d u j u d ) \xd7 f ^ u ( u , \u2211 | ( j , m ) | \u2265 1 | ( k , p ) | \u2212 1 h j m ( t ) u j \u03f5 m , \u03f5 ; t ; k , p ) ] | u = 0 , \u03f5 = 0 .$
- The reduced dynamics on $ W \u03f5 ( E , t )$ is obtained by restricting the $u$-component of system (26) to $ W \u03f5 ( E , t )$, which yieldsHere, $ Q u\u2208 C d \xd7 n$ is a matrix whose $j$th row is $ e ^ j/ ( e ^ j \u22c5 e j )$ for $j=1,\u2026d$, where $ e ^ j$ is the $j$th unit left eigenvector of $P$ corresponding to its unit right eigenvector $ e j$. Equivalently, $ Q u$ is composed of the first $d$ rows of $ P \u2212 1$.$ u \u02d9 = A u u + Q u [ f 0 ( x \u03f5 \u2217 ( t ) + P ( u h \u03f5 ( u , t ) ) ) + \u03f5 f ~ 1 ( x \u03f5 \u2217 ( t ) + P ( u h \u03f5 ( u , t ) ) , t ) + A x \u03f5 \u2217 ( t ) \u2212 x \u02d9 \u03f5 \u2217 ( t ) ] .$
- Let $ ( \xi , \eta ) T= P \u2212 1x$ denote coordinates aligned with the subspace $E$ and the direct sum of the remaining eigenspaces, respectively, emanating from the original $x=0$ fixed point of system (1) for $\u03f5=0$. In these coordinates, the reduced dynamics (38) becomes$ \xi \u02d9 = A u \xi + Q u [ f 0 ( x \u03f5 \u2217 ( t ) + P ( \xi \u2212 Q u x \u03f5 \u2217 ( t ) h \u03f5 ( \xi \u2212 Q u x \u03f5 \u2217 ( t ) , t ) ) ) + \u03f5 f ~ 1 ( x \u03f5 \u2217 ( t ) + P ( \xi \u2212 Q u x \u03f5 \u2217 ( t ) h \u03f5 ( \xi \u2212 Q u x \u03f5 \u2217 ( t ) , t ) ) , t ) ] .$

See Appendix B 2.

#### (Related results)

**(Related results)**

Pötzsche and Rassmussen^{42} derive Taylor approximations for a broad class of invariant manifolds emanating from fixed points of non-autonomous ODE. Due to the generality of that setting, the resulting formulas are less explicit than ours, assume global boundedness on the nonlinear terms, and also assume that the origin remains a fixed point for all times (and hence disallow additive external forcing). Nevertheless, under further assumptions and appropriate reformulations, those formulas should ultimately yield results equivalent to ours when applied to SSMs. Although it does not directly overlap with our results here, we also mention related work by Pötzsche and Rassmussen^{43} on the computation of classic and strong stable/unstable manifolds, as well as center-stable/unstable manifolds, for non-autonomous difference equations.

#### (Extension to discontinuous forcing)

**(Extension to discontinuous forcing)**

The derivatives in the definition of (37) can be evaluated using the same higher-order chain-rule formula of Constantine and Savits^{44} that we applied in the proof of Theorem 2. Using these formulas reveal that $ h k p(t)$ can be formally computed as long as $ f 1$ and its $x$-derivatives at $x=0$ remain uniformly bounded in $t$. Therefore, the formal expansion we have derived for the SSM $ W \u03f5 ( E , t )$ is valid for small enough $\u03f5$ even if $ f 1$ is discontinuous in the time $t$. We illustrate the continued accuracy of these expansions for discontinuous forcing on an example in Sec. IV A 1.

#### (Relation to the case of periodic or quasi-periodic forcing)

**(Relation to the case of periodic or quasi-periodic forcing)**

The recursively defined coefficient vectors $ M k p(s, h j m(s))$ defined in Eq. (37) are explicitly time dependent but otherwise satisfy the same formulas as the terms on the right-hand side of the invariance equations arising from autonomous SSM calculations. Consequently, the recursive formulas originally implemented in SSMTool by Jain *et al.*^{18} for autonomous SSM calculations apply here without modification. The difference is that those autonomous $ M k$ coefficients are used in solving linear algebraic systems of equations for $ h k$, as opposed to the non-autonomous $ M k p$ coefficients are used in linear ODEs for $ h k p(t)$.

#### (Accuracy of asymptotic SSM formulas)

**(Accuracy of asymptotic SSM formulas)**

The accuracy of the reduced-order dynamics (38) increases with increasing $ | ( k , p ) |$. There is no general convergence result under $N\u2192\u221e$ in the reduced dynamics (38), but such convergence is known if $f(x,t)$ is analytic and has periodic or quasi-periodic time dependence (Cabré *et al.*^{11} and Haro and de la Llave^{24}). For those types of time dependencies, the $ O ( \u03f5 )$ term in the reduced dynamics (38) is the same as that used for time-periodic and time-quasiperiodic SSMs in earlier publications (see, e.g., Breunung and Haller^{45} and Cenedese *et al.*^{20}).

^{45}Ponsioen

*et al.*,

^{15}Jain and Haller,

^{16}and Cenedese

*et al.*,

^{19}). In that case, the leading-order correction to the unforced reduced dynamics is simply the third line in (40), which is just the projection of the restriction of the forcing term to the autonomous SSM onto the spectral subspace $E$. If, in addition, the forcing term only depends on time (i.e., $ \u2202 x f ~ 1\u22610$), as is often the case in structural vibrations, then a leading-order approximation for small-amplitude motions on the SSM $ W \u03f5 ( E , t )$ becomes

Alternative parametrizations of $ W \u03f5 ( E , t )$ beyond the graph-style parametrization worked out here in detail are also possible (see Haller and Ponsioen,^{1} Haro *et al.*,^{27} and Vizzaccaro *et al.*^{35}). One option is the normal-form style parameterization which simultaneously brings the reduced dynamics on $ W \u03f5 ( E , t )$ to normal form. Normal form transformations can also be applied separately once $ W \u03f5 ( E , t )$ is located (see, e.g., Cenedese *et al.*^{19}).

## III. ADIABATIC SSMs

### A. Setup

*critical manifold*in the language of geometric singular perturbation theory (Fenichel

^{46}).

### B. Existence and computation of an adiabatic anchor trajectory

As in the non-autonomous case treated in Sec. II, we first discuss the continued existence of an anchor trajectory that acts as a continuation of the critical manifold $ L 0$.

#### (Existence and computation of anchor trajectory for adiabatic SSMs)

**(Existence and computation of anchor trajectory for adiabatic SSMs)**

See Appendix C 1.

### C. Existence and computation of an adiabatic SSM

*adiabatic SSM*anchored along the slow invariant manifold $ L \u03f5$ that is spanned by $ x \u03f5(\alpha )$ for small $\u03f5>0$. To this end, we assume that the first $k$ eigenvalues $ \lambda 1(\alpha ),\u2026, \lambda k(\alpha )$ in the list (48) have negative real parts and there is a nonzero spectral gap between $ Re \lambda k(\alpha )$ and $ Re \lambda k + 1(\alpha )$,

^{1}). As a consequence,

^{38}and Eldering

*et al.*

^{41}More specifically, $ M 0$ is a $\rho $-normally attracting, non-compact, inflowing-invariant manifold with a boundary, as shown in Fig. 2 for $\u03f5=0$ .

We have the following main results on the persistence of the manifold $ M 0$ in the form of an adiabatic SSM for small enough $\u03f5$.

#### (Existence of adiabatic SSMs)

**(Existence of adiabatic SSMs)**

Assume that $ x 0(\alpha )$ and $f(x,\alpha )$ have $r$ continuous derivatives that are uniformly bounded in $\alpha \u2208 R$ in a small, closed neighborhood of $ M 0$. Let $m=min ( r , \rho )$. Then, for $\u03f5>0$ small enough, there exists a persistent invariant manifold *(adiabatic SSM anchored along $ L \u03f5$), denoted* $ M \u03f5$, that is $ C m$ diffeomorphic to $ M 0$, has $m$ uniformly bounded derivatives and is $ O ( \epsilon )$ $ C 1$ -close to $ M 0$. Furthermore, $ L \u03f5\u2282 M \u03f5$ holds.

See Appendix B 2.

We show the geometry of the adiabatic SSM, $ M \u03f5$, in Fig. 2 for $\u03f5>0$.

#### (Applicability to mixed-mode adiabatic SSMs)

**(Applicability to mixed-mode adiabatic SSMs)**

For the purposes of constructing $ M \u03f5$, we had to assume in Eq. (53) that the spectrum of $A(\alpha )$ lies in the negative complex half plane, i.e., $E(\alpha )$ contains only stable directions. In the terminology of Haller *et al.*,^{12} this means that $W ( E ( \alpha ) )$ has to be a like-mode SSM for all values of $\alpha $ for Theorem 6 to apply. This enables us to select an inflowing-invariant, normally attracting manifold $ M ~ 0$ in our proof to which the wormhole construct in Proposition B1 of Eldering *et al.*^{41} is applicable. There is every indication that $ M \u03f5$ also exists under the weaker assumption (46), which only requires $A(\alpha )$ to have no eigenvalues on the imaginary axes, and hence allows $W ( E ( \alpha ) )$ to be a mixed-mode SSM that contains both stable and unstable directions. In this case, however, an $ M ~ 0$ with an inflowing or overflowing boundary cannot be chosen and hence Proposition B1 of Eldering *et al.*^{41} would need a technical extension that we will not pursue here, although it appears doable. We simply note that the asymptotic formulas we will derive for $ M \u03f5$ are formally valid under the weaker assumption (46) as well [i.e., for mixed-mode $W ( E ( \alpha ) )$], as long as $E(\alpha )$ is normally hyperbolic, i.e., attracts solutions at a rate that is uniformly stronger than the contraction rates inside $E(\alpha )$. We will evaluate and confirm these formulas in an example with a mixed-mode adiabatic SSMs in Sec. IV B 2.

As in the general non-autonomous case treated in Sec. II, the adiabatic $ M \u03f5$ will generally persist for larger values of $\u03f5$ and will also be smoother than $ C m$ under additional nonresonance conditions. In the following, we will provide a numerically implementable recursive scheme for computing $ M \u03f5$, assuming that it exists and is smooth enough.

#### (Computation of adiabatic SSMs)

**(Computation of adiabatic SSMs)**

*6*exists in the transformed system (56) for all $\u03f5\u2208 [ 0 , \u03f5 \u2217 ]$ . Assume further that $ M \u03f5$ is $N$-times continuously differentiable and the nonresonance conditions,

Then, for all $\u03f5\u2208 [ 0 , \u03f5 \u2217 ]$.

- The adiabatic SSM, $ M \u03f5$, admits a formal asymptotic expansionThe functions $ h k p(\alpha )$ are uniformly bounded in $\alpha \u2208 R$ for all $ ( k , p )$ and can be computed recursively from their initial conditions,$ M \u03f5= { ( u , v , \alpha ) \u2208 U \u2282 R n : v = h \u03f5 ( u , \alpha ) = \u2211 | ( k , p ) | \u2265 1 N h k p ( \alpha ) u k \u03f5 p + o ( | u | q \u03f5 N \u2212 q )}.$via the formula$ h 0 p ( \alpha ) \u2261 0 , p \u2265 0 ; h k 0 ( \alpha ) \u2261 h k ( \alpha ) , k \u2208 N d ; h k 0 ( \alpha ) \u2261 h k 0 ( \alpha ) \u2261 h k ( \alpha ) = 0 , | k | = 1 ; h k ( \u2212 1 ) ( \alpha ) := 0 ,$where$ h k p(\alpha )= A k \u2212 1(\alpha ) [ [ h k ( p \u2212 1 ) ] \u2032 ( \alpha ) \u2212 M k p ( \alpha , h j m ) ], | ( j , m ) |< | ( k , p ) |,$Specifically, for $\u03f5=0$, the coefficients in expansion (59) for $ M 0$ can be computed as$ A k ( \alpha ) = diag [ \lambda \u2113 ( \alpha ) \u2212 \u2211 j = 1 d k j \lambda j ( \alpha ) ] \u2113 = d + 1 n \u2208 C ( n \u2212 d ) \xd7 ( n \u2212 d ) , M k p ( \alpha , h j m ) = \u2202 | ( k , p ) | \u2202 u 1 k 1 \u2026 \u2202 u d k d \u2202 \u03f5 p [ f ^ v ( u , \u2211 | ( j , m ) | \u2265 1 | ( k , p ) | \u2212 1 h j m ( \alpha ) u j \u03f5 m , \u03f5 ; \alpha ; k , p ) ( h 1 j m ( \alpha ) j 1 u j u 1 \cdots h 1 j m ( \alpha ) j d u j u d \vdots \ddots \vdots h n \u2212 d j m ( \alpha ) j 1 u j u 1 \cdots h n \u2212 d j m ( \alpha ) j d u j u d ) \u2212 \u2211 | ( j , m ) | \u2265 1 | ( k , p ) | \u2212 1 \u03f5 p ( h 1 j m ( \alpha ) j 1 u j u 1 \cdots h 1 j m ( \alpha ) j d u j u d \vdots \ddots \vdots h n \u2212 d j m ( \alpha ) j 1 u j u 1 \cdots h n \u2212 d j m ( \alpha ) j d u j u d ) f ^ u ( u , \u2211 | ( j , m ) | \u2265 1 | ( k , p ) | \u2212 1 h j m ( \alpha ) u j \u03f5 m , \u03f5 ; \alpha ; k , p ) ] | u = 0 , \u03f5 = 0 .$$ h k(\alpha )\u2261 h k 0(\alpha )=\u2212 A k \u2212 1(\alpha ) M k(\alpha , h j 0), | j |< | k |.$

*The reduced dynamics on*$ M \u03f5$*is given by*$ u \u02d9 = Q u ( \alpha ) [ f ( x \u03f5 ( \alpha ) + P ( \alpha ) ( u h \u03f5 ( u , \alpha ) ) , \alpha ) \u2212 \u03f5 x \u03f5 \u2032 ( \alpha ) ] + \u03f5 Q u \u2032 ( \alpha ) P ( \alpha ) ( u h \u03f5 ( u , \alpha ) ) , \alpha \u02d9 = \u03f5 .$*Here, $ Q u ( \alpha )\u2208 C d \xd7 n$ is a matrix whose $j$th row is $ e ^ j(\alpha )/ ( e ^ j ( \alpha ) \u22c5 e j ( \alpha ) )$ for $j=1,\u2026d$, where $ e ^ j(\alpha )$ is the $j$th unit left eigenvector of $P(\alpha )$ corresponding to its unit right eigenvector $ e j(\alpha )$. Equivalently, $ Q u(\alpha )$ is composed of the first $d$ rows of $ P \u2212 1(\alpha )$*.

See Appendix C 3.

### D. Special case: Adiabatic SSM under additive slow forcing

We now assume that the slow time dependence in system (43) arises from slow (but generally not small) additive forcing. This is a reasonable assumption, for instance, in the setting of external control forces acting on a highly damped structure, whose unforced decay rate to an equilibrium is much faster than the rate at which the forcing changes in time. Another relevant setting is very slow (quasi-static) forcing in structural dynamics.

These simplifications imply that the matrix $A( x 0(\alpha ))$, the column matrix $P( x 0(\alpha ))$ of the right eigenvectors of $A( x 0(\alpha ))$ and the row matrix $ Q u( x 0(\alpha ))$ of the first $d$, appropriately normalized left eigenvectors of $A( x 0(\alpha ))$ can now be determined solely from the unforced part $ f 0(x)$ of the right-hand side of system (69) along any parametrized path $ x 0(\alpha )$ defined implicitly by Eq. (71). Similarly, an inspection of the formulas (66) defining the slow SSM, $ M 0$, reveals that the graph $ h 0(u, x 0(\alpha ))$ of $ M 0$ can be computed along the path $ x 0(\alpha )$ solely based on the knowledge of $ f 0(x)$; only the path itself depends on the specific forcing term $ f 1(\alpha )$.

*a priori*compute the parametric family of frozen-time reduced models,

*a priori*along different spatial locations $q$ of the configuration space under the application of static forces $ f 1 v(\alpha )$ that keep the mechanical system at equilibrium ( $v=0$) at those $q(\alpha )$ locations. This is a significant simplification over the general setting of Eq. (74), because it only requires the construction of the model family (74) over a codimension- $n$ subset $ U q= { ( q , v ) \u2208 R n : v = 0}$ of the full phase space $ R n$ of system (69).

This simplified reduced model construction is expected to provide further improvement over current autonomous-SSM-based control strategies that use only a single member of the model family (74), computed at fixed point $ x 0 \u2217$ of the unforced system (see Alora *et al.*^{36,37}). This single model is then used along the full path $ x 0(\alpha )$ to approximate the leading-order model (73). Clearly, this approximation will generally only be valid when the instantaneous equilibrium path $ x 0(\alpha )$ remains close to the unperturbed fixed point $ x 0 \u2217$. We will explore the application of the ideas described in this section to robot control in other upcoming publications.

## IV. EXAMPLES

We formulated our results in Secs. II and III for fully non-dimensionalized systems for which smallness or slowness can simply be imposed by selecting a single perturbation parameter $\u03f5$ small enough. In specific physical examples, such a non-dimensionalization can be done in multiple ways and is often cumbersome to carry out for large systems. One nevertheless needs to assess whether the external forcing is *de facto* smaller or slower than the magnitude or the speed of variation of internal forces along trajectories. While such a consideration has been largely ignored in applications of perturbation results to physical problems, here we will propose and apply heuristic physical measures that help assessing the magnitude and the slowness of the perturbation.

*forcing weakness*$ r w$ and

*forcing speed*$ r s$ as

In a practical setting, we deem a specific external forcing weak if $ r w\u226a1$ or slow if $ r s\u226a1$. This is the range in which the perturbation-based invariant manifold techniques used in proving our relevant theorems can reasonably expected to apply. We have noted, however, that the SSMs obtained from these techniques are normally hyperbolic and hence persist under further increases in $ r w$ and $ r s$. Indeed, we will show that our formal asymptotic SSM expansions tend to remain valid for $ r w>1$ and and $ r s>1$. In some cases, we also find these formulas to have predictive power for $ r w, r s\u226b1$.

All MATLAB scripts used in producing the results for the examples below can be publicly accessed at Kaundinya and Haller.^{47}

### A. Example 1: Chaotically forced cart systems

#### 1. Weak forcing

We consider the two-degree-of-freedom (2 DOF) mechanical system from Haller and Ponsioen,^{1} placed now on a moving cart subject to external chaotic shaking, as shown in Fig. 3(a). In this example, the fixed point of the unforced system is asymptotically stable. As a result, under external forcing, the SSM attached to a unique nearby anchor trajectory will be of like-mode type.

We will consider two different cases of forcing by scaling the magnitude of the chaotic signal $g(t)$ in two different ways. In the first case, we will have $max | F ( t ) |=0.06( N)$, which gives the non-dimensional forcing weakness $ r w=0.16$ from the first formula in Eq. (80). This forcing scheme can, thus, be considered weak albeit not very small. Our subsequent choice of $max | F ( t ) |=3( N)$ gives $ r w=7.8$, which is definitely outside the range of small forcing.

Utilizing that the forcing signal $g(t)$ vanishes for $t<0$, we evaluate the integrals in (82) via the trapezoidal integration scheme of MATLAB. For sufficient accuracy, we use forcing values evaluated at $ 10 4$ equally spaced times to ensure accurate numerical integration. The terms $ x \nu (t)$ in the expansion for $ x \u2217(t)$ are then spline-interpolated in time to obtain an expression for $ x \u2217(t)$ for arbitrary $t\u2208 [ 0 , 500 ]$.

^{47}

To assess the performance of the truncated SSM-reduced model (84) for $N=5$, we first plot the normalized mean trajectory error computed from ten different initial conditions for $max | |F | |=0.06 ( N )$ in Fig. 4(a). We compute the dependence of this error on $max | |F | |$ in Fig. 4(b). As expected, the errors grow with $max | |F | |$ but remains an order of magnitude less than $max | |F | |$.

Probing larger forcing amplitudes, we still obtain accurate reduced-order models up to $max | |F | |=3 ( N )$. As an illustration, Fig. 5 uses the center-of-mass coordinate $ x c$ to show our 5th order asymptotic approximation using Theorem 1 for the anchor trajectory of the chaotic SSM (blue), a simulation from the full model for a trajectory in the 2D SSM converging to the anchor point (black), and a prediction from the 2D SSM-reduced model obtained from Theorem 4 for the same trajectory (red).

In Fig. 6 (Multimedia available online), we also plot snapshots of the chaotic 2D SSM, its anchor trajectory, a trajectory from the full-order model, and the prediction for this trajectory from the SSM-reduced order model at three different times. Note how both the full trajectory and the predicted model trajectory synchronize with the anchor trajectory of the SSM (see Fig. 6). Finally, we would like to illustrate that the improper integrals in our asymptotic formulas for anchor trajectories and SSMs indeed converge and give correct results even for discontinuous forcing, as noted in Remark 6. To this end, we make the forcing discontinuous at ten random points in time, as shown in Fig. 7(b). The same formulas for the anchor trajectory, its attached SSM, and the SSM-reduced order model remain formally applicable and continue to give accurate approximations even for the relatively large forcing amplitude of $max | |F | |=3 ( N )$, as seen in Fig. 7.

First, Fig. 8(a) shows for $ \gamma f=0 ( N / m 3 )$ (i.e., for the linear limit of the right-most spring), the $ O ( 1 )$ and $ O ( 11 )$ asymptotic approximations to the anchor trajectory of the chaotic SSM using the results of Theorem 1. Already in this case, the $ O ( 1 )$ approximation is noticeably improved by the $ O ( 11 )$ approximation, but the $ O ( 1 )$ approximation is also close to the actual generalized steady state.

In contrast, Fig. 8(b) shows a case in which much smaller forcing is applied but the right-most spring is nonlinear with $ \gamma f=0.5( N / m 3)$. In this case, more significant error develops between the actual anchor trajectory and its $ O ( 1 )$, while the $ O ( 11 )$ approximation remains indistinguishably close to the actual anchor trajectory. This example, therefore, underlines the need for higher approximations to the anchor point in the presence of stronger nonlinearities. As the SSM is attached to the anchor trajectory, the accuracy of an SSM-reduced model also depends critically on this refined approximation.

#### 2. Slow forcing

For the same cart system shown in Fig. 3(a), we select the physical parameters $ m 1= m 2=1( kg)$, $ M f=2( kg)$, $k= k f=1( N / m)$, and $c= c f=0.3( Ns / m)$. We now slow down the previously applied chaotic forcing by replacing the external forcing vector in Eq. (81) with $ F ( \alpha )= ( 0 , 0 , N s g ( \alpha ) ) T$. We select this forcing to be of large amplitude by letting $ N s=10$ and consider the phase range $\alpha \u2208[0,6]$. In our analysis, we will start with the minimal forcing speed $\u03f5= \alpha \u02d9=0.001$, for which obtain the non-dimensional forcing speed $ r s=0.72$ from the second formula in (80). This qualifies as moderately slow external forcing relative to the speed of variation of internal forces along unperturbed trajectories. In contrast, for the maximal forcing speed $\u03f5= \alpha \u02d9=0.008$, we will consider, we obtain $ r s=7.26$, which can no longer be considered slow. Again, of particular interest will be how our asymptotic formulas for adiabatic SSMs and their anchor points perform in the latter case, which is clearly faster than what is normally considered slowly varying in perturbation theory.

Once $ x 0(\alpha )$ is available numerically, we evaluate the recursive formula (51) up to order $N=3$ to obtain

We set $N=3$ and solve for all the $\alpha $ dependent coefficients in the expansion of the 2D slow chaotic SSM by solving a recursive linear algebraic equations (64) order by order. We find that all 28 coefficients that describe the SSM up to cubic order are nonzero. An important numerical step used in this task is to perform differentiation with respect to $\alpha $. For this purpose, we first need to re-orient the unit eigenvectors originally returned by MATLAB for specific values of $\alpha $ to make these eigenvectors smooth functions of $\alpha $. We then perform a central finite differencing which includes four adjacent points in the $\alpha $ direction. We finally employ the Savitzky–Golay filtering function of MATLAB to smoothen the results obtained from finite differencing an already finite-differenced signal.

To illustrate the ultimate accuracy of these numerical procedures, Table I shows the normalized mean trajectory error for predictions from two different SSM-reduced models. We see that under increasing $\u03f5$ (i.e., faster forcing), the errors still remain an order-of-magnitude smaller than $\u03f5$.

. | $O(\u03f5,j=1, |m |=3)$ . | $O ( \u03f5 3 , j + | m | = 3 )$ . |
---|---|---|

$\u03f5=0.001$ (r_{s} = 0.72) | 9.8 × 10^{−5} | 7.2 × 10^{−6} |

$\u03f5=0.008$(r_{s} = 7.26) | 6.4 × 10^{−3} | 7.8 × 10^{−4} |

$\u03f5=0.010$ (r_{s} = 8.66) | 1.1 × 10^{−2} | 2.9 × 10^{−3} |

. | $O(\u03f5,j=1, |m |=3)$ . | $O ( \u03f5 3 , j + | m | = 3 )$ . |
---|---|---|

$\u03f5=0.001$ (r_{s} = 0.72) | 9.8 × 10^{−5} | 7.2 × 10^{−6} |

$\u03f5=0.008$(r_{s} = 7.26) | 6.4 × 10^{−3} | 7.8 × 10^{−4} |

$\u03f5=0.010$ (r_{s} = 8.66) | 1.1 × 10^{−2} | 2.9 × 10^{−3} |

As we did for weak forcing, we track the center of mass coordinate $ x c$ to test the accuracy of our asymptotic formulas up to order $N=3$ for $\u03f5=0.008$ in Fig. 9. Note how this cubic-order approximation of the SSM-reduced model already tracks the full solution closely. Finally, as in the weakly forced case, we also plot in Fig. 10 snapshots of the evolution of the slow anchor trajectory, the attached chaotic SSM, a simulated full-order trajectory, and its prediction from the same initial condition based on the SSM-reduced dynamics (Multimedia available online). These snapshots further confirm the accuracy of our analytic approximation formulas, now in several coordinate directions.

### B. Example 2: Chaotically forced bumpy rail

#### 1. Weak forcing

The physical setup of our second example, a particle moving on a shaken rail with a bump in the middle, is shown in Fig. 11(a). A major difference from our previous mechanical example is that the origin of the unforced system is now an unstable, saddle-focus-type fixed point. Therefore, under weak chaotic forcing, a nearby chaotic SSM of mixed-mode type is expected to arise attached to a saddle-focus-type chaotic anchor trajectory.

with for the coordinate vector $ x= ( x c , x ) T$. The parameter $\beta $ controls the height of the bumpy rail and $a$ the distance to the wells. The parameter $c$ models linear damping of the mass $m$ due to the rail. The velocity-dependent nonlinearity $ f( x, x \u02d9)$ is the result of expanding the exact equations of motion with respect to the height of the rail and keeping only the leading-order nonlinearities.

We set $m=1( kg)$, $ M f=4( kg)$, $ k f=1( N / m)$, $ c f=c=0.3( Ns / m)$, $a=0.3( m)$, $\beta = 1 5 a 3$, and $g=9.8( m / s 2)$. The forcing signal $g(t)$ will be the same as in the cart example. We scale the chaotic forcing signal in a way that the non-dimensional forcing weakness measure defined in formula (80) gives $ r w=44.6$ when evaluated on the phase space region containing the trajectories used in our analysis. Therefore, the external force here is large relative to the internal forces of the system.

Without listing the details here, we also perform a similar calculation for the two attracting anchor trajectories created by the chaotic forcing from the two asymptotically stable fixed points of the unperturbed system that lie at the two bottom points of the bumpy rail. These two stable anchor trajectories will represent two chaotic attractors for the forced system. Trajectories initialized near the unstable anchor trajectory on its attached 2D mixed-mode SSM will converge to one of these chaotic attractors. This is indeed observable in Fig. 11(b), in which we launch several trajectories near the unstable anchor trajectory (red) within the 2D mixed-mode SSM approximated up to cubic order. These trajectories then converge to one of the two predicted attracting chaotic anchor trajectories (blue) which we have computed up to cubic order of accuracy from our formulas.

The instability of the anchor trajectory perturbing from the origin makes it challenging to verify the accuracy of our local expansion (64) for the 2D mixed-mode SSM attached to this trajectory. We can nevertheless verify the local accuracy for these formulas by tracking nearby trajectories launched from the predicted 2D SSM and confirm in Fig. 12 (Multimedia available online) that those trajectories evolve together with the predicted SSM until ejected from a neighborhood of the unstable anchor trajectory.

#### 2. Slow forcing

In the same bumpy rail problem but now with height $\beta = 1 10 3$, we rescale the forcing in time to make it slowly varying. The slow forcing vector is given by $ F ( \alpha )= ( N s g ( \alpha ) 0 )$ with $ N s=3$ and the range of interest is $\alpha \u2208[0,6].$ We will perform simulations for the forcing speed $\u03f5= \alpha \u02d9=0.01$, which yields the non-dimensional forcing speed measure $ r s=1.22$, signaling moderately fast forcing.

As in the weakly forced version of this example, we plot the predictions from the 2D SSM-reduced order model and compare it to simulations of the full system in Fig. 13(a) (Multimedia available online). We again see close agreement of the reduced dynamics with the true solution in the vicinity of the slow, unstable anchor trajectory, prior to the convergence of the trajectories to the two stable, chaotic anchor trajectories arising from the two stable unforced equilibria along the rail. Finally, as in our first example, we show snapshots of the projected 2D slow, mixed-mode SSM along with some trajectories of its restricted dynamics. These again stay close locally to trajectories of the full system that are launched from the same initial conditions, thereby confirming the expectations put forward in Remark 9.

## V. CONCLUSIONS

In this paper, we have derived existence results for spectral submanifolds (SSMs) in non-autonomous dynamical systems that are either weakly non-autonomous or slowly varying (adiabatic). While our exact proofs for these SSMs results only cover weak enough or slow enough time dependence, the invariant manifolds we obtain are structurally stable and hence are persistent away from these limits. Under further nonresonance conditions, the manifolds and their reduced dynamics also admit asymptotic expansions up to any finite order for which we present recursive formulas. These formulas suggest that one of the non-autonomous SSMs is as smooth as the original dynamical system, just as primary SSMs are known to be in the autonomous case reviewed in Sec. I. MATLAB implementations of our recursive formulas for the examples treated here are available from Kaundinya and Haller.^{47}

In our results, weakly forced SSMs are guaranteed to exist under uniformly bounded forcing. This restriction is expected because unbounded perturbations to an autonomous limit of a dynamical system will generally wipe out any structure identified in that limit. For adiabatic SSMs, the uniformity of the forcing magnitude is replaced by the uniformity of the strength of hyperbolicity along the anchor trajectory of the SSM in the limit of frozen time.

In the weakly forced case, the SSMs emanate from unique, uniformly bounded hyperbolic anchor trajectories. Our expansions for these anchor trajectories are of independent interest as they provide formal approximations of the *generalized steady state* response of any nonlinear system subject to moderate but otherwise arbitrary time-dependent forcing. In finite-element simulations under aperiodic forcing, such steady states have been assumed to exist but their computation tends to involve lengthy simulations of randomly chosen initial conditions until they have reached a (somewhat vaguely defined) statistical steady state.

Even with today’s computational power, such simulations are still prohibitively long due to the small damping of typical structural materials, the high degrees of freedom of the models used, and the costly evaluations of nonlinearities arising from complex geometries and multi-physics. With the explicit recursive formulas we obtained here, one can now directly approximate these generalized steady states without lengthy simulations.

For these asymptotic formulas to converge, the forcing does not actually have to be uniformly bounded: it may grow temporally unbounded in a compact neighborhood of the origin as long as the improper integrals in the statements of Theorems 2 and 4 converge. To the best of our knowledge, such explicit, readily computable asymptotic expansions for time-dependent steady states have been unavailable in the literature even for time-periodic or time-quasiperiodic forcing, let alone time-aperiodic forcing. Weak or slow time-periodic and time-quasi-periodic forcing is also covered by these formulas as special cases.

Our existence results for SSMs assume temporal smoothness for the dynamical system and hence do not strictly cover the case of random forcing with continuous paths. Our approximation formulas, however, are formally applicable (i.e., the improper integrals in them converge) for much broader forcing types, including mildly exponentially growing or temporally discontinuous forcing. Therefore, the results in this paper provide formal nonlinear model-reduction formulas that can be implemented as approximations for all forcing types arising in practice.

The simple examples considered here already illustrate the ability of our asymptotic formulas to work under various forcing types, including chaotic and discontinuous forcing. Our weakly non-autonomous SSM theory is expected to be helpful in model reduction for structural vibration problems in which so far time-periodic and time-quasiperiodic SSMs have been constructed. An application of our adiabatic SSM theory is already under way in the model-predictive control of soft robots, wherein the target trajectories of the robot generally have a much slower time scale than the highly damped robot itself.

Going from the simple, illustrative examples treated here to finite-element problems with temporally aperiodic forcing will require no further theoretical development. Indeed, the existence and smoothness results derived in this paper are valid in arbitrarily high (finite) dimensions. However, a computational reformulation will be required to make these results directly applicable to second-order mechanical systems without the need to convert them to first-order systems with diagonalized linear parts. The development of such a reformulation using the approach of Jain and Haller^{16} is currently under way.

## ACKNOWLEDGMENTS

We are grateful to Christian Pötzsche for very helpful and detailed explanations and references on the continuity of the dichotomy spectrum for non-autonomous differential equations. We are also grateful to Martin Rasmussen and Peter Kloeden for helpful references on non-autonomous invariant manifolds.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**George Haller:** Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Methodology (equal); Supervision (lead); Writing – original draft (lead); Writing – review & editing (lead). **Roshan S. Kaundinya:** Methodology (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The codes required to generate these numerical results are openly available at https://github.com/haller-group/Aperiodic-SSMs, Ref. 47.

### APPENDIX A: PROOF OF THEOREMS 1 AND 2

#### 1. Proof of Theorem 1: Existence and regularity of a hyperbolic anchor trajectory

^{29}reformulated from Hartman,

^{48}for a general non-autonomous ODE of the form

*exponential dichotomy*, i.e., there exists a fundamental matrix solution $\Phi (t)$ of (A2), constants $K,\kappa >0$ and a projection $P$ (with $ P 2=P)$ such that

^{29}Furthermore, $ x \u2217(t;p)$ is continuous in the parameter $p$. [If $f$ is smooth in its arguments, then the continuity of $ x \u2217(t;p)$ in $p$ can be strengthened to smooth dependence on $p$ using more general results for the persistence of non-compact, normally hyperbolic invariant manifolds; see Eldering.

^{40}Such a strengthened version, however, would not be specific enough about the norm of $ | x \u2217 ( t ; p ) |$ to the extent given in (A5).] Note that $ x \u2217(t;p)$ takes over the role of $x=0$ as a unique, uniformly bounded solution under nonzero $f(x,t,p)$. The result of Palmer

^{29}quoted above makes no assumption on $f(x,t,p)$ containing only nonlinear terms, and hence $f(x,t,p)$ is allowed to be an arbitrary perturbation to the linear ODE (A2), as long as it is uniformly bounded and uniformly Lipschitz.

The uniform boundedness conditions (A4)–(A5) will not be satisfied in realistic applications. Nevertheless, one can still use the above results in such applications using smooth cutoff functions (see, e.g., Fenichel^{38}) that make $f(x,t,p)$ vanish for all $t\u2208 B \delta $ outside a small ball $ B \delta \u2282 R n$. The latter classic tool from invariant manifold theory is, however, only applicable here if the predicted hyperbolic trajectory $ x \u2217(t;p)$ lies inside $ B \delta $, where the cutoff right-hand side and the original right-hand side of (A1) still coincide.

Statement (ii) of Theorem 1 follows from the smooth dependence of $ x \u2217(t;p)$ on the parameter vector $p$, as shown by Palmer.^{29}

#### 2. Proof of Theorem 2: Approximation of the anchor trajectory *x**(*t*)

^{44}for higher-order derivatives of composite functions. To recall the general form of this formula, we first consider a general composite function $H: R p\u2192 R$, defined as

$ | \mu |< | \nu |$,

$ | \mu |= | \nu |$ and $ \mu 1< \nu 1$, or

$ | \mu |= | \nu |$ and $ \mu 1= \nu 1,\u2026 \mu k= \nu k, \mu k + 1< \nu k + 1$ for some $1\u2264k<p$.

^{44}prove the following multi-variate version of the Faá di Bruno formula:

### APPENDIX B: PROOFS OF THEOREMS 3 AND 4

#### 1. Proof of Theorem 3: Existence of non-autonomous spectral submanifolds

^{49}do apply. Their main condition stated in our context is

$ | \u2202 x f 0 ( x \u2217 ( t ) ) |$ is uniformly bounded,

$ | \u2202 x f 1 ( x , t ) |$ and $ | \u2202 x 2 f 1 ( x , t ) |$ are uniformly bounded in the $ B \delta $ neighborhood of $x=0$,

$ | g ( y , t ) |$ is uniformly bounded in a neighborhood of $y=0$, and

$ \u2202 y | g ( y , t ) |$ is uniformly bounded in a neighborhood of $y=0$.

First, note that condition (a) is satisfied because $ x \u2217(t)$ is uniformly bounded. Second, since $ x \u2217(t)$ stays in the $ B \delta $ ball under the conditions (11), we obtain that $ | \u2202 x f 0 ( x \u2217 ( t ) + y ) \u2212 \u2202 x f 0 ( x \u2217 ( t ) ) |$ in uniformly bounded in a compact neighborhood of $y=0.$ Therefore, for condition (B7) to hold [and hence to ensure that (c)–(d) to hold], it remains to require that $ | \u2202 x f 1 ( x \u2217 ( t ) + y ) \u2212 \u2202 x f 1 ( x \u2217 ( t ) ) |$ remain uniformly bounded in a neighborhood of $y=0$. By the mean value theorem, this holds if $ | \u2202 x 2 f 1 ( x , t ) |$ remains uniformly bounded in $ B \delta $, which is just condition (b) above. Therefore, in addition to our assumptions in (11), we need to assume condition (20) for (a)–(d) to hold. In summary, under the conditions (11) and (20), the local invariant manifold results of Kloeden and Rassmussen^{49} are applicable near the $Y=0$ fixed point of system (B4)–(B5).

^{49}apply to general non-autonomous systems of differential equations of the form (B4). These results in turn build on classic result of Sacker and Sell

^{50}that establish the existence of a dichotomy spectrum $\Sigma $ for $ A(t)$ that consists up to $n+1$ disjoint closed intervals of the form

^{49}for each $j=1,\u2026,m\u22121$, there exist two non-autonomous invariant manifolds $ W j \xb1(t)$ for system (B4) such that

$ W j \xb1(t)$ contain the $Y=0$ fixed point of system (B4).

- In a neighborhood $U\u2282 R n + 1,$ the manifolds $ W j \xb1(t)$ can described as graphs with the help of $ C 0$ functions $ w j \xb1 : U \xd7 R \u2192 R n + 1$ such that$ W j \xb1 ( t ) = { Y = Z + w j \xb1 ( Z , t ) \u2208 U : Z \u2208 range ( P \xb1 j ( t ) ) , w j \xb1 ( Z , t ) \u2208 range ( P \u2213 j ( t ) )} .$
$ w j \xb1(Z,t)$ is uniformly $o ( | Z | )$, i.e., $ lim u \u2192 0 \Vert w j \xb1 ( Z , t ) \Vert \Vert Z \Vert =0$, uniformly in $t$.

- Ifholds for a positive integer $ m j +$, then $ W j +(t)$ is of class $ C m j +$. Similarly, if$ m j + \kappa j +< \kappa j \u2212$holds for a positive integer $ m j \u2212$ , then $ W j \u2212(t)$ is of class $ C m j \u2212$.$ \kappa j +< m j \u2212 \kappa j \u2212$
If $ \kappa j +<0$ then for all $\gamma > \kappa j +$, we have $ sup t \u2265 0 \Vert Y ( t , t 0 , Y 0 ) \Vert e \u2212 \gamma t<\u221e$ for all $ Y 0\u2208 W j +( t 0)\u2229U$ in a small enough neighborhood $U$ of $Y=0$. Similarly, if $ \kappa j \u2212>0$ then for all $\gamma < \kappa j \u2212$, we have $ sup t \u2265 0 \Vert Y ( t , t 0 , Y 0 ) \Vert e \u2212 \gamma t<\u221e$ for all $ Y 0\u2208 W j +( t 0)\u2229U$ in a small enough neighborhood $U$ of $Y=0$.

^{49}assume no hyperbolicity for the $Y=0$ fixed point of system (B4), which makes them applicable to the $Y=0$ fixed point of the extended system (B4). Specifically, the dichotomy spectrum $\Sigma $ of the matrix $ A(t)$ defined in Eq. (B5) is discrete and given by

^{49}point out, the definition (B13) yields a hierarchy of invariant manifolds, which in turn implies the hierarchy of invariant manifolds shown in Table II for the original system (1) arising from this construct for the extended system.

$ W 1 +(t)$ | ⊂ | $ W 2 +(t)$ | ⊂ | ⋅⋅ ⋅ | ⊂ | $ W n \u2212 1 +(t)$ | ⊂ | $ R n + 1$ |

$\u222a$ | $\u222a$ | $\u222a$ | $\u222a$ | |||||

$ W 2 , 1(t)$ | ⊂ | ⋅⋅ ⋅ | ⊂ | $ W n \u2212 1 , 1(t)$ | ⊂ | $ W 1 \u2212(t)$ | ||

$\u222a$ | $\u222a$ | |||||||

$\ddots $ | ⋮ | ⋮ | ||||||

$ W n \u2212 1 , n \u2212 2(t)$ | ⊂ | $ W n \u2212 2 \u2212(t)$ | ||||||

$\u222a$ | ||||||||

$ W n \u2212 1 \u2212(t)$ |

$ W 1 +(t)$ | ⊂ | $ W 2 +(t)$ | ⊂ | ⋅⋅ ⋅ | ⊂ | $ W n \u2212 1 +(t)$ | ⊂ | $ R n + 1$ |

$\u222a$ | $\u222a$ | $\u222a$ | $\u222a$ | |||||

$ W 2 , 1(t)$ | ⊂ | ⋅⋅ ⋅ | ⊂ | $ W n \u2212 1 , 1(t)$ | ⊂ | $ W 1 \u2212(t)$ | ||

$\u222a$ | $\u222a$ | |||||||

$\ddots $ | ⋮ | ⋮ | ||||||

$ W n \u2212 1 , n \u2212 2(t)$ | ⊂ | $ W n \u2212 2 \u2212(t)$ | ||||||

$\u222a$ | ||||||||

$ W n \u2212 1 \u2212(t)$ |

The $\u03f5=$constant slices of $ W i +(t)$, $ W j \u2212(t)$, and $ W i , j(t)$ then provide similar invariant manifolds $ W i +(t;\u03f5)$, $ W j \u2212(t;\u03f5)$, and $ W i , j(t;\u03f5)$ of one less dimension along the $ x \u2217(t;\u03f5)$ trajectory of the original system (1) for small enough with $\u03f5>0$ under the assumed scaling (B3). The degree of smoothness of $ W i +(t;\u03f5)$, $ W j \u2212(t;\u03f5)$, and $ W i , j(t;\u03f5)$ in $x$ and $\u03f5$ coincides with the general degree of smoothness that can be inferred from the above construct for their extended counterparts, $ W j +(t)$, $ W j \u2212(t)$, and $ W i , j(t)$. Specifically, the following non-autonomous invariant manifolds can be inferred from Table II for system (1):

$ W i +(t;\u03f5)$: the time-dependent smooth continuations of the fastest decaying $i$ modes for any $1\u2264i\u2264n\u22121$. These manifolds satisfy (B11) for arbitrary large $m$ and hence are as smooth as the original dynamical system. They are generally $ C 0$ in $\u03f5$ which can be improved to $ C r$ for pseudo-stable manifolds, i.e., for all $ W i +(t;\u03f5)$ with $ \mu i>0$.

$ W j \u2212(t;\u03f5)$: the time-dependent smooth continuations of the slowest decaying (or even growing) $n\u22121\u2212j$ modes. If $ \mu j \u2212 1/ \mu j> m j \u2212$ for some positive integer $ m j \u2212$, then these manifolds satisfy (B12) with $ m j \u2212 1$ and hence are of smoothness class $ C m j \u2212$. For all $j$ with $ \mu j<0$, $ W j \u2212(t;\u03f5)$ is also $ C m j \u2212$ smooth in $\u03f5$. Any manifold $ W j \u2212(t;\u03f5)$ with $ \mu j>0$ (fast unstable manifolds) can only be concluded to be $ C 0$ in $\u03f5$ from this argument.

$ W i , j(t;\u03f5)$; the continuation of any (like-mode or mixed-mode) spectral subspace, that is, $ r i , j$-normally hyperbolic. In that case, $ W i , j(t;\u03f5)$ is also of smoothness class $ C r i , j$ in $\u03f5$.

#### 2. Proof of Theorem 4: Approximation of the SSM, *W* (*E*, *t*), anchored at $ x \u03f5 \u2217(t)$

*a. Invariance equation*

^{1}). For $\u03f5>0$ small, under the conditions of statement (i) of Theorem 3, this primary SSM persists in the form of an (generally non-unique) invariant manifold

^{1}).

Also note that for the smooth persistence of $ W 0(E)$ as $ W \u03f5 ( E , t )$, a possible 1:1 resonance between the an eigenvalue inside $E$ and another one outside $E$ has to be excluded in order to secure the normal hyperbolicity of $ W 0(E)$. This is reflected by the lowering of the lower index in the summation in Eq. (33) from $2$ to $1$ relative to what one requires for the existence of $ W 0 ( E )$.

*b. Structure of the invariance equation*

*c. Solution of the invariance equation*

Next, we note that the nonresonance assumption (18) implies that the coefficient matrix $ A k$ of the homogeneous part of Eq. (B29) is nonsingular. We additionally require now this homogenous part to admit a hyperbolic fixed point, which leads to the stronger nonresonance condition (32) of Theorem 3. This is the same non-resonance condition that arises in the work of Haro and de la Llave^{24} for quasiperiodic forcing, which is a subset of the general forcing class we are considering here.