and hence is not exactly of the form (A20). To address this issue, we also
observe that

$\displaystyle\frac{d^{\nu}}{d\epsilon^{\nu}}H\left(0\right)$
$\displaystyle=\frac{d^{\nu}}{d\epsilon^{\nu}}H_{0}\left(0\right)+\nu\frac{d^{\nu-1}}{d\epsilon^{\nu-1}}H_{1}\left(0\right),$
$\displaystyle H_{0}\left(\epsilon\right)$
$\displaystyle=f_{0q}\left(\sum_{j=1}^{\nu-1}\epsilon^{j}\xi_{j}^{1}(\tau),\ldots,\sum_{j=1}^{\nu-1}\epsilon^{j}\xi_{j}^{n}(\tau)\right),$
$\displaystyle H_{1}\left(\epsilon\right)$
$\displaystyle=\tilde{f}_{1q}\left(\sum_{j=1}^{\nu-1}\epsilon^{j}\xi_{j}^{1}(\tau),\ldots,\sum_{j=1}^{\nu-1}\epsilon^{j}\xi_{j}^{n}(\tau),\tau\right),$

and hence $H_{0}\left(\epsilon\right)$ and $H_{1}\left(\epsilon\right)$ are
individually of the form (A20). Applied to these two functions, the formula
(A21) simplifies to

$\displaystyle\frac{d^{\nu}}{d\epsilon^{\nu}}H_{0}\left(0\right)$
$\displaystyle=\left.D_{\epsilon}^{\nu}f_{0q}\left(\sum_{j=1}^{\nu-1}\epsilon^{j}\xi_{j}^{1}(\tau),\ldots,\sum_{j=1}^{\nu-1}\epsilon^{j}\xi_{j}^{n}(\tau)\right)\right|_{\epsilon=0}$
$\displaystyle=\sum_{1\leq\left|\boldsymbol{\gamma}\right|\leq\nu}D_{x}^{\boldsymbol{\gamma}}f_{q}\left(0\right)\sum_{s=1}^{\nu}\sum_{p_{s}\left(\nu,\boldsymbol{\gamma}\right)}\nu!\prod_{j=1}^{s}\frac{\prod_{i=1}^{n}\left[\left(\ell_{j}\right)!\xi_{\ell_{j}}^{i}(\tau)\right]^{k_{ji}}}{\left(\mathbf{k}_{j}\right)!\left[\left(\ell_{j}\right)!\right]^{\left|\mathbf{k}_{j}\right|}}$
$\displaystyle=\sum_{1\leq\left|\boldsymbol{\gamma}\right|\leq\nu}D_{x}^{\boldsymbol{\gamma}}f_{0q}\left(0\right)\sum_{s=1}^{\nu}\sum_{p_{s}\left(\nu,\boldsymbol{\gamma}\right)}\nu!\prod_{j=1}^{s}\frac{\prod_{i=1}^{n}\left[\xi_{\ell_{j}}^{i}(\tau)\right]^{k_{ji}}}{\prod_{i=1}^{n}k_{ji}!},\quad
q=1,\ldots,n,$
$\displaystyle\frac{d^{\nu-1}}{d\epsilon^{\nu-1}}H_{1}\left(0\right)$
$\displaystyle=\left.D_{\epsilon}^{\nu-1}\tilde{f}_{1q}\left(\sum_{j=1}^{\nu-1}\epsilon^{j}\xi_{j}^{1}(\tau),\ldots,\sum_{j=1}^{\nu-1}\epsilon^{j}\xi_{j}^{n}(\tau),\tau\right)\right|_{\epsilon=0}$
$\displaystyle=\sum_{1\leq\left|\boldsymbol{\gamma}\right|\leq\nu-1}D_{x}^{\mathbf{\boldsymbol{\gamma}}}\tilde{f}_{1q}\left(0,\tau\right)\sum_{s=1}^{\nu-1}\sum_{p_{s}\left(\nu-1,\boldsymbol{\gamma}\right)}\left(\nu-1\right)!\prod_{j=1}^{s}\frac{\prod_{i=1}^{n}\left[\left(\ell_{j}\right)!\xi_{\ell_{j}}^{i}(\tau)\right]^{k_{ji}}}{\left(\mathbf{k}_{j}\right)!\left[\left(\ell_{j}\right)!\right]^{\left|\mathbf{k}_{j}\right|}}.$

Substitution of these expressions into formula (A19) then gives the final
recursive formulas

$\displaystyle\xi_{\nu}(t)$
$\displaystyle=\sum_{1\leq\left|\boldsymbol{\gamma}\right|\leq\nu}\,\sum_{s=1}^{\nu}\sum_{p_{s}\left(\nu,\boldsymbol{\gamma}\right)}\int_{-\infty}^{t}e^{A\left(t-\tau\right)}\left[\frac{\partial^{\left|\boldsymbol{\gamma}\right|}f_{0}\left(0\right)}{\partial
x_{1}^{\gamma_{1}}\cdots\partial
x_{n}^{\gamma_{n}}}\prod_{j=1}^{s}\frac{\prod_{i=1}^{n}\left[\xi_{\ell_{j}}^{i}(\tau)\right]^{k_{ji}}}{\prod_{i=1}^{n}k_{ji}!}\right]d\tau$
(A23)
$\displaystyle+\sum_{1\leq\left|\boldsymbol{\gamma}\right|\leq\nu-1}\,\sum_{s=1}^{\nu-1}\sum_{p_{s}\left(\nu-1,\boldsymbol{\gamma}\right)}\int_{-\infty}^{t}e^{A\left(t-\tau\right)}\left[\frac{\partial^{\left|\boldsymbol{\gamma}\right|}\tilde{f}_{1}\left(0,\tau\right)}{\partial
x_{1}^{\gamma_{1}}\cdots\partial
x_{n}^{\gamma_{n}}}\prod_{j=1}^{s}\frac{\prod_{i=1}^{n}\left[\xi_{\ell_{j}}^{i}(\tau)\right]^{k_{ji}}}{\prod_{i=1}^{n}k_{ji}!}\right]d\tau,\quad\nu\geq
1,$

where $k_{ji}$ is the $i^{th}$ component of the integer vector
$\mathbf{k}_{j}\in\mathbb{N}^{n}-\left\\{\mathbf{0}\right\\}$ appearing in the
index set

$p_{s}\left(\nu,\boldsymbol{\gamma}\right)=\left\\{\left(\mathbf{k}_{1},\ldots,\mathbf{k}_{s},\ell_{1},\ldots,\ell_{s}\right):\mathbf{k}_{i}\in\mathbb{N}^{n}-\left\\{\mathbf{0}\right\\},\ell_{i}\in\mathbb{N},0<\ell_{1}<\cdots<\ell_{s},\sum_{i=1}^{s}\mathbf{k}_{i}=\boldsymbol{\gamma},\sum_{i=1}^{s}\left|\mathbf{k}_{i}\right|\ell_{i}=\nu\right\\}.$

We have also used the notational convection
$\frac{\partial^{\left|\boldsymbol{\gamma}\right|}}{\partial
x_{1}^{\gamma_{1}}\cdots\partial x_{n}^{\gamma_{n}}}=I$ for
$\boldsymbol{\gamma}=0$.

Substitution of (A23) into the expansion (A14) gives

$\displaystyle x^{*}(t)$ $\displaystyle=\sum_{\nu\geq 1}x_{\nu}(t),$ (A24)
$\displaystyle x_{\nu}(t)=$
$\displaystyle\sum_{1\leq\left|\boldsymbol{\gamma}\right|\leq\nu}\,\sum_{s=1}^{\nu}\sum_{p_{s}\left(\nu,\boldsymbol{\gamma}\right)}\int_{-\infty}^{t}e^{A\left(t-\tau\right)}\left[\frac{\partial^{\left|\boldsymbol{\gamma}\right|}f_{0}\left(0\right)}{\partial
x_{1}^{\gamma_{1}}\cdots\partial
x_{n}^{\gamma_{n}}}\frac{\prod_{i=1}^{n}\left[x_{\ell_{j}}^{i}(\tau)\right]^{k_{ji}}}{\prod_{i=1}^{n}k_{ji}!}\right]d\tau$
$\displaystyle+\sum_{1\leq\left|\boldsymbol{\gamma}\right|\leq\nu-1}\,\sum_{s=1}^{\nu-1}\sum_{p_{s}\left(\nu-1,\boldsymbol{\gamma}\right)}\int_{-\infty}^{t}e^{A\left(t-\tau\right)}\left[\frac{\partial^{\left|\boldsymbol{\gamma}\right|}f_{1}\left(0,\tau\right)}{\partial
x_{1}^{\gamma_{1}}\cdots\partial
x_{n}^{\gamma_{n}}}\prod_{j=1}^{s}\frac{\prod_{i=1}^{n}\left[x_{\ell_{j}}^{i}(\tau)\right]^{k_{ji}}}{\prod_{i=1}^{n}k_{ji}!}\right]d\tau,\quad\nu\geq
1.$

These results cover anchor trajectories of like-mode SSMs of stable hyperbolic
fixed points but do not cover anchor trajectories for mixed-mode SSMs of such
fixed points or anchor trajectories of like-mode SSMs of unstable hyperbolic
fixed points.

We now extend formula (A23) to mixed-mode SSMs by weakening the stability
assumption (A18) back to our general hyperbolicity assumption (6). In this
case, the matrix $A$ has an exponential dichotomy, as described by the
inequalities (10). The dichotomy exponent $\kappa>0$ can be selected as in
(A9). We use the matrix $T$ defined in (12) in the coordinate transformation

$x=Ty,\quad
y=\left(y^{s},y^{u}\right)^{\mathrm{T}}\in\mathbb{R}^{s}\times\mathbb{R}^{u},$

which brings system system (1) to the form

$\displaystyle\dot{y}^{s}$ $\displaystyle=A^{s}y^{s}+f^{s}(y,t),$
$\displaystyle\dot{y}^{u}$ $\displaystyle=A^{u}y^{u}+f^{u}(y,t),$

with $\left(f^{s}(y,t),f^{u}(y,t)\right)^{\mathrm{T}}=T^{-1}f(Ty,t)$ and

$\mathrm{Re}\left[\mathrm{spect}\left(A^{s}\right)\right]\subset\left(-\infty,0\right),\quad\mathrm{Re}\left[\mathrm{spect}\left(A^{u}\right)\right]\subset\left(0,\infty\right).$
(A25)

In these new coordinates, with the notation

$\left(\begin{array}[]{c}\hat{\xi}_{\nu}^{s}(t)\\\
\hat{\xi}_{\nu}^{u}(t)\end{array}\right)=T^{-1}\xi_{\nu}(t),$

the system of ODEs (A16) becomes

$\displaystyle\dot{\hat{\xi}}_{\nu}^{s}(t)$
$\displaystyle=A^{s}\hat{\xi}_{\nu}^{s}(t)+\frac{1}{\nu!}\left.D_{\epsilon}^{\nu}f^{s}\left(T^{-1}x_{\epsilon}(t),t\right)\right|_{\epsilon=0},$
$\displaystyle\dot{\hat{\xi}}_{\nu}^{u}(t)$
$\displaystyle=A^{u}\hat{\xi}_{\nu}^{u}(t)+\frac{1}{\nu!}\left.D_{\epsilon}^{\nu}f^{u}\left(T^{-1}x_{\epsilon}(t),t\right)\right|_{\epsilon=0},$
(A26)

whose solutions can be written as

$\displaystyle\hat{\xi}_{\nu}^{s}(t)$
$\displaystyle=e^{A^{s}\left(t-t_{0}^{s}\right)}\hat{\xi}_{\nu}^{s}(t_{0}^{s})+\frac{1}{\nu!}\int_{t_{0}^{s}}^{t}e^{A^{s}\left(t-\tau\right)}\left.D_{\epsilon}^{\nu}f^{s}\left(T^{-1}x_{\epsilon}(\tau),\tau\right)\right|_{\epsilon=0}d\tau,$
$\displaystyle\hat{\xi}_{\nu}^{u}(t)$
$\displaystyle=e^{A^{u}\left(t-t_{0}^{u}\right)}\hat{\xi}_{\nu}^{u}(t_{0}^{u})+\frac{1}{\nu!}\int_{t_{0}^{u}}^{t}e^{A^{u}\left(t-\tau\right)}\left.D_{\epsilon}^{\nu}f^{u}\left(T^{-1}x_{\epsilon}(\tau),\tau\right)\right|_{\epsilon=0}d\tau.$
(A27)

Based on the spectral properties of $A^{s}$ and $A^{u}$ listed in (A25) and
the uniform boundedness of $\hat{\xi}_{\nu}^{s,u}(t_{0}^{s,u})$, we can take
the limits $t_{0}^{s}\to-\infty$ and $t_{0}^{u}\to+\infty$ to obtain

$\displaystyle\hat{\xi}_{\nu}^{s}(t)$
$\displaystyle=\frac{1}{\nu!}\int_{-\infty}^{t}e^{A^{s}\left(t-\tau\right)}\left.D_{\epsilon}^{\nu}f^{s}\left(T^{-1}x_{\epsilon}(\tau),\tau\right)\right|_{\epsilon=0}d\tau,$
$\displaystyle\hat{\xi}_{\nu}^{u}(t)$
$\displaystyle=-\frac{1}{\nu!}\int_{t}^{\infty}e^{A^{u}\left(t-\tau\right)}\left.D_{\epsilon}^{\nu}f^{u}\left(T^{-1}x_{\epsilon}(\tau),\tau\right)\right|_{\epsilon=0}d\tau.$
(A28)

Therefore, introducing the Green’s function (13), we can write the solution
$\xi_{\nu}(t)$ in the $y$ coordinates as

$\xi_{\nu}(t)=\frac{1}{\nu!}\int_{-\infty}^{\infty}G(t-\tau)\left.D_{\epsilon}^{\nu}f\left(x_{\epsilon}(\tau),\tau\right)\right|_{\epsilon=0}d\tau,\quad\nu\geq
1.$

This then leads to the final formula (15) for a general hyperbolic trajectory
$x^{*}(t)$, in analogy with formula (A24) for an asymptotically stable
hyperbolic trajectory.

##  Appendix B Proofs of Theorems 3 and 4

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

Under the conditions of Theorem 1, we can shift coordinates by letting

$y=x-x^{*}(t),$

which transforms system (1) to the form

$\displaystyle\dot{y}$
$\displaystyle=\left[A+\partial_{x}f(x^{*}(t),t)\right]y+g(y,t),$ (B1)

where

$\displaystyle g(y,t)$
$\displaystyle=f(x^{*}(t)+y)-f(x^{*}(t))-\partial_{x}f(x^{*}(t),t)y=\mathcal{O}\left(\left|y\right|^{2}\right).$
(B2)

Next we introduce a perturbation parameter $0\leq\epsilon<\delta$ and the
scalings

$x=\epsilon\xi,\quad y=\epsilon\eta=\epsilon\xi-\epsilon\xi^{*}(t),\quad
f_{0}(x)=\epsilon^{2}\tilde{f}_{0}\left(\xi;\epsilon\right),\quad
f_{1}(x,t)=\epsilon\tilde{f}_{1}(x,t).$ (B3)

These scalings imply

$\displaystyle f_{0}(x^{*}(t)+y)$
$\displaystyle=\epsilon^{2}\tilde{f}_{0}(\xi^{*}(t)+\eta;\epsilon),\quad
D_{x}f_{0}(x^{*}(t))=\epsilon D_{\xi}\tilde{f}_{0}(\xi^{*}(t);\epsilon),$
$\displaystyle\partial_{x}f_{1}(x^{*}(t),t)$
$\displaystyle=\epsilon\partial_{x}\tilde{f}_{1}(x^{*}(t),t),\quad
g(y,t)=\epsilon^{2}\tilde{g}(\eta,t;\epsilon).$

With these expressions, we can rewrite eq. (B1) as

$\displaystyle\dot{\eta}$ $\displaystyle=\left[A+\epsilon
D_{\xi}\tilde{f}_{0}(\xi^{*}(t);\epsilon)+\epsilon\partial_{x}\tilde{f}_{1}(\epsilon\xi^{*}(t),t)\right]\eta+\epsilon\tilde{g}(\eta,t;\epsilon),\quad\tilde{g}(\eta,t)=\mathcal{O}\left(\left|\eta\right|^{2}\right),$
$\displaystyle\dot{\epsilon}$ $\displaystyle=0.$

We can further rewrite these equations as

$\displaystyle\dot{Y}$
$\displaystyle=\mathcal{A}(t)Y+\mathcal{G}(Y,t),\qquad\mathcal{G}(Y,t)=\mathcal{O}\left(\left|Y\right|^{2}\right),\qquad
Y=\left(\begin{array}[]{c}\eta\\\
\epsilon\end{array}\right)\in\mathbb{R}^{n+1},$ (B6)
$\displaystyle\mathcal{A}(t)$ $\displaystyle=\left(\begin{array}[]{cc}A&0\\\
0&0\end{array}\right),\qquad\mathcal{G}(Y,t)=\left(\begin{array}[]{c}\epsilon\partial_{x}\tilde{f}_{0}(\xi^{*}(t);\epsilon)\eta+\epsilon\partial_{x}\tilde{f}_{1}(\epsilon\xi^{*}(t),t)\eta+\epsilon\tilde{g}(\eta,t)\\\
0\end{array}\right).$ (B11)

While linearization results do not apply to eq. (B6) due to the non-
hyperbolicity of the origin, the invariant manifold results of Kloeden and
Rassmussen [29] do apply. Their main condition stated in our context is

$\mathcal{G}(0,t)=0,\qquad\lim_{Y\to
0}\sup_{t\in\mathbb{R}}\left|D_{Y}\mathcal{G}(Y,t)\right|=0,$ (B12)

of which the first one is already satisfied and the second one requires the
local uniform boundedness of the first derivatives of $\mathcal{G}(Y,t)$ in a
neighborhood of $Y=0$. By inspection of formula (B11), we see that the second
condition in (B12) is satisfied if the following four conditions are
fulfilled:

(a) $\left|\partial_{x}f_{0}(x^{*}(t))\right|$ is uniformly bounded

(b) $\left|\partial_{x}f_{1}(x,t)\right|$ and
$\left|\partial_{x}^{2}f_{1}(x,t)\right|$ are uniformly bounded in the
$B_{\delta}$ neighborhood of $x=0$

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

(d) $\partial_{y}\left|g(y,t)\right|$ is uniformly bounded in a neighborhood
of $y=0$.

Given the definition of $g(y,t)$ in (B2), conditions (c)-(d) are satisfied if
$\left|\partial_{x}f(x^{*}(t),t)\right|$ and

$\left|\partial_{x}f(x^{*}(t)+y,t)-\partial_{x}f(x^{*}(t),t)\right|$ (B13)

are uniformly bounded in $B_{\delta}$.

First, note that condition (a) is satisfied because $x^{*}(t)$ is uniformly
bounded. Second, since $x^{*}(t)$ stays in the $B_{\delta}$ ball under the
conditions (11), we obtain that
$\left|\partial_{x}f_{0}(x^{*}(t)+y)-\partial_{x}f_{0}(x^{*}(t))\right|$ in
uniformly bounded in a compact neighborhood of $y=0.$ Therefore, for condition
(B13) to hold (and hence to ensure that (c)-(d) to hold), it remains to
require that
$\left|\partial_{x}f_{1}(x^{*}(t)+y)-\partial_{x}f_{1}(x^{*}(t))\right|$
remain uniformly bounded in a neighborhood of $y=0$. By the mean value
theorem, this holds if $\left|\partial_{x}^{2}f_{1}(x,t)\right|$ 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 [29] are applicable
near the $Y=0$ fixed point of system (B6)-(B11).

Specifically, under conditions (11) and (20), the results in Section 6.3 of
Kloeden and Rassmussen [29] apply to general non-autonomous systems of
differential equations of the form (B6). These results in turn build on
classic result of Sacker and Sell [43] that establish the existence of a
dichotomy spectrum $\Sigma$ for $\mathcal{A}(t)$ that consists up to $n+1$
disjoint closed intervals of the form

$\Sigma=\left[a_{1},b_{1}\right]\cap\ldots\cap\left[a_{m},b_{m}\right],\quad
m\leq m+1.$ (B14)

(By the definition of the dichotomy spectrum, the linear system of ODE
$\dot{Y}=$ $\left(\mathcal{A}(t)-\lambda I\right)Y$ admits no exponential
dichotomy for any $\lambda\in\Sigma$.) Assuming $m\geq 2$, one can therefore
select constants $\kappa_{j}^{+},\kappa_{-j}^{-}\in\mathbb{R}$ for
$j=1,\ldots,m-1$ from the gaps among the closed spectral subintervals in (B14)
as follows:

$b_{j}<\kappa_{j}^{+}<\kappa_{j}^{-}<a_{j+1},\quad j=1,\ldots,m-1,$

such that for appropriate constants $K>0$ and projection maps
$P_{\pm}^{j}(t_{0})\in\mathbb{R}^{(n+1)\times\left(n+1\right)}$ with
$P_{\pm}^{j}(t_{0})P_{\pm}^{j}(t_{0})=P_{\pm}^{j}(t_{0})$, the normalized
fundamental matrix solution $\Phi\left(t,t_{0}\right)$ of
$\dot{Y}=\mathcal{A}(t)Y$satisfies

$\displaystyle\left\|\Phi\left(t,t_{0}\right)P_{-}^{j}(t_{0})\right\|$
$\displaystyle\leq Ke^{\kappa_{j}^{+}(t-t_{0})},\quad t\geq t_{0},$
$\displaystyle\left\|\Phi\left(t,t_{0}\right)P_{+}^{j}(t_{0})\right\|$
$\displaystyle\leq Ke^{\kappa_{j}^{-}(t-t_{0})},\quad t\leq t_{0},$ (B15)

for all $j=1,\ldots,m-1$.

Then, by Theorem 6.10 and Remark 6.6 of Kloeden and Rassmussen [29], for each
$j=1,\ldots,m-1$, there exist two non-autonomous invariant manifolds
$W_{j}^{\pm}(t)$ for system (B6) such that:

(i)

$\mathcal{W}_{j}^{\pm}(t)$ contain the $Y=0$ fixed point of system (B6).

(ii)

In a neighborhood $U\subset\mathbb{R}^{n+1},$ the manifolds
$\mathcal{W}_{j}^{\pm}(t)$ can described as graphs with the help of $C^{0}$
functions $w_{j}^{\pm}\colon U\times\mathbb{\mathbb{R}}\to\mathbb{R}^{n+1}$
such that

$\mathcal{W}_{j}^{\pm}(t)=\left\\{Y=Z+w_{j}^{\pm}(Z,t)\in
U:\,\,Z\in\mathrm{range}\left(P_{\pm}^{j}(t)\right),\,\,\,\,w_{j}^{\pm}(Z,t)\in\mathrm{range}\left(P_{\mp}^{j}(t)\right)\right\\}.$
(B16)

(iii)

$w_{j}^{\pm}(Z,t)$ is uniformly $o$$\left(\left|Z\right|\right)$, i.e.,
$\lim_{u\to 0}\frac{\left\|w_{j}^{\pm}(Z,t)\right\|}{\left\|Z\right\|}=0$,
uniformly in $t$.

(iv)

If

$m_{j}^{+}\kappa_{j}^{+}<\kappa_{j}^{-}$ (B17)

holds for a positive integer $m_{j}^{+}$, then $\mathcal{W}_{j}^{+}(t)$ is of
class $C^{m_{j}^{+}}.$ Similarly, if

$\kappa_{j}^{+}<m_{j}^{-}\kappa_{j}^{-}$ (B18)

holds for a positive integer $m_{j}^{-}$ , then $\mathcal{W}_{j}^{-}(t)$ is of
class $C^{m_{j}^{-}}.$

(v)

If $\kappa_{j}^{+}<0$ then for all $\gamma>\kappa_{j}^{+}$, we have
$\sup_{t\geq 0}\left\|Y(t,t_{0},Y_{0})\right\|e^{-\gamma t}<\infty$ for all
$Y_{0}\in\mathcal{W}_{j}^{+}(t_{0})\cap U$ in a small enough neighborhood $U$
of $Y=0$. Similarly, if $\kappa_{j}^{-}>0$ then for all
$\gamma<\kappa_{j}^{-}$, we have $\sup_{t\geq
0}\left\|Y(t,t_{0},Y_{0})\right\|e^{-\gamma t}<\infty$ for all
$Y_{0}\in\mathcal{W}_{j}^{+}(t_{0})\cap U$ in a small enough neighborhood $U$
of $Y=0$.

Based on these results, further non-autonomous invariant manifolds can be
obtained by letting

$\mathcal{W}_{i,j}(t)=\mathcal{W}_{i}^{+}(t)\cap\mathcal{W}_{j}^{-}(t),\quad
1\leq j<i\leq m-1.$ (B19)

Note that Theorem 6.10 and Remark 6.6 of Kloeden and Rassmussen [29] assume no
hyperbolicity for the $Y=0$ fixed point of system (B6), which makes them
applicable to the $Y=0$ fixed point of the extended system (B6). Specifically,
the dichotomy spectrum $\Sigma$ of the matrix $\mathcal{A}(t)$ defined in eq.
(B11) is discrete and given by

$\Sigma=\left\\{\mu_{1},\mu_{2},\ldots,\mu_{c}=0,\ldots,\mu_{m+1}\right\\}=\mathrm{Re}\left[\mathrm{spect}\left(A\right)\right]\cup\left\\{0\right\\}.$

As Kloeden and Rassmussen [29] point out, the definition (B19) yields a
hierarchy of invariant manifolds, which in turn implies the hierarchy of
invariant manifolds shown in Table 1 for the original system (1) arising from
this construct for the extended system.

$\mathcal{W}_{1}^{+}(t)$ | $\subset$ | $\mathcal{W}_{2}^{+}(t)$ | $\subset$ | $\cdots$ | $\subset$ | $\mathcal{W}_{n-1}^{+}(t)$ | $\subset$ | $\mathbb{R}^{n+1}$
---|---|---|---|---|---|---|---|---
|  | $\cup$ |  | $\cup$ |  | $\cup$ |  | $\cup$
|  | $\mathcal{W}_{2,1}(t)$ | $\subset$ | $\cdots$ | $\subset$ | $\mathcal{W}_{n-1,1}(t)$ | $\subset$ | $\mathcal{W}_{1}^{-}(t)$
|  |  |  |  |  | $\cup$ |  | $\cup$
|  |  |  | $\ddots$ |  | $\vdots$ |  | $\vdots$
|  |  |  |  |  | $\mathcal{W}_{n-1,n-2}(t)$ | $\subset$ | $\mathcal{W}_{n-2}^{-}(t)$
|  |  |  |  |  |  |  | $\cup$
|  |  |  |  |  |  |  | $\mathcal{W}_{n-1}^{-}(t)$
Table 2: Hierarchy of invariant manifolds for the extended system (B6).

For any index $j\leq c$, all invariant manifolds $\mathcal{W}_{i,j}^{+}(t)$
defined in formula (B19) are graphs over a spectral subspace that contains the
$\mu_{c}=0$ center direction and hence they are as smooth in $\epsilon$ as
they are in the other variables along that spectral subspace. This smoothness
property is in turn shared by all invariant manifolds of the form

$\displaystyle\mathcal{W}_{i,c}(t)$
$\displaystyle=\mathcal{W}_{i}^{+}(t)\cap\mathcal{W}_{c}^{-}(t),\quad 1\leq
c<i\leq m-1.$

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

(A) $W_{i}^{+}(t;\epsilon)$: the time-dependent smooth continuations of the
fastest decaying $i$ modes for any $1\leq i\leq n-1$. These manifolds satisfy
(B17) for arbitrary large $m$ and hence are as smooth as the original
dynamical system. They are generally $C^{0}$ in $\epsilon$ which can be
improved to $C^{r}$ for pseudo-stable manifolds, i.e., for all
$W_{i}^{+}(t;\epsilon)$ with $\mu_{i}>0$.

(B) $W_{j}^{-}(t;\epsilon)$: the time-dependent smooth continuations of the
slowest decaying (or even growing) $n-1-j$ modes. If
$\mu_{j-1}/\mu_{j}>m_{j}^{-}$ for some positive integer $m_{j}^{-}$, then
these manifolds satisfy (B18) with $m_{j}^{-1}$ and hence are of smoothness
class $C^{m_{j}^{-}}$. For all $j$ with $\mu_{j}<0$, $W_{j}^{-}(t;\epsilon)$
is also $C^{m_{j}^{-}}$ smooth in $\epsilon$. Any manifold
$W_{j}^{-}(t;\epsilon)$ with $\mu_{j}>0$ (fast unstable manifolds) can only be
concluded to be $C^{0}$ in $\epsilon$ from this argument.

(C) $W_{i,j}(t;\epsilon)$; 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;\epsilon)$ is also of smoothness class $C^{r_{i,j}}$ in $\epsilon$.

###  B.2 Proof of Theorem 4: Approximation of the SSM, $W\left(E,t\right)$,
anchored at $x_{\epsilon}^{*}(t)$

####  B.2.1 Invariance equation

Next, we derive an approximation for the non-autonomous SSMs anchored to the
hyperbolic trajectory $x^{*}(t)$. We perform a linear change of coordinates

$\left(\begin{array}[]{c}u\\\
v\end{array}\right)=P^{-1}\left(x-x_{\epsilon}^{*}(t)\right),$

where $P=[e_{1},\ldots,e_{n}]\in\mathbb{C}^{n}$ contains the complex
eigenvectors corresponding to the ordered eigenvalues (5) of $A$, and
$\left(u,v\right)\in\mathbb{C}^{d}\times\mathbb{C}^{n-d}$. Rewriting the
scaled version (A13) of system (1) in these complex coordinates, we obtain

$\displaystyle\left(\begin{array}[]{c}\dot{u}\\\ \dot{v}\end{array}\right)$
$\displaystyle=P^{-1}\left(\dot{x}-\dot{x}_{\epsilon}^{*}(t)\right)=P^{-1}\left[Ax+f(x,t)-\dot{x}_{\epsilon}^{*}(t)\right]$
(B22) $\displaystyle=P^{-1}AP\left(\begin{array}[]{c}u\\\
v\end{array}\right)+P^{-1}Ax_{\epsilon}^{*}(t)+P^{-1}\left[f\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
v\end{array}\right),t\right)-\dot{x}_{\epsilon}^{*}(t)\right]$ (B27)
$\displaystyle=\left(\begin{array}[]{cc}A^{u}&0\\\
0&A^{v}\end{array}\right)\left(\begin{array}[]{c}u\\\
v\end{array}\right)+\hat{f}(u,v,\epsilon;t),$ (B32)

where

$\displaystyle\hat{f}(u,v,\epsilon;t)$
$\displaystyle=P^{-1}\left[f_{0}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
v\end{array}\right)\right)+\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
v\end{array}\right),t\right)+Ax_{\epsilon}^{*}(t)-\dot{x}_{\epsilon}^{*}(t)\right],$
(B37)

and

$\displaystyle A^{u}=\left(\begin{array}[]{ccc}\lambda_{1}&0&0\\\
0&\ddots&0\\\ 0&0&\lambda_{d}\end{array}\right),\quad
A^{v}=\left(\begin{array}[]{ccc}\lambda_{d+1}&0&0\\\ 0&\ddots&0\\\
0&0&\lambda_{n}\end{array}\right).$

Note that under the scaling (A13), the anchor trajectory $x^{*}(t)$ becomes
$\epsilon$-dependent, which is reflected by our modified notation
$x_{\epsilon}^{*}(t)$ for the same trajectory.

By definition, we have

$x_{0}^{*}(t)\equiv 0,$

and

$f_{0}\left(x_{\epsilon}^{*}(t)\right)+Ax_{\epsilon}^{*}(t)-\dot{x}_{\epsilon}^{*}(t)+\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t),t\right)\equiv
0,$

therefore

$D_{\epsilon}^{p}\left[f_{0}\left(x_{\epsilon}^{*}(t)\right)+Ax_{\epsilon}^{*}(t)-\dot{x}_{\epsilon}^{*}(t)+\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t),t\right)\right]=D_{\epsilon}^{p}\hat{f}(0,0,\epsilon;t)\equiv
0,\quad p\in\mathbb{N},\quad t\in\mathbb{R}.$ (B38)

By the definition of $\hat{f}$ in (B37) and by formula (B38) with $p=1$ , we
also have

$\hat{f}(0,0,0;t)=0,\quad D_{u}\hat{f}(0,0,0;t)=0,\quad
D_{v}\hat{f}(0,0,0;t)=0,\quad D_{\epsilon}\hat{f}(0,0,0;t)=0.$ (B39)

For $\epsilon=0$, under the nonresonance conditions (17), we have a unique,
primary spectral submanifold

$W_{0}\left(E\right)=\left\\{\left(u,v\right)\in
U\subset\mathbb{R}^{n}:\,v=h_{0}(u)=\sum_{\left|\mathbf{k}\right|\geq
0}h^{\mathbf{k}}u^{\mathbf{k}}\right\\},$ (B40)

where $h_{0}(u)=\mathcal{O}\left(\left|u\right|^{2}\right)\in C^{r}$ defines
the primary SSM as a smooth graph over $E$ with a quadratic tangency to $E$ at
$x=0$ (see Haller and Ponsioen [18]). For $\epsilon>0$ small, under the
conditions of statement (i) of Theorem 3, this primary SSM persist in the form
of a (generally non-unique) invariant manifold

$W_{\epsilon}\left(E,t\right)=\left\\{\left(u,v\right)\in
U\subset\mathbb{R}^{n}:\,v=h_{\epsilon}(u,t)=\sum_{\left|\mathbf{k}\right|,p\geq
0}h^{\mathbf{k}p}(t)u^{\mathbf{k}}\epsilon^{p}\right\\},$ (B41)

with coefficients $h^{\mathbf{k}p}(t)$ that are uniformly bounded in $t$
within a $B_{\delta}$ ball around the origin $x=0$ for
$0\leq\epsilon\leq\delta$, and with the notation

$u^{\mathbf{k}}=u_{1}^{k_{1}}\cdot\ldots\cdot u_{d}^{k_{d}}.$

Note that this expansion is only valid up to the order of smoothness of the
SSM. For general members of the $W_{\epsilon}\left(E,t\right)$, this
smoothness can only guaranteed to be class $C^{\rho}$, but exceptional members
may admit higher-order expansions, just as primary SSMs do in the case of
autonomous, time-periodic and time-quasiperiodic SSMs (see Haller and Ponsioen
[18]).

Also note that for the smooth persistence of $W_{0}(E)$ as
$W_{\epsilon}\left(E,t\right)$, a possible $1\colon 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. (62) from $2$ to $1$
relative to what one requires for the existence of $W_{0}\left(E\right)$.

Differentiating the definition (B41) of the invariant manifold
$W_{\epsilon}\left(E,t\right)$ in time and using the system of ODEs (B32), we
obtain

$\displaystyle\dot{v}$
$\displaystyle=D_{u}h_{\epsilon}(u,t)\dot{u}+D_{t}h_{\epsilon}(u,t)$
$\displaystyle=D_{u}h_{\epsilon}(u,t)\left[A^{u}u+\hat{f}^{u}(u,h_{\epsilon}(u,t),t;\epsilon)\right]+D_{t}h_{\epsilon}(u,t)$
(B42)

At the same time, substitution into (B32) gives

$\displaystyle\dot{v}$
$\displaystyle=A^{v}h_{\epsilon}(u,t)+\hat{f}^{v}(u,h_{\epsilon}(u,t),t;\epsilon).$
(B43)

Comparing (B42) and (B43) give the invariance PDE satisfied by
$W_{\epsilon}\left(E,t\right)$:

$D_{u}h_{\epsilon}\left[A^{u}u+\hat{f_{0}^{u}}+\epsilon\hat{f_{1}^{u}}\right]+D_{t}h_{\epsilon}=A^{v}h_{\epsilon}+\hat{f_{0}^{v}}+\epsilon\hat{f_{1}^{v}}.$
(B44)

As the origin $(u,v)=(0,0)$ is a fixed point of system (B32) for all $t$ and
$\epsilon$, we must have $h(0,t;\epsilon)\equiv 0$ for all $\epsilon\geq 0$
small and all $t\in\mathbb{R}$. This, in turn, implies,

$h^{\mathbf{0}p}(t)\equiv 0,\quad t\in\mathbb{R},\quad p\in\mathbb{N}.$ (B45)

Furthermore, for $\epsilon=0$, the non-autonomous SSM
$W_{\epsilon}\left(E,t\right)$ becomes the autonomous SSM $W(E)$, and hence we
must have

$h^{\mathbf{k}0}(t)\equiv h^{\mathbf{k}},\quad
t\in\mathbb{R},\quad\mathbf{k}\in\mathbb{N}^{d},\qquad
h^{\mathbf{k}0}(t)\equiv h^{\mathbf{k}}=0,\quad\left|\mathbf{k}\right|=1,$
(B46)

with the time-independent Taylor series coefficients $h^{\mathbf{k}}$ in the
expansion for $W_{0}(E)$ in eq. (B40).

####  B.2.2 Structure of the invariance equation

Substitution of (B41) into the invariance equation (B44) gives

$\sum_{\left|\mathbf{k}\right|,p\geq
0}h^{\mathbf{k}p}(t)\epsilon^{p}D_{u}u^{\mathbf{k}}\left[A^{u}u+\hat{f}^{u}\left(u,\sum_{\left|\mathbf{k}\right|,p\geq
0}h^{\mathbf{k}p}(t)u^{\mathbf{k}}\epsilon^{p},\epsilon,t\right)\right]+\sum_{\left|\mathbf{k}\right|,p\geq
0}\dot{h}^{\mathbf{k}p}(t)u^{\mathbf{k}}\epsilon^{p}$

$=\sum_{\left|\mathbf{k}\right|\geq 0,p\geq
0}A^{v}h^{\mathbf{k}p}(t)u^{\mathbf{k}}\epsilon^{p}+\hat{f}^{v}\left(u,\sum_{\left|\mathbf{k}\right|,p\geq
0}h^{\mathbf{k}p}(t)u^{\mathbf{k}}\epsilon^{p},\epsilon,t\right).$ (B47)

We observe that

$\displaystyle
A^{v}h^{\mathbf{k}p}(t)u^{\mathbf{k}}-h^{\mathbf{k}p}(t)D_{u}u^{\mathbf{k}}A^{u}u$
$\displaystyle=\left(\begin{array}[]{ccc}\lambda_{d+1}&0&0\\\ 0&\ddots&0\\\
0&0&\lambda_{n}\end{array}\right)\left(\begin{array}[]{c}h_{1}^{\mathbf{k}p}\\\
\vdots\\\ h_{n-d}^{\mathbf{k}p}\end{array}\right)u^{\mathbf{k}}$ (B54)
$\displaystyle-\left(\begin{array}[]{c}h_{1}^{\mathbf{k}p}\\\ \vdots\\\
h_{n-d}^{\mathbf{k}p}\end{array}\right)\left(\begin{array}[]{ccc}\frac{k_{1}u^{\mathbf{k}}}{u_{1}}&\cdots&\frac{k_{d}u^{\mathbf{k}}}{u_{d}}\end{array}\right)\left(\begin{array}[]{c}\lambda_{1}u_{1}\\\
\vdots\\\ \lambda_{d}u_{d}\end{array}\right)$ (B62)
$\displaystyle=\left(\begin{array}[]{c}\lambda_{d+1}h_{1}^{\mathbf{k}p}\\\
\vdots\\\
\lambda_{n}h_{n-d}^{\mathbf{k}p}\end{array}\right)u^{\mathbf{k}}-\left(\begin{array}[]{ccc}h_{1}^{\mathbf{k}p}\frac{k_{1}u^{\mathbf{k}}}{u_{1}}&\cdots&h_{1}^{\mathbf{k}p}\frac{k_{d}u^{\mathbf{k}}}{u_{d}}\\\
\vdots&\ddots&\vdots\\\
h_{n-d}^{\mathbf{k}p}\frac{k_{1}u^{\mathbf{k}}}{u_{1}}&\cdots&h_{n-d}^{\mathbf{k}p}\frac{k_{d}u^{\mathbf{k}}}{u_{d}}\end{array}\right)\left(\begin{array}[]{c}\lambda_{1}u_{1}\\\
\vdots\\\ \lambda_{d}u_{d}\end{array}\right)$ (B72)
$\displaystyle=\left(\begin{array}[]{c}\lambda_{d+1}h_{1}^{\mathbf{k}p}\\\
\vdots\\\
\lambda_{n}h_{n-d}^{\mathbf{k}p}\end{array}\right)u^{\mathbf{k}}-\left(\begin{array}[]{c}h_{1}^{\mathbf{k}p}\sum_{j=1}^{d}k_{j}\lambda_{j}\\\
\vdots\\\
h_{n-d}^{\mathbf{k}p}\sum_{j=1}^{d}k_{j}\lambda_{j}\end{array}\right)u^{\mathbf{k}}$
(B79)
$\displaystyle=\left(\begin{array}[]{ccc}\lambda_{d+1}-\sum_{j=1}^{d}k_{j}\lambda_{j}&\cdots&0\\\
\vdots&\ddots&\vdots\\\
0&\cdots&\lambda_{n}-\sum_{j=1}^{d}k_{j}\lambda_{j}\end{array}\right)h^{\mathbf{k}p}(t)u^{\mathbf{k}}.$
(B83)

Therefore, the invariance equation (B47) can be rewritten as

$\sum_{\left|\mathbf{k}\right|,p\geq
0}\dot{h}^{\mathbf{k}p}(t)u^{\mathbf{k}}\epsilon^{p}=\sum_{\left|\mathbf{k}\right|,p\geq
0}A_{\mathbf{k}}h^{\mathbf{k},p}(t)u^{\mathbf{k}}\epsilon^{p}+\sum_{\left|\mathbf{k}\right|,p\geq
0}M^{\mathbf{k}p}(t,h^{\mathbf{j}m})u^{\mathbf{k}}\epsilon^{p},$ (B84)

where

$\displaystyle A_{\mathbf{k}}$
$\displaystyle=\mathrm{diag}\left[\lambda_{\ell}-\sum_{j=1}^{d}k_{j}\lambda_{j}\right]_{\ell=d+1}^{n}\in\mathbb{C}^{\left(n-d\right)\times\left(n-d\right)},$
$\displaystyle\sum_{\left|\mathbf{k}\right|,p\geq
0}M^{\mathbf{k}p}(t,h^{\mathbf{j}m})u^{\mathbf{k}}\epsilon^{p}$
$\displaystyle=M(u,h^{\mathbf{j}m},\epsilon;t)=\hat{f}^{v}\left(u,\sum_{\left|\mathbf{k}\right|,p\geq
0}h^{\mathbf{k}p}(t)u^{\mathbf{k}}\epsilon^{p},\epsilon;t\right)$ (B85)
$-\sum_{\left|\mathbf{k}\right|,p\geq
0}\epsilon^{p}\left(\begin{array}[]{ccc}h_{1}^{\mathbf{k}p}(t)\frac{k_{1}u^{\mathbf{k}}}{u_{1}}&\cdots&h_{1}^{\mathbf{k}p}(t)\frac{k_{d}u^{\mathbf{k}}}{u_{d}}\\\
\vdots&\ddots&\vdots\\\
h_{n-d}^{\mathbf{k}p}(t)\frac{k_{1}u^{\mathbf{k}}}{u_{1}}&\cdots&h_{n-d}^{\mathbf{k}p}(t)\frac{k_{d}u^{\mathbf{k}}}{u_{d}}\end{array}\right)\hat{f}^{u}\left(u,\sum_{\left|\mathbf{k}\right|,p\geq
0}h^{\mathbf{k}p}(t)u^{\mathbf{k}}\epsilon^{p},\epsilon;t\right).$

Recall from eq. (B39) that $\hat{f}^{v}$ and $\hat{f}^{u}$ vanish for
$u,v,\epsilon=0$ and have no $\mathcal{O}\left(\left|v\right|\right)$ terms.
Therefore, when $\sum_{\left|\mathbf{k}\right|,p\geq
0}h^{\mathbf{k}p}(t)u^{\mathbf{k}}\epsilon^{p}$ is substituted into the terms
in the Taylor expansion of $\hat{f}^{v}\left(u,v,\epsilon;t\right)$, then it
is multiplied in each case by at least the first power of $u$ or at least the
first power of $\epsilon.$ (When $\sum_{\left|\mathbf{k}\right|,p\geq
0}h^{\mathbf{k}p}(t)u^{\mathbf{k}}\epsilon^{p}$ is substituted into terms of
$\mathcal{O}\left(\left|v\right|^{2}\right)$ then it is multiplied by itself,
but $h^{\mathbf{0}0}(t)\equiv 0$, and hence the lowest order term multiplying
$\sum_{\left|\mathbf{k}\right|,p\geq
0}h^{\mathbf{k}p}(t)u^{\mathbf{k}}\epsilon^{p}$ will be again at least the
first power of $u$ or at least the first power of $\epsilon.$) As a result,
$M^{\mathbf{k}p}(h^{\mathbf{j}m},t)$ will only depend on $h^{\mathbf{j}m}$
that are lower in order, i.e.,

$D_{h^{\mathbf{j}m}}M^{\mathbf{k}p}(h^{\mathbf{j}m},t)=0,\quad\left|(\mathbf{j},m)\right|\geq\left|(\mathbf{k},p)\right|.$

Equating coefficients of equal powers of $u$ in eq. (B84), we therefore obtain
the recursively solvable linear system of inhomogeneous linear ODEs

$\dot{h}^{\mathbf{k}p}(t)=A_{\mathbf{k}}h^{\mathbf{k}p}(t)+M^{\mathbf{k}p}(t,h^{\mathbf{j}m}),\quad\left|(\mathbf{j},m)\right|<\left|(\mathbf{k},p)\right|,$
(B86)

with $M^{\mathbf{k}p}(t,h^{\mathbf{j}m})$ defined in formula (66). So far, we
know from eqs. (B45) and (B46) that

$h^{\mathbf{0}p}(t)\equiv 0,\quad t\in\mathbb{R},\quad p\in\mathbb{N},\qquad
h^{\mathbf{k}0}(t)\equiv h^{\mathbf{k}},\quad
t\in\mathbb{R},\quad\mathbf{k}\in\mathbb{N}^{d},$ (B87)

which is a homogeneous system of ODEs in (B86) for
$\left|\mathbf{k}\right|>0$.

####  B.2.3 Solution of the invariance equation

A first observation about solving the family (B86) of linear ODEs is that for
arbitrary $\left|\mathbf{k}\right|$ and $p=0$, we obtain from eqs. (B46) and
(B86) the algebraic equations

$\displaystyle 0$
$\displaystyle=A_{\mathbf{k}}h^{\mathbf{k}0}+M^{\mathbf{k}0}(h^{\mathbf{j}0}).$

We know from classic SSM theory that this algebraic system of equations is a
recursively solvable linear system of equations, because
$\left|\mathbf{j}\right|<\left|\mathbf{k}\right|$ holds for all
$M^{\mathbf{k}0}(h^{\mathbf{j}0})$, and hence

$h^{\mathbf{k}0}=-A_{\mathbf{k}}^{-1}M^{\mathbf{k}0}(h^{\mathbf{j}0}),\quad\left|\mathbf{j}\right|<\left|\mathbf{k}\right|.$
(B88)

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

Under the nonresonance condition (61), all diagonal elements of
$A_{\mathbf{k}}$ have nonzero real parts. We can therefore uniquely split
$A_{\mathbf{k}}$ into the sum of a diagonal matrix $A_{\mathbf{k}}^{-}$ which
contains the stable eigenvalues of $A_{\mathbf{k}}$ as well as zeros, and a
diagonal matrix $A_{\mathbf{k}}^{+}$ which contains the unstable eigenvalues
of $A_{\mathbf{k}}$ as well as zeros:

$\displaystyle A_{\mathbf{k}}$
$\displaystyle=A_{\mathbf{k}}^{-}+A_{\mathbf{k}}^{+},$ $\displaystyle
A_{\mathbf{k}}^{-}(j,j)$
$\displaystyle=\left\\{\begin{array}[]{ccc}A_{\mathbf{k}}(j,j),&&\mathrm{Re}A_{\mathbf{k}}(j,j)<0,\\\
\\\ 0,&&\mathrm{otherwise},\end{array}\right.\quad
A_{\mathbf{k}}^{+}(j,j)=\left\\{\begin{array}[]{ccc}A_{\mathbf{k}}(j,j),&&\mathrm{Re}A_{\mathbf{k}}(j,j)>0,\\\
\\\ 0,&&\mathrm{otherwise}.\end{array}\right.$

We define the corresponding splitting for the vectors $h^{\mathbf{k}p}(t)$ and
$M^{\mathbf{k}p}(t,h^{\mathbf{j}m}(t))$ as

$\displaystyle h^{\mathbf{k}p}$
$\displaystyle=h^{\mathbf{k}p-}+h^{\mathbf{k}p+},\quad$ $\displaystyle
h_{j}^{\mathbf{k}p-}$
$\displaystyle=\left\\{\begin{array}[]{ccc}h_{j}^{\mathbf{k}p},&&\mathrm{Re}A_{\mathbf{k}}(j,j)<0,\\\
\\\ 0,&&\mathrm{otherwise},\end{array}\right.\quad
h_{j}^{\mathbf{k}p+}=\left\\{\begin{array}[]{ccc}h_{j}^{\mathbf{k}p},&&\mathrm{Re}A_{\mathbf{k}}(j,j)>0,\\\
\\\ 0,&&\mathrm{otherwise},\end{array}\right.$ $\displaystyle M^{\mathbf{k}p}$
$\displaystyle=M^{\mathbf{k}p-}+M^{\mathbf{k}p+},$ $\displaystyle
M_{j}^{\mathbf{k}p-}$
$\displaystyle=\left\\{\begin{array}[]{ccc}M_{j}^{\mathbf{k}p},&&\mathrm{Re}A_{\mathbf{k}}(j,j)<0,\\\
\\\ 0,&&\mathrm{otherwise},\end{array}\right.\quad
M_{j}^{\mathbf{k}p+}=\left\\{\begin{array}[]{ccc}M_{j}^{\mathbf{k}p},&&\mathrm{Re}A_{\mathbf{k}}(j,j)>0,\\\
\\\ 0,&&\mathrm{otherwise}.\end{array}\right..$

We then then split the expression for the solution of system (B86) as

$\displaystyle h^{\mathbf{k}p-}(t)$
$\displaystyle=e^{A_{\mathbf{k}}^{-}\left(t-t_{0}\right)}h^{\mathbf{k}p-}(t_{0})+\int_{t_{0}}^{t}e^{A_{\mathbf{k}}^{-}\left(t-s\right)}M^{\mathbf{k}p-}\left(s,h^{\mathbf{j}m}(s)\right)\,ds,$
$\displaystyle h^{\mathbf{k}p+}(t)$
$\displaystyle=e^{A_{\mathbf{k}}^{+}\left(t-t_{0}\right)}h^{\mathbf{k}p+}(t_{0})+\int_{t_{0}}^{t}e^{A_{\mathbf{k}}^{+}\left(t-s\right)}M^{\mathbf{k}p+}\left(s,h^{\mathbf{j}m}(s)\right)\,ds.$
(B89)

If $h^{\mathbf{k}p}(t)$ is a uniformly bounded solution of system (B86) for
all forward and backward times, then using the signs of the nonzero diagonal
entries of the matrices $A_{\mathbf{k}}^{\pm}$, we obtain from the formulas
(B89) that

$\displaystyle h^{\mathbf{k}p-}(t)$
$\displaystyle=\lim_{t_{0}\to-\infty}\left[e^{A_{\mathbf{k}}^{-}\left(t-t_{0}\right)}h^{\mathbf{k}p-}(t_{0})+\int_{t_{0}}^{t}e^{A_{\mathbf{k}}^{-}\left(t-s\right)}M^{\mathbf{k}p-}\left(s,h^{\mathbf{j}m}(s)\right)\,ds\right]$
$\displaystyle=\int_{-\infty}^{t}e^{A_{\mathbf{k}}^{-}\left(t-s\right)}M^{\mathbf{k}p-}(s,h^{\mathbf{j}m}(s))\,ds,$
$\displaystyle h^{\mathbf{k}p-}(t)$
$\displaystyle=\lim_{t_{0}\to+\infty}\left[e^{A_{\mathbf{k}}^{+}\left(t-t_{0}\right)}h^{\mathbf{k}+}(t_{0})+\int_{t_{0}}^{t}e^{A_{\mathbf{k}}^{+}\left(t-s\right)}M^{\mathbf{k}p+}\left(s,h^{\mathbf{j}m}(s)\right)\,ds\right]$
$\displaystyle=-\int_{t}^{+\infty}e^{A_{\mathbf{k}}^{+}\left(t-s\right)}e^{A_{\mathbf{k}}^{+}\left(t-s\right)}M^{\mathbf{k}p+}\left(s,h^{\mathbf{j}m}(s)\right)\,ds,$

which gives

$\displaystyle h^{\mathbf{k}p}(t)$
$\displaystyle=h^{\mathbf{k}p-}(t)+h^{\mathbf{k}p-}(t)=\int_{-\infty}^{t}e^{A_{\mathbf{k}}^{-}\left(t-s\right)}M^{\mathbf{k}p-}(s,h^{\mathbf{j}m}(s))\,ds$
$\displaystyle-\int_{t}^{+\infty}e^{A_{\mathbf{k}}^{+}\left(t-s\right)}e^{A_{\mathbf{k}}^{+}\left(t-s\right)}M^{\mathbf{k}p+}\left(s,h^{\mathbf{j}m}(s)\right)\,ds.$

Therefore, introducing the Green’s function as the diagonal matrix
$G_{\mathbf{k}}(t)\in\mathbb{C}^{\left(n-d\right)\times\left(n-d\right)}$
defined in formula (56), we can write the unique globally bounded solution of
eq. (B86) in the form

$h^{\mathbf{k}p}(t)=\int_{-\infty}^{\infty}G_{\mathbf{k}}(t-s)M^{\mathbf{k}p}(s,h^{\mathbf{j}m}(s))\,ds,\quad\left|(\mathbf{j},m)\right|<\left|(\mathbf{k},p)\right|,\quad
t\in\mathbb{R},$ (B90)

with $M^{\mathbf{k}p}(t,h^{\mathbf{j}m}(t))$ defined in formula (66). As we
have already seen in eqs. B87) and (B88), we specifically have

$h^{\mathbf{0}p}(t)\equiv 0,\quad t\in\mathbb{R},\quad p\in\mathbb{N},\qquad
h^{\mathbf{k}0}(t)\equiv-
A_{\mathbf{k}}^{-1}M^{\mathbf{k}0}(h^{\mathbf{j}0}),\quad\left|\mathbf{j}\right|<\left|\mathbf{k}\right|,\quad
h^{\mathbf{0}}=0.$ (B91)

####  B.2.4 Reduced dynamics

We obtain the reduced dynamics on $W_{\epsilon}\left(E,t\right)$ by
restricting the $u$-component of system (B32) to
$W_{\epsilon}\left(E,t\right)$:

$\dot{u}=A^{u}u+\hat{f}^{u}\left(u,h_{\epsilon}\left(u,t\right)\right).$ (B92)

We obtain the form (67) of the reduced dynamics by substituting the definition
of $\hat{f}^{u}$ into eq. (B92), using the defining relationship

$\dot{x}_{\epsilon}^{*}(t)=Ax_{\epsilon}^{*}(t)+f_{0}\left(x_{\epsilon}^{*}(t)\right)+\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t),t\right)$

of the anchor trajectory $x_{\epsilon}^{*}(t)$, and noting that the rows of
$P^{-1}$ are the appropriately normalized left eigenvectors of $A$, and hence
$Q_{u}\in\mathbb{C}^{d\times n}$ in formula (67) contains the first $d$ rows
of $P^{-1}$. The end result is then:__

$\displaystyle\dot{u}$
$\displaystyle=A^{u}u+\hat{f}^{u}\left(u,h_{\epsilon}(u,t),\epsilon;t\right)$
$\displaystyle=A^{u}u+Q_{u}\left[f_{0}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
h_{\epsilon}(u,t)\end{array}\right)\right)+\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
h_{\epsilon}(u,t)\end{array}\right),t\right)+Ax_{\epsilon}^{*}(t)-\dot{x}_{\epsilon}^{*}(t)\right]$
$\displaystyle=A^{u}u$
$\displaystyle+Q_{u}\left[f_{0}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
h_{\epsilon}(u,t)\end{array}\right)\right)-f_{0}\left(x_{\epsilon}^{*}(t)\right)+\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
h_{\epsilon}(u,t)\end{array}\right),t\right)-\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t),t\right)\right],$

as claimed in statement (ii) of Theorem 4.

The $\left(u,v\right)$ coordinates are measured from the perturbed anchor
trajectory $x^{*}(t)$. We can also express the reduced dynamics on
$W_{\epsilon}\left(E,t\right)$ in coordinates that are aligned with the
subspace $E$ and emanate from the original $x=0$ fixed point of system (1) for
$\epsilon=0.$ Let

$\left(\begin{array}[]{c}\xi\\\ \eta\end{array}\right)=P^{-1}x,$

which gives

$\left(\begin{array}[]{c}u\\\
v\end{array}\right)=P^{-1}\left(x-x_{\epsilon}^{*}(t)\right)=\left(\begin{array}[]{c}\xi\\\
\eta\end{array}\right)-P^{-1}x_{\epsilon}^{*}(t),$ (B93)

which implies

$\xi=u+Q_{u}x_{\epsilon}^{*}(t),\quad\eta=v+Q_{v}x_{\epsilon}^{*}(t)\qquad
P^{-1}=\left[\begin{array}[]{c}Q_{u}\\\ Q_{v}\end{array}\right],\quad
Q_{u}\in\mathbb{C}^{d\times n},\quad Q_{u}\in\mathbb{C}^{(n-d)\times n}.$

Note that

$P^{-1}AP=\left[\begin{array}[]{cc}A^{u}&0_{d\times n}\\\
0_{\left(n-d\right)\times n}&A^{v}\end{array}\right],$

which implies

$\left[\begin{array}[]{c}Q_{u}A\\\
Q_{v}A\end{array}\right]=P^{-1}A=\left[\begin{array}[]{cc}A^{u}&0_{d\times
n}\\\ 0_{\left(n-d\right)\times
n}&A^{v}\end{array}\right]P^{-1}=\left[\begin{array}[]{cc}A^{u}&0_{d\times
n}\\\ 0_{\left(n-d\right)\times
n}&A^{v}\end{array}\right]\left[\begin{array}[]{c}Q_{u}\\\
Q_{v}\end{array}\right]=\left[\begin{array}[]{c}A^{u}Q_{u}\\\
A^{v}Q_{v}\end{array}\right],$

yielding the identity

$Q_{u}A=A^{u}Q_{u}.$ (B94)

Given that __

$\displaystyle\dot{u}$
$\displaystyle=A^{u}u+Q_{u}\left[f_{0}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
h_{\epsilon}(u,t)\end{array}\right)\right)+\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
h_{\epsilon}(u,t)\end{array}\right),t\right)+Ax_{\epsilon}^{*}(t)-\dot{x}_{\epsilon}^{*}(t)\right],$
(B99)

we obtain

$\displaystyle\dot{\xi}$
$\displaystyle=\dot{u}+Q_{u}\dot{x}_{\epsilon}^{*}(t)$
$\displaystyle=A^{u}u+Q_{u}\left[f_{0}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
h_{\epsilon}(u,t)\end{array}\right)\right)+\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
h_{\epsilon}(u,t)\end{array}\right),t\right)+Ax_{\epsilon}^{*}(t)\right]$
$\displaystyle=A^{u}\left(\xi-
Q_{u}x_{\epsilon}^{*}(t)\right)+Q_{u}Ax_{\epsilon}^{*}(t)$
$\displaystyle+Q_{u}\left[f_{0}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
h_{\epsilon}(u,t)\end{array}\right)\right)+\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}u\\\
h_{\epsilon}(u,t)\end{array}\right),t\right)\right]$
$\displaystyle=A^{u}\xi+\left[Q_{u}A-A^{u}Q_{u}\right]x_{\epsilon}^{*}(t)$
$\displaystyle+Q_{u}\left[f_{0}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}\xi-
Q_{u}x_{\epsilon}^{*}(t)\\\ h_{\epsilon}(\xi-
Q_{u}x_{\epsilon}^{*}(t),t)\end{array}\right)\right)+\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}\xi-
Q_{u}x_{\epsilon}^{*}(t)\\\ h_{\epsilon}(\xi-
Q_{u}x_{\epsilon}^{*}(t),t)\end{array}\right),t\right)\right]$
$\displaystyle=A^{u}\xi+Q_{u}f_{0}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}\xi-
Q_{u}x_{\epsilon}^{*}(t)\\\ h_{\epsilon}(\xi-
Q_{u}x_{\epsilon}^{*}(t),t)\end{array}\right)\right)$
$\displaystyle+Q_{u}\epsilon\tilde{f}_{1}\left(x_{\epsilon}^{*}(t)+P\left(\begin{array}[]{c}\xi-
Q_{u}x_{\epsilon}^{*}(t)\\\ h_{\epsilon}(\xi-
Q_{u}x_{\epsilon}^{*}(t),t)\end{array}\right),t\right)$

where we have used the identity (B94). This completes the proof of statement
(iii) of Theorem 4.

##  Appendix C Proof of Theorems 5-7

###  C.1 Proof of Theorem 5: Expansion for the hyperbolic slow manifold
$\mathcal{L}_{\epsilon}$

We now assume that $f$ and $x_{0}^{\prime}(\alpha)$ have $r\geq 1$ uniformly
bounded derivatives over a closed neighborhood of $x_{0}(\alpha)$ for all
$\alpha$. Then, by the results of Eldering [12] on the persistence of non-
compact, normally hyperbolic invariant manifolds, a unique and uniformly
bounded slow manifold $\mathcal{L}_{\epsilon}$ exists near $\mathcal{L}_{0}$
for $\epsilon$ small enough. These results are applicable because the
$\mathcal{O}\left(\epsilon\right)$ perturbation term in (107) is uniformly
bounded in the $C^{1}$ norm. The slow manifold $\mathcal{L}_{\epsilon}$ is
$C^{r}$-diffeomorphic to $\mathcal{L}_{0}$ and is as smooth in the $\epsilon$
parameter as system (107). As a consequence, we have an asymptotic expansion

$\mathcal{L}_{\epsilon}=\left\\{(x,\alpha)\in\mathbb{R}^{n}\times\mathbb{R}\colon\,\,\,\,\,x=x_{\epsilon}(\alpha)=\sum_{j=0}^{r}\epsilon^{j}x_{j}(\alpha)+o\left(\epsilon^{r}\right)\right\\},$
(C1)

where the functions $x_{j}(\alpha)$ are uniformly bounded in $\alpha$. By the
definition of $x_{\epsilon}(\alpha)$ and by eq. (110), we have

$D_{\epsilon}^{p}\left[f\left(x_{\epsilon}(\alpha),\alpha\right)-\epsilon
x_{\epsilon}^{\prime}(\alpha)\right]\equiv 0,\quad
p\in\mathbb{N},\quad\alpha\in\mathbb{R}.$ (C2)

Substituting the expansion (C1) into (107), we obtain

$\sum_{j\geq 0}\epsilon^{j+1}x_{j}^{\prime}(\alpha)=f\left(\sum_{j\geq
0}\epsilon^{j}x_{j}(\alpha),\alpha\right).$

Comparison of equal powers of $\epsilon$ in the last equation gives

$x_{\nu-1}^{\prime}(\alpha)=\left.\frac{1}{\nu!}\frac{\partial^{\nu}}{\partial\epsilon^{\nu}}f\left(\sum_{j\geq
0}\epsilon^{j}x_{j}(\alpha),\alpha\right)\right|_{\epsilon=0},\quad\nu\geq 1.$
(C3)

Specifically, we have

$\epsilon
x_{0}^{\prime}(\alpha)+\epsilon^{2}x_{1}^{\prime}(\alpha)+\epsilon^{3}x_{2}^{\prime}(\alpha)+\ldots=f\left(x_{0}(\alpha),\alpha\right)+\epsilon
Dfx_{1}(\alpha)+\left.\frac{1}{2}\frac{\partial^{2}}{\partial\epsilon^{2}}f\left(\sum_{p\geq
0}\epsilon^{p}x_{p}(\alpha),\alpha\right)\right|_{\epsilon=0}\epsilon^{2}+\ldots.$

Note that

$\displaystyle\frac{\partial^{2}}{\partial\epsilon^{2}}f\left(\sum_{j\geq
0}\epsilon^{j}x_{j}(\alpha),\alpha\right)$
$\displaystyle=\frac{\partial}{\partial\epsilon}\left[D_{x}f\left(\sum_{j\geq
0}\epsilon^{j}x_{j}(\alpha),\alpha\right)\sum_{j\geq
0}j\epsilon^{j-1}x_{j}(\alpha)\right]$
$\displaystyle=D_{x}^{2}f\left(\sum_{j\geq
0}\epsilon^{j}x_{j}(\alpha),\alpha\right)\otimes\left(\sum_{j\geq
1}j\epsilon^{j-1}x_{j}(\alpha)\right)\otimes\left(\sum_{j\geq
1}j\epsilon^{j-1}x_{j}(\alpha)\right)$ $\displaystyle+D_{x}f\left(\sum_{j\geq
0}\epsilon^{j}x_{j}(\alpha),\alpha\right)\sum_{j\geq
2}j\left(j-1\right)\epsilon^{j-2}x_{j}(\alpha),$

therefore, we have

$\left.\frac{1}{2}\frac{\partial^{2}}{\partial\epsilon^{2}}f\left(\sum_{j\geq
0}\epsilon^{j}x_{j}(\alpha),\alpha\right)\right|_{\epsilon=0}=\frac{1}{2}D_{x}^{2}f\left(x_{0}(\alpha),\alpha\right)\otimes
x_{1}(\alpha)\otimes
x_{1}(\alpha)+D_{x}f\left(x_{0}(\alpha),\alpha\right)x_{2}(\alpha)$

which gives formulas (115) as solutions up to second order. We note that
$A^{-1}(\alpha)=\left[D_{x}f\left(x_{0}(\alpha),\alpha\right)\right]^{-1}$ is
known to exist by the hyperbolicity assumption (109).

To obtain the coefficients in eq. (C1) up to an arbitrary order $\nu$, we need
to solve the recursive set of algebraic equations (C3) for $x_{j}(\alpha)$. To
solve these equations, we again recall the multi-variate Faá di Bruno formula
(A21), for which we now have $p=1$, wherein we have $H=f_{q}$,
$g^{i}=\sum_{p\geq 0}\epsilon^{p}x_{p}^{i}(\alpha)$,
$\mathbf{x}=\epsilon\in\mathbb{R},$ $\mathbf{x}^{0}=0\in\mathbb{R}$,
$\mathbf{\boldsymbol{\nu}}=\nu\in\mathbb{N}$,
$\boldsymbol{\ell}_{i}=\ell_{i}\in\mathbb{N}$, and
$\mathbf{y}=x\in\mathbb{R}^{n}$. Therefore, for

$H(\epsilon)=f_{q}\left(\sum_{j\geq
0}\epsilon^{j}x_{j}^{1}(\alpha),\ldots,\sum_{j\geq
0}\epsilon^{j}x_{j}^{n}(\alpha),\alpha\right),\quad q=1,\ldots,n,$

formula (A21) gives

$\displaystyle\frac{d^{\nu}}{d\epsilon^{\nu}}H\left(0\right)$
$\displaystyle=\left.D_{\epsilon}^{\nu}f_{q}\left(\sum_{j\geq
0}\epsilon^{j}x_{j}^{1}(\alpha),\ldots,\sum_{j\geq
0}\epsilon^{j}x_{j}^{n}(\alpha),\alpha\right)\right|_{\epsilon=0}$
$\displaystyle=\sum_{1\leq\left|\boldsymbol{\gamma}\right|\leq\nu}D_{x}^{\mathbf{\boldsymbol{\gamma}}}f_{q}\left(x_{0}(\alpha),\alpha\right)\sum_{s=1}^{\nu}\sum_{p_{s}\left(\nu,\boldsymbol{\gamma}\right)}\nu!\prod_{j=1}^{s}\frac{\prod_{i=1}^{n}\left[\left(\ell_{j}\right)!x_{\ell_{j}}^{i}(\alpha)\right]^{k_{ji}}}{\left(\mathbf{k}_{j}\right)!\left[\left(\ell_{j}\right)!\right]^{\left|\mathbf{k}_{j}\right|}}$
$\displaystyle=\sum_{1\leq\left|\boldsymbol{\gamma}\right|\leq\nu}D_{x}^{\boldsymbol{\gamma}}f_{q}\left(x_{0}(\alpha),\alpha\right)\sum_{s=1}^{\nu}\sum_{p_{s}\left(\nu,\boldsymbol{\gamma}\right)}\nu!\prod_{j=1}^{s}\frac{\prod_{i=1}^{n}\left[x_{\ell_{j}}^{i}(\alpha)\right]^{k_{ji}}}{\prod_{i=1}^{n}k_{ji}!}$
$\displaystyle=\nu!\left[A(\alpha)x_{\nu}(\alpha)\right]_{q}$
$\displaystyle+\sum_{1<\left|\boldsymbol{\gamma}\right|\leq\nu}D_{x}^{\boldsymbol{\gamma}}f_{q}\left(x_{0}(\alpha),\alpha\right)\sum_{s=1}^{\nu}\sum_{p_{s}\left(\nu,\boldsymbol{\gamma}\right)}\nu!\prod_{j=1}^{s}\frac{\prod_{i=1}^{n}\left[x_{\ell_{j}}^{i}(\alpha)\right]^{k_{ji}}}{\prod_{i=1}^{n}k_{ji}!},\quad
q=1,\ldots,n.$

Using this last expression in (C3), we obtain formula (114).

###  C.2 Existence of the adiabatic SSM $\mathcal{M}_{\epsilon}$

The classic persistence results of Fenichel [14] assume compactness for the
underlying manifold and hence do not guarantee the smooth persistence of
$\mathcal{M}_{0}$ for $\epsilon>0$. To conclude the persistence of
$\mathcal{M}_{0}$, we first employ the “wormhole” construct from Proposition
B1 of Eldering et al. [13] that extends $\mathcal{M}_{0}$ smoothly over its
boundary so that it becomes a subset of a $\rho$-normally hyperbolic, normally
attracting, class $C^{r}$ invariant manifold $\mathcal{\tilde{M}}_{0}$ without
boundary. Under this extension, the stable foliation of
$W^{s}\left(\mathcal{M}_{0}\right)$ coincides with that part of the stable
foliation of $W^{s}\left(\mathcal{\tilde{M}}_{0}\right)$.

Due to the exclusion of a $1:1$ resonance (see eq. (136)) between eigenvalues
inside and outside $\mathcal{\tilde{M}}_{0}$, , the non-compact, boundaryless
extended manifold $\mathcal{\tilde{M}}_{0}$ is a normally attracting invariant
manifold. Its persistence can then be concluded for small enough $\epsilon>0$
from related results by Eldering [12] as long as assumption (117) holds for
the same $\rho$, uniformly in $\alpha$. In addition, $\mathcal{M}_{0}$ and
$f(x,\alpha)$ must have $r$ derivatives that are uniformly bounded in $\alpha$
in a small neighborhood of $\mathcal{M}_{0}.$ For $\epsilon>0$ small enough,
we then obtain a unique (for a given choice of the wormhole construct)
persistent invariant manifold $\mathcal{M}_{\epsilon}$ that is diffeomorphic
to $\mathcal{M}_{0}$, has $r$ uniformly bounded derivatives and is
$\mathcal{O}\left(\varepsilon\right)$ $C^{1}$ -close to $\mathcal{M}_{0}$. The
smoothness class of $\mathcal{M}_{\epsilon}$ is $C^{m}$ where
$m=\min\left(r,\rho\right)$.

We do not obtain uniqueness from any of these constructs as they all involve
modifications of the vector field. But the initial conditions are anyway with
probability zero on an SSM that is smoother than the smoothness implied by the
spectral gap. In other words, there is an inherent non-uniqueness in the
choice of $\mathcal{M}_{\epsilon}$ as a $\rho$-normally hyperbolic invariant
manifold tangent to $E(\alpha)$ for each $\alpha$, as there are infinitely
many different choices for $W\left(E(\alpha)\right)$ to begin with.

Despite the non-uniqueness of $\mathcal{M}_{\epsilon}$, all persisting
manifolds $\mathcal{M}_{\epsilon}$ must contain the unique, persisting
continuation $\mathcal{L}_{\epsilon}$ of $\mathcal{L}_{0}$. The reason is that
by the results of Eldering [12], $\mathcal{L}_{\epsilon}$ is unique, uniformly
bounded and lies fully in a small, inflowing neighborhood of
$\mathcal{M}_{\epsilon}$ whose size is $\mathcal{O}(1)$ in $\epsilon$. Points
on $\mathcal{L}_{\epsilon}$ would then be mapped by the inverse flow map
outside that inflowing neighborhood, unless they are contained in
$\mathcal{M}_{\epsilon}$. Therefore,
$\mathcal{L}_{\epsilon}\subset\mathcal{M}_{\epsilon}$ must hold.

###  C.3 Computation of the adiabatic SSM $\mathcal{M}_{\epsilon}$

Based on formula (C2),

$\displaystyle\hat{f}(u,v,\epsilon;\alpha)$
$\displaystyle=P^{-1}\left(\alpha\right)\left[f\left(x_{\epsilon}(\alpha)+P\left(\alpha\right)\left(\begin{array}[]{c}u\\\
v\end{array}\right),\alpha\right)-A\left(\alpha\right)P\left(\alpha\right)\left(\begin{array}[]{c}u\\\
v\end{array}\right)-\epsilon x_{\epsilon}^{\prime}(\alpha)-\epsilon
P^{\prime}\left(\alpha\right)\left(\begin{array}[]{c}u\\\
v\end{array}\right)\right]$ (C10)
$\displaystyle=P^{-1}\left(\alpha\right)\left[f\left(x_{\epsilon}(\alpha),\alpha\right)+D_{x}f\left(x_{\epsilon}(\alpha),\alpha\right)P\left(\alpha\right)\left(\begin{array}[]{c}u\\\
v\end{array}\right)\right.$ (C13)
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left.-A\left(\alpha\right)P\left(\alpha\right)\left(\begin{array}[]{c}u\\\
v\end{array}\right)-\epsilon x_{\epsilon}^{\prime}(\alpha)-\epsilon
P^{\prime}\left(\alpha\right)\left(\begin{array}[]{c}u\\\
v\end{array}\right)+\mathcal{O}\left(\left|u\right|^{2},\left|u\right|\left|v\right|,\left|v\right|^{2}\right)\right]$
(C18)
$\displaystyle=P^{-1}\left(\alpha\right)\left[\left[D_{x}f\left(x_{\epsilon}(\alpha),\alpha\right)-D_{x}f\left(x_{0}(\alpha),\alpha\right)\right]P\left(\alpha\right)\left(\begin{array}[]{c}u\\\
v\end{array}\right)\right.$ (C21)
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left.+\mathcal{O}\left(\left|u\right|^{2},\left|u\right|\left|v\right|,\left|v\right|^{2},\epsilon\left|u\right|,\epsilon\left|v\right|\right)\right]$
(C22)
$\displaystyle=P^{-1}\left(\alpha\right)\left[\mathcal{O}\left(\left|u\right|^{2},\left|u\right|\left|v\right|,\left|v\right|^{2},\epsilon\left|u\right|,\epsilon\left|v\right|\right)\right].$
(C23)

Under the assumptions (109) and (116), and by the definition of $\hat{f}$ in
(133), we have

$\hat{f}(0,0,0;\alpha)=0,\quad D_{u}\hat{f}(0,0,0;\alpha)=0,\quad
D_{v}\hat{f}(0,0,0;\alpha)=0,\quad D_{\epsilon}\hat{f}(0,0,0;\alpha)=0.$ (C24)

Specifically, in the $(u,v,\alpha)$ coordinates, the perturbed slow manifold
$\mathcal{L}_{\epsilon}$ satisfies

$\mathcal{L}_{\epsilon}=\left\\{(u,v,\alpha)\in\mathbb{R}^{n}\times\mathbb{R}\colon\,\,\,\,\,u=0,\quad
v=0\right\\}.$

We note the similarity between formulas (125)-(C24) and the setting of eqs.
(B32)-(B39) for general non-autonomous SSMs. Based on this similarity, we will
follow the same strategy here that we employed to compute invariant manifolds
in system (B32). Specifically, we will seek the perturbed invariant manifold
$\mathcal{M}_{\epsilon}$ in system (125) in the form of the asymptotic
expansion (137).

On the one hand, differentiating the definition of the invariant manifold
$\mathcal{M}_{\epsilon}$ from eq. (137) with respect to $t$ and using the ODE
(32), we obtain

$\displaystyle\dot{v}$
$\displaystyle=D_{u}h_{\epsilon}(u,\alpha)\dot{u}+\epsilon
D_{\alpha}h_{\epsilon}(u,\alpha)$
$\displaystyle=D_{u}h_{\epsilon}(u,\alpha)\left[A^{u}(\alpha)u+\hat{f}^{u}(u,h_{\epsilon}(u,\alpha),\alpha;\epsilon)\right]+\epsilon
D_{\alpha}h_{\epsilon}(u,\alpha).$ (C25)

On the other hand, substitution of the definition of the invariant manifold
$\mathcal{M}_{\epsilon}$ into (125) gives

$\displaystyle\dot{v}$
$\displaystyle=A^{v}h_{\epsilon}(u,\alpha)+\hat{f}^{v}(u,h_{\epsilon}(u,\alpha),\alpha;\epsilon).$
(C26)

Comparing (C25) and (C26) we obtain

$D_{u}h_{\epsilon}\left[A^{u}u+\hat{f^{u}}\right]+\epsilon
D_{\alpha}h_{\epsilon}=A^{v}h_{\epsilon}+\hat{f^{v}}.$ (C27)

####  C.3.1 Structure and solution of the invariance equation

Substitution of (137) into (C27) gives

$\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}h^{\mathbf{k}p}(\alpha)\epsilon^{p}D_{u}u^{\mathbf{k}}\left[A^{u}(\alpha)u+\hat{f}^{u}\left(u,\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}h^{\mathbf{k}p}(\alpha)u^{\mathbf{k}}\epsilon^{p},\alpha;\epsilon\right)\right]+\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}\left[h^{\mathbf{k}p}\right]^{\prime}(\alpha)u^{\mathbf{k}}\epsilon^{p+1}$

$=\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}A^{v}h^{\mathbf{k}p}(\alpha)u^{\mathbf{k}}\epsilon^{p}+\hat{f}^{v}\left(u,\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}h^{\mathbf{k}p}(\alpha)u^{\mathbf{k}}\epsilon^{p},\alpha;\epsilon\right).$
(C28)

As we did in the general non-autonomous case, we observe that

$A^{v}(\alpha)h^{\mathbf{k}p}(\alpha)u^{\mathbf{k}}-h^{\mathbf{k}p}(\alpha)D_{u}u^{\mathbf{k}}A^{u}(\alpha)u$
$=\left(\begin{array}[]{ccc}\lambda_{d+1}(\alpha)-\sum_{j=1}^{d}k_{j}\lambda_{j}(\alpha)&\cdots&0\\\
\vdots&\ddots&\vdots\\\
0&\cdots&\lambda_{n}(\alpha)-\sum_{j=1}^{d}k_{j}\lambda_{j}(\alpha)\end{array}\right)h^{\mathbf{k}p}(\alpha)u^{\mathbf{k}}.$
(C29)

Therefore, the invariance equation (C28) can be rewritten as

$\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}\left[h^{\mathbf{k}p}\right]^{\prime}(\alpha)u^{\mathbf{k}}\epsilon^{p+1}=\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}A_{\mathbf{k}}(\alpha)h^{\mathbf{k},p}(\alpha)u^{\mathbf{k}}\epsilon^{p}+\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}M^{\mathbf{k}p}(\alpha,h^{\mathbf{j}m})u^{\mathbf{k}}\epsilon^{p},$ (C30)

where

$\displaystyle A_{\mathbf{k}}$
$\displaystyle(\alpha)=\mathrm{diag}\left[\lambda_{\ell}(\alpha)-\sum_{j=1}^{d}k_{j}\lambda_{j}(\alpha)\right]_{\ell=d+1}^{n}\in\mathbb{C}^{\left(n-d\right)\times\left(n-d\right)},$
$\displaystyle\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}M^{\mathbf{k}p}(\alpha,h^{\mathbf{j}m})u^{\mathbf{k}}\epsilon^{p}$
$\displaystyle=M(u,\alpha,h^{\mathbf{j}m},\epsilon)=\hat{f}^{v}\left(u,\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}h^{\mathbf{k}p}(\alpha)u^{\mathbf{k}}\epsilon^{p},\alpha;\epsilon\right)$
(C31) $-\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}\epsilon^{p}\left(\begin{array}[]{ccc}h_{1}^{\mathbf{k}p}(\alpha)\frac{k_{1}u^{\mathbf{k}}}{u_{1}}&\cdots&h_{1}^{\mathbf{k}p}(\alpha)\frac{k_{d}u^{\mathbf{k}}}{u_{d}}\\\
\vdots&\ddots&\vdots\\\
h_{n-d}^{\mathbf{k}p}(\alpha)\frac{k_{1}u^{\mathbf{k}}}{u_{1}}&\cdots&h_{n-d}^{\mathbf{k}p}(\alpha)\frac{k_{d}u^{\mathbf{k}}}{u_{d}}\end{array}\right)\hat{f}^{u}\left(u,\sum_{\left|\left(\mathbf{k},p\right)\right|\geq
1}h^{\mathbf{k}p}(\alpha)u^{\mathbf{k}}\epsilon^{p},\alpha;\epsilon\right).$

Recall from eq. (C24) that $\hat{f}^{v}$ and $\hat{f}^{u}$ vanish for
$u,v,\epsilon=0$ and have no $\mathcal{O}\left(\left|v\right|\right)$ terms.
Then, as we concluded in the general non-autonomous case,
$M^{\mathbf{k}p}(\alpha,h^{\mathbf{j}m})$ will only depend on
$h^{\mathbf{j}m}$ that are lower in order, i.e.,

$D_{h^{\mathbf{j}m}}M^{\mathbf{k}p}(\alpha,h^{\mathbf{j}m})u^{\mathbf{k}}=0,\quad\left|(\mathbf{j},m)\right|\geq\left|(\mathbf{k},p)\right|.$
(C32)

Equating coefficients of equal powers of $u$ in eq. (C30), we then obtain the
system of equations

$\left[h^{\mathbf{k}(p-1)}\right]^{\prime}(\alpha)=A_{\mathbf{k}}(\alpha)h^{\mathbf{k}p}(\alpha)+M^{\mathbf{k}p}(\alpha,h^{\mathbf{j}m}),\quad\left|(\mathbf{j},m)\right|<\left|(\mathbf{k},p)\right|.$
(C33)

This is a recursively defined set of linear algebraic equations for
$h^{\mathbf{k}p}(\alpha)$, which can be uniquely solved as long as
$A_{\mathbf{k}}(\alpha)$ is nonsingular, i.e., the non-resonance conditions
(135) are satisfied. In that case, the recursive solution of (C33) starts from

$h^{\mathbf{0}p}(\alpha)\equiv 0,\quad p\geq 0,\qquad
h^{\mathbf{k}0}(\alpha)\equiv
h^{\mathbf{k}}(\alpha),\quad\mathbf{k}\in\mathbb{N}^{d},\qquad
h^{\mathbf{k}0}(\alpha)\equiv
h^{\mathbf{k}}(\alpha)=0,\quad\left|\mathbf{k}\right|=1,$ (C34)

and takes the form (139), with the quantities in eq. (140) obtained from
formulas (C31) using the relation (C32).

####  C.3.2 Reduced dynamics

To obtain the form of the reduced dynamics on the adiabatic SSM
$\mathcal{M}_{\epsilon}$, we consider the $u$ component of the transformation
formula (118),

__

$u=Q_{u}\left(\alpha\right)\left(x-x_{\epsilon}(\alpha)\right),$ (C35)

where the rows of $Q_{u}\left(\alpha\right)\in\mathbb{C}^{d\times n}$ are the
appropriately scaled unit left eigenvectors of $A(\alpha)$ corresponding to
its first $d$ (right) eigenvectors. This scaling is specified in statement
(ii) of the theorem, ensuring that the rows of $Q_{u}\left(\alpha\right)$
coincide with the first $d$ rows of $P^{-1}(\alpha)$.

Differentiation of (C35) with respect to time gives __

$\displaystyle\dot{u}$ $\displaystyle=\epsilon
Q_{u}^{\prime}\left(\alpha\right)\left(x-x_{\epsilon}(\alpha)\right)+Q_{u}\left(\alpha\right)\left(\dot{x}-\epsilon
x_{\epsilon}^{\prime}(\alpha)\right)$ $\displaystyle=\epsilon
Q_{u}^{\prime}\left(\alpha\right)P(\alpha)\left(\begin{array}[]{c}u\\\
v\end{array}\right)+Q_{u}\left(\alpha\right)\left(f(x,\alpha)-\epsilon
x_{\epsilon}^{\prime}(\alpha)\right)$
$\displaystyle=Q_{u}\left(\alpha\right)\left(f\left(x_{\epsilon}(\alpha)+P(\alpha)\left(\begin{array}[]{c}u\\\
v\end{array}\right),\alpha\right)-\epsilon
x_{\epsilon}^{\prime}(\alpha)\right)+\epsilon
Q_{u}^{\prime}\left(\alpha\right)P(\alpha)\left(\begin{array}[]{c}u\\\
v\end{array}\right).$

Restricting this last formula to the graph $v=h_{\epsilon}(u,\alpha)$ then
proves formula (147) in statement (ii) of the theorem.

## References

  * Alora et al. [2023a] J.I. Alora, M. Cenedese, E. Schmerling, G. Haller, and M. Pavone.  Data-driven spectral submanifold reduction for nonlinear optimal control of high-dimensional robots.  _Proc. IEEE Int. Conf. Robotics and Autom. (ICRA), London_ , pages 2627–2633, 2023a.
  * Alora et al. [2023b] J.I. Alora, M. Cenedese, E. Schmerling, G. Haller, and M. Pavone.  Practical deployment of spectral submanifold reduction for optimal control of high-dimensional systems.  _Proc. IFAC World Congress_ , page to appear, 2023b.
  * Axås et al. [2022] J. Axås, M. Cenedese, and G. Haller.  Fast data-driven model reduction for nonlinear dynamical systems.  (https://doi.org/10.1007/s11071-022-08014-0), 2022.
  * Breunung and Haller [2018] T. Breunung and G. Haller.  Explicit backbone curves from spectral submanifolds of forced-damped nonlinear mechanical systems.  _Proc. Royal Soc. A_ , 474:20180083, 2018.
  * Cabré et al. [2003] X. Cabré, E. Fontich, and R. de la Llave.  The parameterization method for invariant manifolds i: Manifolds associated to non-resonant subspaces.  _Indiana Univ. Math. J._ , 52(2):283–328, 2003\.
  * Carr [1982] J. Carr.  _Applications of Center Manifold Theory_.  Springer-Verlag, New York, 1982.
  * Cenedese et al. [2021] M. Cenedese, J. Axås, and G. Haller.  SSMLearn.  _https://github.com/haller-group/SSMLearn_ , 2021.
  * Cenedese et al. [2022a] M. Cenedese, J. Axås, B. Bäuerlein, K. Avila, and G. Haller.  Data-driven modeling and prediction of non-linearizable dynamics via spectral submanifolds.  _Nature Comm._ , 198, 2022a.
  * Cenedese et al. [2022b] M. Cenedese, J. Axås, H. Yang, M. Eriten, and G. Haller.  Data-driven nonlinear model reduction to spectral submanifolds in mechanical systems.  _Phil. Trans. R. Soc. A_ , 380:20210194, 2022b.
  * Constantine and Savits [1996] G. M. Constantine and T. H. Savits.  A Multivariate Faa di Bruno formula with applications.  _Trans. AMS_ , 348(2):503–520, 1996.
  * Cuong et al. [2019] L.V. Cuong, T.S. Doan, and S. Siegmund.  A Sternberg theorem for nonautonomous differential equations.  _J. Dyn. Diff. Equat._ , 31:1279–1299, 2019.
  * Eldering [2013] J. Eldering.  _Normally Hyperbolic Invariant Manifolds - The Noncompact Case_.  Atlantis Press, Paris, 2013.
  * Eldering et al. [2018] J. Eldering, M Kvalheim, and S. Revzen.  Global linearization and fiber bundle structure of invariant manifolds.  _Nonlinearity_ , 31:4202–4245, 2018.
  * Fenichel [1971] N. Fenichel.  Persistence and smoothness of invariant manifolds for flows.  _Indiana Univ. Math. J._ , 21(3):193–226, 1971\.
  * Fenichel [1979] N. Fenichel.  Geometric singular perturbation theory for ordinary differential equations.  _J. Diff. Eqs._ , 31(1):53–98, 1979.
  * Fontich et al. [2006] E. Fontich, R de la Llave, and P. MartÍn.  Invariant pre-foliations for non-resonant non-uniformly hyperbolic systems.  _Trans. AMS_ , 358(3):1317–1345, 2006.
  * Guckenheimer and Holmes [1983] J. Guckenheimer and P. Holmes.  _Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields_.  Springer, New York, 1983.
  * Haller and Ponsioen [2016] G. Haller and S. Ponsioen.  Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction.  _Nonlinear Dyn._ , 86(3):1493–1534, 2016.
  * Haller et al. [2023] G. Haller, B. Kaszás, A. Liu, and J. Axås.  Nonlinear model reduction to fractional and mixed-mode spectral submanifolds.  _Chaos_ , 33(6):063138, 2023.
  * Haro and de la Llave [2006] A. Haro and R. de la Llave.  A parameterization method for the computation of invariant tori and their whiskers in quasi-periodic maps: rigorous results.  _J. Differential Eqs._ , 228(2):530–579, 2006\.
  * Haro et al. [2016] A. Haro, M. Canadell, J.-L. Figueras, A. Luque, and J.M. Mondelo.  _The Parameterization Method for Invariant Manifolds: from Rigorous Results to Effective Computations_.  Springer, New York, 2016.
  * Hartman [1964] P. Hartman.  _Ordinary Differential Equations_.  John Wiley, New York, 1964.
  * Hirsch et al. [1977] M.W Hirsch, C.C. Pugh, and M. Shub.  _Invariant Manifolds. Lecture Notes in Mathematics, vol. 583_.  Springer-Verlag, Berlin-New York, 1977.
  * Jain and Haller [2022] S. Jain and G. Haller.  How to compute invariant manifolds and their reduced dynamics in high-dimensional finite-element models?  _Nonlin. Dyn._ , 107:1417–1450, 2022.
  * Jain et al. [2023] S. Jain, T. Thurner, M. Li, and G. Haller.  SSMTool 2.3: Computation of invariant manifolds and their reduced dynamics in high-dimensional mechanics problems.  _https://doi.org/10.5281/zenodo.4614201_ , pages 1417–1450, 2023\.
  * Kaszás et al. [2022] B. Kaszás, M. Cenedese, and G. Haller.  Dynamics-based machine learning of transitions in couette flow.  _Phys. Rev. Fluids._ , 7:L082402, 2022.
  * Kaundinya and Haller [2023] R.S. Kaundinya and G. Haller.  Matlab code for aperiodic SSM examples.  _https://github.com/haller-group/Aperiodic-SSMs.git_ , 2023.
  * Kerschen et al. [2009] G. Kerschen, M. Peeters, J.C. Golinval, and A.F. Vakakis.  Nonlinear normal modes, Part I: A useful framework for the structural dynamicist.  _Mech. Sys. Signal Proc._ , 23(1):170–194, 2009\.  ISSN 0888-3270.
  * Kloeden and Rassmussen [2011] P.E. Kloeden and M. Rassmussen.  _Nonautonomous Dynamical Systems_.  AMS, Providence, 2011.
  * Li et al. [2022a] M. Li, S. Jain, and G. Haller.  Nonlinear analysis of forced mechanical systems with internal resonance using spectral submanifolds, Part I: Periodic response and forced response curve.  _Nonlin.Dyn._ , pages https://doi.org/10.1007/s11071–022–07714–x, 2022a.
  * Li et al. [2022b] M. Li, S. Jain, and G. Haller.  Nonlinear analysis of forced mechanical systems with internal resonance using spectral submanifolds, Part II: Bifurcation and quasi-periodic response.  _Nonlin.Dyn._ , pages https://doi.org/10.1007/s11071–022–07476–6, 2022b.
  * Mikhlin and Avramov [2023] Y. Mikhlin and K.V. Avramov.  Nonlinear Normal Modes of Vibrating Mechanical Systems: 10 Years of Progress.  _Applied Mech. Rev._ , pages 1–57, 09 2023.
  * Mikhlin and Avramov [2011] Y. V. Mikhlin and K.V. Avramov.  Nonlinears Normal Modes for Vibrating Mechanical Systems. Review of Theoretical Developments.  _Applied Mech. Rev,_ , 63(6):060802, 06 2011.
  * Opreni et al. [2023] A. Opreni, A. Vizzaccaro, and A. Frangi.  High-order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to generic forcing terms and parametrically excited systems.  _Nonlin. Dyn._ , 111:5401–5447, 2023.
  * Palmer [1973] K.J. Palmer.  Generalization of Hartman’s linearization theorem.  _J. Math. Anal. App._ , 41:753–758, 1973.
  * Ponsioen et al. [2018] S. Ponsioen, T. Pedergnana, and G. Haller.  Automated computation of autonomous spectral submanifolds for nonlinear modal analysis.  _J. Sound Vibr._ , 420:269–295, 2018.
  * Ponsioen et al. [2019] S. Ponsioen, T. Pedergnana, and G. Haller.  Analytic prediction of isolated forced response curves from spectral submanifolds.  _Nonlin. Dyn._ , 98:2755–2773, 2019.
  * Ponsioen et al. [2020] S. Ponsioen, S. Jain, and G. Haller.  Model reduction to spectral submanifolds and forced-response calculation in high-dimensional mechanical systems.  _J. Sound Vibr._ , 488:115640, 2020.
  * Pötzsche [2010] C. Pötzsche.  _Geometric Theory of Discrete Nonautonomous Dynamical Systems_ , volume 2002.  Springer–Verlag, Berlin, Heidelberg, 01 2010.
  * Pötzsche and Rassmussen [2009] C. Pötzsche and M. Rassmussen.  Computation of nonautonomous invariant and inertial manifolds.  _Numer. Math._ , 112:449–483, 2009.
  * Pötzsche and Rassmussen [2016] C. Pötzsche and M. Rassmussen.  Taylor approximation of integral manifolds.  _J. Dyn. Diff. Equat._ , 18(2):427–461, 2016\.
  * Roberts [2015] A.J. Roberts.  _Model Emergent Dynamics in Complex Systems_.  SIAM, Providence, RI, 2015.
  * Sacker and Sell [1978] R.J. Sacker and G.R. Sell.  A spectral theory for linear differential systems.  _J. Diff. Eqs._ , 27:320–358, 1978.
  * Shaw and Pierre [1993] S.W. Shaw and C. Pierre.  Normal modes for non-linear vibratory systems.  _J. Sound Vibr._ , 164(1):85–124, 1993.
  * Son [2017] D. T. Son.  On analyticity for Lyapunov exponents of generic bounded linear random dynamical systems, 2017.
  * Thurner et al. [2023] T. Thurner, G. Haller, and S. Jain.  Nonautonomous spectral submanifolds for model reduction of nonlinear mechanical systems under parametric resonance.  _Chaos (submitted)_ , 2023.
  * Touzé et al. [2021] C. Touzé, A. Vizzaccaro, and O. Thomas.  Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques.  _Nonlinear Dyn._ , 105:1141–1190, 2021.
  * Vizzaccaro et al. [2022] A. Vizzaccaro, A. Opreni, L. Salles, A. Frangi, and C. Touzé.  High order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to large amplitude vibrations and uncovering of a folding point.  _Nonlin. Dyn._ , 110:525–571, 2022.
  * Vizzaccaro et al. [2023] A. Vizzaccaro, G. Goat, A. Frangi, and C. Touzé.  Direct parametrisation of invariant manifolds for generic non-autonomous systems including superharmonic resonances.  _arXiv:2306.09860_ , 2023.
  * Yomdin [1988] Y. Yomdin.  _Nonautonomous linearization. in Dynamical Systems (College Park, MD, 1986–87), Lecture Notes in Math., vol. 1342_ , volume 384, pages 718–726.  Springer, Berlin, 1988.