# Adjoint-Based Projections for Uncertainty Quantification near Stochastically
Perturbed Limit Cycles and Tori††thanks: A shortened treatment of the
covariance problem for limit cycles, including the first example in Section
2.2, was published previously in [1]. Additionally, the first author’s
doctoral dissertation (https://www.ideals.illinois.edu/items/124667) contains
an early draft of the content of this manuscript.

Zaid Ahsan Department of Mechanical Science and Engineering, University of
Illinois at Urbana-Champaign, Urbana, IL 61801 (email<EMAIL_ADDRESS>Harry Dankowicz22footnotemark: 2 Christian Kuehn Faculty of Mathematics,
Technical University of Munich, Boltzmannstr. 3, 85748 Garching bei München,
Germany (email<EMAIL_ADDRESS>
###### Abstract

This paper presents a new boundary-value problem formulation for quantifying
uncertainty induced by the presence of small Brownian noise near transversally
stable periodic orbits (limit cycles) and quasiperiodic invariant tori of the
deterministic dynamical systems obtained in the absence of noise. The
formulation uses adjoints to construct a continuous family of transversal
hyperplanes that are invariant under the linearized deterministic flow near
the limit cycle or quasiperiodic invariant torus. The intersections with each
hyperplane of stochastic trajectories that remain near the deterministic cycle
or torus over intermediate times may be approximated by a Gaussian
distribution whose covariance matrix can be obtained from the solution to the
corresponding boundary-value problem. In the case of limit cycles, the
analysis improves upon results in the literature through the explicit use of
state-space projections, transversality constraints, and symmetry-breaking
parameters that ensure uniqueness of the solution despite the lack of
hyperbolicity along the limit cycle. These same innovations are then
generalized to the case of a quasiperiodic invariant torus of arbitrary
dimension. In each case, a closed-form solution to the covariance boundary-
value problem is found in terms of a convergent series. The methodology is
validated against the results of numerical integration for two examples of
stochastically perturbed limit cycles and one example of a stochastically
perturbed two-dimensional quasiperiodic invariant torus in $\mathbb{R}^{2}$,
$\mathbb{R}^{2}\times S^{1}$, and $\mathbb{R}^{2}\times S^{1}$, respectively,
for which explicit expressions may be found for the associated covariance
functions using the proposed series solutions. Finally, an implementation of
the covariance boundary-value problem in the numerical continuation package
coco is applied to analyze the small-noise limit near a two-dimensional
quasiperiodic invariant torus in a nonlinear deterministic dynamical system in
$\mathbb{R}^{4}$ that does not support closed-form analysis. Excellent
agreement with numerical evidence from stochastic time integration shows the
potential for using deterministic continuation techniques to study the
influence of stochastic perturbations for both autonomous and periodically
excited deterministic vector fields.

##  1 Introduction

Asymptotically stable equilibria [3, 17, 35], limit cycles [9, 13, 31],
invariant curves [6], and transversally stable quasiperiodic invariant tori
[7, 8, 23] of deterministic dynamical systems continue to organize the local
flow over intermediate times also with the addition of noise of small
intensity. There results metastable dynamics near the deterministic limit sets
with stochastic trajectories characterized over such time scales by
approximately stationary Gaussian probability distributions of intersections
with hyperplanes transversal to the local vector field. Knowledge of these
distributions allows the investigator to predict the noise-induced spread of
uncertainty about the deterministic limit set [3, 10, 19, 20] and the
likelihood of transitions between different metastable objects [22, 26, 27,
38, 39, 41, 42].

Given a stochastic differential equation that appends noisy excitation to an
underlying deterministic dynamical system, numerical integration of solution
trajectories is a commonly deployed brute-force approach for the
characterization of noise-induced uncertainty, especially for global analysis
[25]. While computationally expensive and necessarily limited to particular
choices of parameter values and noise models, the results of such analysis
provide an important reality check for theoretical predictions [10]. For local
analysis—of primary interest here—a more affordable alternative is the direct
computation of covariance ellipsoids [4, 37] in each hyperplane through the
solution of an associated Lyapunov equation [11, 12, 24] for a covariance
tensor function. Indeed, as shown by Kuehn [28] (see also [2, 29, 30]) in the
case of equilibria, such a computation may be embedded within a numerical
continuation framework in order to investigate the parameter dependence of the
metastable dynamics and the probabilities of noise-induced transitions between
metastable objects (see, e.g., [5]). One goal of this paper is to extend such
a methodology, compatible with a numerical parameter continuation framework,
also to the case of limit cycles and transversally stable quasiperiodic
invariant tori.

To achieve such an outcome, we expand on our preliminary treatment of limit
cycles in [1] to formulate a smooth, regular boundary-value problem, such that
pairings of parameter values, deterministic limit sets, and associated
covariance tensor functions span a smooth manifold in a suitably defined
variable space. Such a formulation then supports simultaneous continuation of
limit sets and covariance functions using well-established software packages,
e.g., auto [16] and coco [14], including the possibility of monitoring or
constraining the geometry of the covariance ellipsoids during parameter
variations. In this paper, we rely on an implementation in coco to enable
numerical validation of the theoretical predictions, and make the
corresponding code freely available at
https://github.com/hdankowicz/covariance-bvp2023-scripts. Although we leave
the utility of such a combined framework largely to future work, a natural
application is to problems in optimal robust design of noise-disturbed
engineering systems.

Several features distinguish our treatment from that in most of the existing
literature. Foremost among these is the construction of transversal
hyperplanes using projection operators defined in terms of the solutions to
appropriately defined adjoint boundary-value problems with integral
constraints as described in [15]. Notably, this construction is invariant
under local coordinate transformations. The resultant foliation of hyperplanes
is also invariant under the linearized flow near the limit set, which
contracts the corresponding transversal deviations exponentially fast to $0$.
Neither of these properties hold for the typical construction (for autonomous
systems) of transversal hyperplanes that are chosen arbitrarily [34, 43] or
required to be orthogonal to the local vector field [7, 33], as this
construction is not invariant under local coordinate transformations, the
resultant hyperplanes are not invariant under the linearized flow, and
deviations tangential to the limit set must be eliminated, for example, by
projection onto a local set of coordinates in each hyperplane (as in [33] and
[43]; see also [21] for an analysis of phase diffusion along a limit cycle).
Such a coordinate projection is not required by the proposed construction
which proceeds to make predictions in the original coordinates. Furthermore,
as is demonstrated here through several examples, for non-autonomous,
periodically excited dynamical systems, the proposed construction
automatically generates hyperplanes in stroboscopic sections (the _a priori_
choice in [10]). Finally, in our computational implementation, we benefit from
the existence of automated support for the construction of the adjoint
boundary-value problems using the staged construction paradigm in coco [32].

As a consequence of these choices, we obtain a Lyapunov equation that
conserves the scalar-valued image of the covariance tensor function on any
pair of adjoint vectors. We restrict attention to the intersection of the
zero-level sets of all such images by appending a corresponding set of
boundary conditions and augmenting the Lyapunov equation with a regularizing
damping term that must vanish on solutions to the full covariance boundary-
value problem (see [18, 36] for a general methodology). Using the invariance
of the foliation under the linearized flow and exponential contraction, we
obtain a unique solution for the covariance tensor function in terms of an
exponentially convergent series. The result is a function whose value at each
point of the limit set is a positive semi-definite matrix with as many nonzero
eigenvalues, generically, as the co-dimension of the limit set.

That the analysis for limit cycles may be extended also to the case of
quasiperiodic invariant tori is demonstrated already in [7, 23, 43]. In our
formulation, this extension results in several coupled partial differential
equations on a cylindrical domain with characteristics parallel to the time
axis and with boundary conditions, integral constraints, and regularizing
terms enforcing quasiperiodicity, non-degeneracy, and uniqueness. In our
numerical implementation, we rely on a problem discretization in terms of
Fourier series in the non-temporal variables and continuous, piecewise-
polynomial functions along characteristics, and show how the discretized
adjoint integro-differential boundary-value problems may be obtained from a
variational formulation, as also implemented in coco.

The remainder of this paper is arranged as follows. Section 2 considers the
influence of noise on the dynamics near limit cycles. Rigorous derivations for
systems of arbitrary dimension are there validated against the predictions of
numerically integrated stochastic trajectories for two examples where the
sought covariance functions may be obtained in closed form, viz., an
autonomous, deterministic vector field in $\mathbb{R}^{2}$ obtained from the
Hopf normal form and a harmonically excited, damped, linear oscillator with
dynamics in $\mathbb{R}^{2}\times\mathbb{S}^{1}$. Section 3 generalizes the
construction of families of local projections to the case of transversally
stable quasiperiodic invariant tori of arbitrary dimension in problems of
arbitrary dimension and shows that a unique covariance matrix may again be
associated with each corresponding hyperplane. Here, validation is delivered
first using a two-dimensional torus of a deterministic vector field in
$\mathbb{R}^{2}\times\mathbb{S}^{1}$, for which closed-form analysis of the
covariance function is possible, and then a two-dimensional torus of a
deterministic vector field in $\mathbb{R}^{4}$, which is not amenable to such
analysis and for which an approximation must be obtained from a discretization
of the governing boundary-value problem using the software package coco. This
latter implementation also illustrates how uncertainty quantification may be
appended to or integrated with parameter continuation of limit cycles or
quasiperiodic invariant tori. Appendix A collects important results that
enable automated problem construction and approximate analysis in coco of the
various adjoint boundary-value problems, including a choice of discretization
in terms of finite sets of characteristics in the case of quasiperiodic
invariant tori that can be used analogously for finding the tori in the first
place and analyzing the associated covariance functions. A concluding
discussion in Section 4 reflects on directions for further study.

##  2 Periodic Orbits

In this section, we derive the set of adjoint and covariance boundary-value
problems (previously stated without derivation in [1]) whose combined
solution, for a given limit cycle of an autonomous dynamical system,
characterizes the metastable dynamics in the limit of small noise intensity.
The theoretical predictions are validated against two examples for which
closed-form expressions may be compared to the results of stochastic numerical
integration.

###  2.1 Theoretical derivations

Let $\gamma(t)\in\mathbb{R}^{n}$ be a periodic solution with period $1$ of the
smooth dynamical system defined by the ordinary differential equation (ODE)

$\frac{\textnormal{d}x}{\textnormal{d}t}=:\dot{x}=Tf(x),\quad\mbox{ where
}x:=x(t)\in\mathbb{R}^{n},$ (1)

for some positive scalar111We obtain the form of (1) from a dynamical system
of the form $y^{\prime}(\eta)=f(y(\eta))$ with periodic orbit $\gamma(\eta/T)$
of period $T$ by letting $t=\eta/T$ and $x(t)=y(Tt)$. For an autonomous vector
field, $T$ is not typically known a priori and is, instead, solved for
simultaneously with $\gamma$ by appending a suitable phase condition to the
periodic boundary-value problem in $x$, e.g., $h(x(0))=0$ for some smooth
scalar-valued function $h$. Throughout, we omit explicit dependence on problem
parameters to reduce the notational complexity. $T$. Let $X(t)$ denote the
solution to the corresponding variational initial-value problem

$\dot{X}=T\mathrm{D}f(\gamma)X,\quad X(0)=I_{n},$ (2)

where $\mathrm{D}f$ denotes the Jacobian of $f$ and $I_{n}$ is the identity
matrix of size $n\times n$. By periodicity of $\gamma$, it follows that
$X(t+1)=X(t)X(1)$. Furthermore, $X^{-1}(t)f(\gamma(t))\equiv f(\gamma(0))$,
since equality holds at $t=0$ and differentiation shows that the left-hand
side is constant. As a special case, again by periodicity of $\gamma$, it
follows that $f(\gamma(0))$ is a right nullvector of $X(1)-I_{n}$ and
$f(\gamma(t))$ is a right nullvector of the linear operator

$\Gamma(t):=X(t)\left(X(1)-I_{n}\right)X^{-1}(t).$ (3)

Following [15], assume next that $\gamma$ is hyperbolic in the sense that
there exists a $1$-periodic, smooth family of projections of the form

$Q(t):=X(t)(I_{n}-f(\gamma(0))w^{\mathsf{T}})X^{-1}(t)$ (4)

for some vector $w$, such that $\Gamma(t)$ is a bijection on the image of
$Q(t)$. As shown below, this is equivalent to the standard definition of
hyperbolicity in terms of the eigenvalues of the monodromy matrix $X(1)$.
Importantly, the family $Q(t)$ and the complementary family $I_{n}-Q(t)$
respect the spectral decomposition of $\Gamma(t)+I_{n}$.

####  2.1.1 The adjoint boundary-value problem

The following sequence of lemmas and corollaries culminate in the formulation
of the adjoint boundary-value problem in Lemma 2.5 and a projection of the
local dynamics onto transversal directions relative to the limit cycle in
Corollary 2.10. The discussion is adapted from [15] and organized to set up a
point of departure for the subsequent analysis of noise-induced transversal
dynamics.

######  Lemma 2.1.

The vector $w$ satisfies the normalization condition

$w^{\mathsf{T}}f(\gamma(0))=1.$ (5)

###### Proof.

Since $Q(0)$ is a projection, it is idempotent, i.e., $Q^{2}(0)=Q(0)$ or, by
(4),

$I_{n}-f(\gamma(0))w^{\mathsf{T}}=\left(I_{n}-f(\gamma(0))w^{\mathsf{T}}\right)^{2}.$
(6)

Algebraic manipulation shows that this holds if and only if (5) is satisfied.
∎

######  Lemma 2.2.

Let $\operatorname{Im}Q(t)$ and $\ker Q(t)$ denote the image and kernel of the
projection $Q(t)$. Then,

$\operatorname{Im}Q(t)=\\{X^{{-1},\mathsf{T}}(t)w\\}^{\perp}\mbox{ and }\ker
Q(t)=\mathrm{span}\\{X(t)f(\gamma(0))\\}.$ (7)

###### Proof.

By (4) it suffices to show this for $t=0$. If $v\in\ker Q(0)$, then

$0=\left(I_{n}-f(\gamma(0))w^{\mathsf{T}}\right)v=v-f(\gamma(0))w^{\mathsf{T}}v,$
(8)

i.e., that $v\in\mathrm{span}\\{f(\gamma(0))\\}$. That $\ker
Q(0)=\mathrm{span}\\{f(\gamma(0))\\}$ then follows from Lemma 2.1. Similarly,
if $v\in\operatorname{Im}Q(0)$, then there exists a vector $u$ such that

$v=(I_{n}-f(\gamma(0))w^{\mathsf{T}})u$ (9)

and, by Lemma 2.1,

$w^{\mathsf{T}}v=\left(w^{\mathsf{T}}-w^{\mathsf{T}}f(\gamma(0))w^{\mathsf{T}}\right)u=0,$
(10)

i.e., $\operatorname{Im}Q(0)\in\\{w\\}^{\perp}$. Equality then follows from
the rank-nullity theorem. ∎

######  Lemma 2.3.

The vector $\lambda(t):=X^{{-1},\mathsf{T}}(t)w$ is the unique left nullvector
of $\Gamma(t)$ for which

$\lambda^{\mathsf{T}}(t)f(\gamma(t))=1.$ (11)

###### Proof.

By (3) it suffices to show this for $t=0$, where (11) reduces to (5). To show
that $w$ is a left nullvector of $\Gamma(0)=X(1)-I_{n}$, note that the left
nullvector property $w^{\mathsf{T}}\Gamma(0)=0$ holds if and only if
$w^{\mathsf{T}}X(1)=w^{\mathsf{T}}$. But the latter follows by $1$-periodicity
of $Q$, since $Q(0)=Q(1)$ implies that

$X(1)(I_{n}-f(\gamma(0))w^{\mathsf{T}})X^{-1}(1)=I_{n}-f(\gamma(0))w^{\mathsf{T}},$
(12)

which reduces to $f(\gamma(0))w^{\top}=f(\gamma(0))w^{\top}X(1)$ by direct
algebraic manipulation and use of the fact that
$X(1)f(\gamma(0))=f(\gamma(0))$. Finally, suppose that $\tilde{w}$ is any left
nullvector of $\Gamma(0)$ and let $v$ denote an arbitrary vector in
$\operatorname{Im}Q(0)$. Since $\Gamma(0)$ is onto $\operatorname{Im}Q(0)$,
there exists a $u$ such that $\Gamma(0)u=v$ and, consequently,
$\tilde{w}^{\mathsf{T}}v=\tilde{w}^{\mathsf{T}}\Gamma(0)u=0$. By Lemma 2.2,
$\tilde{w}\in\mathrm{span}\\{w\\}$ and the uniqueness claim follows as well. ∎

######  Corollary 2.4.

The right nullspace of $\Gamma(t)$ consists of vectors tangential to the curve
$s\mapsto\gamma(s)$ at $s=t$.

###### Proof.

By (3) it suffices to show this for $t=0$. By substitution from (3) and (4),
it follows that

$\Gamma(0)Q(0)=\left(X(1)-I_{n}\right)\left(I_{n}-f(\gamma(0))w^{\mathsf{T}}\right)=\Gamma(0),$
(13)

where we again used the fact that $X(1)f(\gamma(0))=f(\gamma(0))$. Now suppose
that $v$ is a right nullvector of $\Gamma(0)$. It follows that
$\Gamma(0)Q(0)v=\Gamma(0)v=0$. Since $\Gamma(0)$ is one-to-one on the image of
$Q(0)$, $Q(0)v=0$ and, by Lemma 2.2, $v\in\mathrm{span}\\{f(\gamma(0))\\}$. ∎

Conversely, suppose that the eigenspace of the monodromy matrix $X(1)$
corresponding to the eigenvalue $1$ is given by
$\mathrm{span}\\{f(\gamma(0))\\}$ and choose $w$ in the construction of $Q(t)$
as the unique left nullvector of $\Gamma(0)$ such that (5) holds. Then, by
(13), $v\in\operatorname{Im}Q(0)$ and $\Gamma(0)v=0$ imply that $v=0$.
Similarly, since

$Q(0)\Gamma(0)=\left(I_{n}-f(\gamma(0))w^{\mathsf{T}}\right)\left(X(1)-I_{n}\right)=\Gamma(0),$
(14)

it follows from the rank-nullity theorem that if $v\in\operatorname{Im}Q(0)$
then there exists a $\tilde{v}\in\operatorname{Im}Q(0)$ such that
$\Gamma(0)\tilde{v}=v$. We conclude that $\Gamma(0)$ is a bijection on the
image of $Q(0)$ and, consequently, that $\Gamma(t)$ is a bijection on the
image of $Q(t)$.

######  Lemma 2.5.

The function $\lambda(t):=X^{{-1},\mathsf{T}}(t)w$ is the unique solution to
the adjoint boundary-value problem

$0=-\dot{\mu}^{\mathsf{T}}-T\mu^{\mathsf{T}}\mathrm{D}f(\gamma),\quad
0=\mu^{\mathsf{T}}(0)-\mu^{\mathsf{T}}(1),\quad
1=\int_{0}^{1}\mu^{\mathsf{T}}f(\gamma(s))\,\mathrm{d}t.$ (15)

###### Proof.

Let $\mu$ denote an arbitrary solution to this boundary-value problem. It
follows that $\mu^{\mathsf{T}}(t)X(t)\equiv w^{\mathsf{T}}$ and
$\mu^{\mathsf{T}}(t)f(\gamma(s))\equiv 1$, since i) differentiation and use of
(2) and (15) shows that the left-hand sides are constant, ii) the integral
condition implies that the second constant equals $1$ and, consequently,
$\mu^{\mathsf{T}}(0)f(\gamma(0))=1$, iii) periodicity shows that
$\mu^{\mathsf{T}}(0)X(1)=\mu^{\mathsf{T}}(0)$ and, by the proof of Lemma 2.3,
that $\mu^{\mathsf{T}}(0)$ is a left nullvector of $\Gamma(0)$ that iv) by
Lemma 2.3 must equal $w$. ∎

######  Corollary 2.6.

In the notation of the previous lemma,

$Q(t)=I_{n}-f(\gamma(s))\lambda^{\mathsf{T}}(t)$ (16)

and projects along $f(\gamma(s))$ onto a co-dimension-one hyperplane
orthogonal to $\lambda(t)$.

######  Lemma 2.7.

The family of operators $P(t)=Q(t)$ is the unique periodic family of
projections onto co-dimension-one hyperplanes which satisfies the invariance
condition

$X(t)P(0)=P(t)X(t)$ (17)

and such that $\Gamma(t)$ is a bijection on the image of $P(t)$.

###### Proof.

For any such family of projections, there exist two periodic functions $\phi$
and $\mu$, such that $P(t)=I_{n}-\phi(t)\mu^{\mathsf{T}}(t)$,
$\mu^{\mathsf{T}}(t)\phi(t)\equiv 1$ and, by invariance, $\phi(t)=X(t)\phi(0)$
and $\mu(t)=X^{-1,\mathsf{T}}(t)\mu(0)$. By direct algebraic manipulation and
periodicity, we obtain

$\Gamma(0)P(0)=\Gamma(0)=P(0)\Gamma(0).$ (18)

It follows from the first equality that $P(0)f(\gamma(0))$ is a nullvector of
$\Gamma(0)$ and, by injectivity of $\Gamma(0)$ on $\operatorname{Im}P(0)$,
that $P(0)f(\gamma(0))=0$. Without loss of generality, we conclude that
$\phi(0)=f(\gamma(0))$. From the second equality, we obtain
$\mu^{T}(0)\Gamma(0)=0$. By surjectivity of $\Gamma(0)$ on
$\operatorname{Im}P(0)$, as in the proof of Lemma 2.3, then $\mu(0)=w$. ∎

######  Corollary 2.8.

The stable and unstable eigenspaces of $\Gamma(t)+I_{n}$ lie in the image of
$Q(t)$.

We proceed to investigate the local geometry on a neighborhood of $\gamma$
characterized in terms of the projections $Q$. These provide a local
decomposition into state-space displacements tangential and transversal to the
periodic orbit. When convenient, we perform formal calculations with
differentials in order to capture leading-order effects.

######  Lemma 2.9.

For $x$ sufficiently close to $\gamma(t)$ for some $t$, there exists a unique
time $\tau$ (with $\tau\approx t$) such that

$x=\gamma(\tau)+x_{\mathrm{tr}},$ (19)

where $\lambda(\tau)^{\mathsf{T}}x_{\mathrm{tr}}=0$ and, consequently,
$Q(\tau)x_{\mathrm{tr}}=x_{\mathrm{tr}}$.

###### Proof.

Let $H(x,\tau)=\lambda^{\mathsf{T}}(\tau)\left(x-\gamma(\tau)\right)$. It
follows that $H(\gamma(t),t)=0$ and

$\partial_{\tau}H(\gamma(t),t)=-T\lambda^{\mathsf{T}}(t)f(\gamma(s))=-T.$ (20)

The claim follows from the implicit function theorem. ∎

######  Corollary 2.10.

For the differential $\mathrm{d}x$, the result of Lemma 2.9 implies the
existence of differentials $\mathrm{d}\tau$ and $\mathrm{d}x_{\mathrm{tr}}$,
such that

$x+\mathrm{d}x=\gamma(\tau+\mathrm{d}\tau)+x_{\mathrm{tr}}+\mathrm{d}x_{\mathrm{tr}},$
(21)

where
$\lambda^{\mathsf{T}}(\tau+\mathrm{d}\tau)(x_{\mathrm{tr}}+\mathrm{d}x_{\mathrm{tr}})=0$
and, consequently,

$Q(\tau+\mathrm{d}\tau)(x_{\mathrm{tr}}+\mathrm{d}x_{\mathrm{tr}})=x_{\mathrm{tr}}+\mathrm{d}x_{\mathrm{tr}}.$
(22)

It follows that

$\mathrm{d}x_{\mathrm{tr}}-T\mathrm{D}f(\gamma(\tau))x_{\mathrm{tr}}\mathrm{d}\tau=Q(\tau)\left(\mathrm{d}x-T\mathrm{D}f(\gamma(\tau))x_{\mathrm{tr}}\mathrm{d}\tau\right).$
(23)

###### Proof.

Eq. (21) implies that

$\mathrm{d}x_{\mathrm{tr}}=\mathrm{d}x-Tf(\gamma(\tau))\mathrm{d}\tau.$ (24)

Similarly, from (22), we obtain

$\mathrm{d}x_{\mathrm{tr}}=\dot{Q}(\tau)x_{\mathrm{tr}}\mathrm{d}\tau+Q(\tau)\mathrm{d}x_{\mathrm{tr}}.$
(25)

From the definition of $Q$ it follows that

$\dot{Q}(\tau)=-T\mathrm{D}f(\gamma(\tau))f(\gamma(\tau))\lambda^{\mathsf{T}}(\tau)+Tf(\gamma(\tau))\lambda^{\mathsf{T}}(\tau)\mathrm{D}f(\gamma(\tau))$
(26)

and, consequently, that

$\dot{Q}(\tau)x_{\mathrm{tr}}=T(I_{n}-Q(\tau))\mathrm{D}f(\gamma(\tau))x_{\mathrm{tr}},$
(27)

since $\lambda^{\mathsf{T}}(\tau)x_{\mathrm{tr}}=0$. Substitution of (24) and
(27) into the right-hand side of (25) yields (23), since
$Q(\tau)f(\gamma(\tau))=0$. ∎

####  2.1.2 The covariance boundary-value problem

In lieu of the deterministic dynamics in (1), now consider the Itô SDE

$\mathrm{d}x=Tf(x)\mathrm{d}t+\sigma\sqrt{T}F(x)\mathrm{d}W_{t}$ (28)

in terms of the noise intensity $\sigma$ and a vector $W_{t}\in\mathbb{R}^{m}$
of independent standard Brownian motions. We consider the limit of small noise
intensity and are concerned only with leading-order expansions in $\sigma$. In
the notation of the preceding lemmas, we continue to calculate with
differentials. Since we are only interested in first-order terms in $\sigma$,
we drop any Itô corrections arising from the Itô-formula. The following
proposition is then relatively straightforward.

######  Proposition 2.11.

For $x$ sufficiently close to $\gamma(t)$ for some $t$, to first order in
$\|x_{\mathrm{tr}}\|$ and $\sigma$,

$\mathrm{d}x_{\mathrm{tr}}=T\mathrm{D}f(\gamma(\tau))x_{\mathrm{tr}}\mathrm{d}\tau+\sigma\sqrt{T}Q(\tau)F(\gamma(\tau))\mathrm{d}W_{\tau}.$
(29)

###### Proof.

For $x_{\mathrm{tr}}=0$ and $\sigma=0$, we have $\tau=t$ and
$\mathrm{d}x_{\mathrm{tr}}=0$. In this limit, (28) and (24) imply that

$0=Tf\left(\gamma(\tau)\right)(\mathrm{d}t-\mathrm{d}\tau),$ (30)

from which it follows that $\mathrm{d}\tau=\mathrm{d}t$ to zeroth order in
$\|x_{\mathrm{tr}}\|$ and $\sigma$. To first order in $\|x_{\mathrm{tr}}\|$
and $\sigma$, (28) then yields

$Q(\tau)\left(\mathrm{d}x-T\mathrm{D}f(\gamma(\tau))x_{\mathrm{tr}}\mathrm{d}\tau\right)=\sigma\sqrt{T}Q(\tau)F(\gamma(\tau))\mathrm{d}W_{\tau}$
(31)

and substitution in (23) yields the claimed result. ∎

The Ornstein-Uhlenbeck process in (29) is a leading-order approximation of the
local dynamics near $\gamma$ that neglects the cumulative effects of higher-
order terms in $\|x_{\mathrm{tr}}\|$ and $\sigma$. Given an initial condition
$x_{\mathrm{tr}}(0)$, (29) may be solved explicitly to yield

$\displaystyle x_{\mathrm{tr}}(\tau)$
$\displaystyle=X(\tau)x_{\mathrm{tr}}(0)+\sigma\sqrt{T}X(\tau)\int_{0}^{\tau}X^{-1}(s)Q(s)F(\gamma(s)))\mathrm{d}W_{s}$
$\displaystyle=X(\tau)x_{\mathrm{tr}}(0)+\sigma\sqrt{T}Q(\tau)X(\tau)\int_{0}^{\tau}G(s)\mathrm{d}W_{s},$
(32)

where $G(s)=X^{-1}(s)F(\gamma(s))$ and the invariance condition (17) with
$P(t)=Q(t)$ was used twice. It follows that
$\lambda^{\mathsf{T}}x_{\mathrm{tr}}\equiv 0$ provided that
$w^{\mathsf{T}}x_{\mathrm{tr}}(0)=0$. The rescaled covariance matrix

$C:=\frac{1}{\sigma^{2}}\mathbb{E}\left[x_{\mathrm{tr}}x_{\mathrm{tr}}^{\mathrm{T}}\right]$
(33)

is then given by

$\displaystyle C(\tau)$
$\displaystyle=X(\tau)C(0)X^{\mathsf{T}}(\tau)+TQ(\tau)X(\tau)\left(\int_{0}^{\tau}G(s)G^{\mathsf{T}}(s)\mathrm{d}s\right)X^{\mathsf{T}}(\tau)Q^{\mathsf{T}}(\tau),$
(34)

where we used the stochastic equalities [25]

$\mathbb{E}\left[\int_{0}^{\tau}Y(s)\mathrm{d}W_{s}\right]=0$ (35)

and

$\mathbb{E}\left[\int_{0}^{\tau}Y(s)\mathrm{d}W_{s}\int_{0}^{\tau}Z^{\mathsf{T}}(s)\mathrm{d}W_{s}\right]=\int_{0}^{\tau}Y(s)Z^{\mathsf{T}}(s)\mathrm{d}s.$
(36)

It follows that $\lambda^{\mathsf{T}}C\lambda\equiv w^{\mathsf{T}}C(0)w$ and
that $\lambda^{\mathsf{T}}C\lambda\equiv 0$ if and only if
$w^{\mathsf{T}}C(0)w=0$.

######  Lemma 2.12.

Let

$\mathcal{I}=\int_{0}^{1}G(s)G^{\mathsf{T}}(s)\mathrm{d}s.$ (37)

Then, for arbitrary $C(0)$ and $k\in\mathbb{N}\cup\\{0\\}$,

$C(k+1)=X(1)C(k)X^{\mathsf{T}}(1)+TQ(0)X(1)\mathcal{I}X^{\mathsf{T}}(1)Q^{\mathsf{T}}(0).$
(38)

###### Proof.

By (34),

$C(k)=X(k)C(0)X^{\mathsf{T}}(k)+TQ(0)X(k)\left(\int_{0}^{k}G(s)G^{\mathsf{T}}(s)\mathrm{d}s\right)X^{\mathsf{T}}(k)Q^{\mathsf{T}}(0).$
(39)

Moreover, by the definition of $G$, we have

$\mathcal{I}=X(k)\left(\int_{k}^{k+1}G(s)G^{\mathsf{T}}(s)\mathrm{d}s\right)X^{\mathsf{T}}(k).$
(40)

The claim follows by substitution. ∎

We proceed to assume that the unstable eigenspace of $X(1)$ is empty, such
that

$X(k)=f(\gamma(0))w^{\mathsf{T}}+\mathcal{O}\left(\mathrm{exp}(-k/\tau_{\mathrm{tr}})\right)$
(41)

for some positive constant $\tau_{\mathrm{tr}}$ and large $k$ (here, the
asymptotics is with respect to large times so that the
$\mathcal{O}(\cdot)$-term becomes negligible as k becomes large).

######  Proposition 2.13.

Given a non-negative constant $b$, there exists an initial condition $C(0)$
with $w^{\mathsf{T}}C(0)w=b$, such that $C(k)=C(0)$ for all $k$.

###### Proof.

The previous lemma shows that $C(1)=C(0)\Rightarrow C(k)=C(0)$ for all $k$.
Suppose that

$\displaystyle C(0)$
$\displaystyle=bf(\gamma(0))f^{\mathsf{T}}(\gamma(0))+T\sum_{n=1}^{\infty}Q(0)X(n)\mathcal{I}X^{\mathsf{T}}(n)Q^{\mathsf{T}}(0).$
(42)

The series converges by the assumption on $X(1)$. By substitution and use of
periodicity and the invariance condition

$X(1)C(0)X^{\mathsf{T}}(1)=T\sum_{n=2}^{\infty}Q(0)X(n)\mathcal{I}X^{\mathsf{T}}(n)Q^{\mathsf{T}}(0)$
(43)

and the claim follows by evaluation of (34) at $\tau=1$. ∎

######  Corollary 2.14.

The function $C(\tau)$ in (34) with $w^{\mathsf{T}}C(0)w=b$ converges to the
periodic function obtained with $C(0)$ given by (42) as
$\tau\rightarrow\infty$.

######  Proposition 2.15.

Let $C_{\mathrm{per}}$ denote the function in (34) obtained with $C(0)$ given
by (42) when $b=0$. It follows that, when $a=0$, $C_{\mathrm{per}}$ is the
unique periodic solution of the periodic Lyapunov equation (cf. [12])

$\displaystyle\dot{C}$
$\displaystyle=T\left(\mathrm{D}f(\gamma)C+C\mathrm{D}f^{\mathsf{T}}(\gamma)+QF(\gamma)F^{\mathsf{T}}(\gamma)Q^{\mathsf{T}}+af(\gamma)f^{\mathsf{T}}(\gamma)\right),$
(44)

for which $w^{\mathsf{T}}C(0)w=0$ and that no periodic solution exists
otherwise.

###### Proof.

It follows from (44) that

$\frac{d}{dt}\lambda^{\mathsf{T}}C\lambda=aT$ (45)

Thus, $C$ is periodic only if $a=0$. The transformation
$C=X\tilde{C}X^{\mathsf{T}}$ then implies that

$\dot{\tilde{C}}=TQ(0)GG^{\mathsf{T}}Q^{\mathsf{T}}(0)$ (46)

and the conclusion follows by integration. ∎

######  Corollary 2.16.

Let $\mathfrak{vec}$ denote the vectorization operator that turns a matrix
into a column vector. Then, $\mathfrak{vec}\left(C_{\mathrm{per}}(0)\right)$
may be obtained from the first $n^{2}$ components of the image of the column
vector

$\displaystyle\begin{pmatrix}\mathfrak{vec}\left((Q(0)X(1)\mathcal{I}X^{\mathsf{T}}(1)Q^{\mathsf{T}}(0)\right)\\\
0\end{pmatrix}$ (47)

under multiplication by the matrix

$\begin{pmatrix}I_{2n}-X(1)\otimes X(1)&f(\gamma(0))\otimes f(\gamma(0))\\\
w^{\mathsf{T}}\otimes w^{\mathsf{T}}&0\end{pmatrix}^{-1}\\\ $

where the inverse exists since no pair of eigenvalues of $X(1)$ are reciprocal
numbers.

###  2.2 Numerical examples

An analytically tractable example (a version of which was considered in [1])
is given by the periodic solution

$\gamma(t)=\begin{pmatrix}\cos 2\pi t\\\ \sin 2\pi t\end{pmatrix}$ (48)

with period $1$ of the dynamical system (1) with

$f(x)=\begin{pmatrix}\vspace{1mm}x_{1}-x_{2}-x_{1}(x_{1}^{2}+x_{2}^{2})\\\
x_{1}+x_{2}-x_{2}(x_{1}^{2}+x_{2}^{2})\end{pmatrix}$ (49)

and $T=2\pi$, such that

$f(\gamma(t))=\begin{pmatrix}-\sin 2\pi t\\\ \cos 2\pi t\end{pmatrix}.$ (50)

Since

$X(t)=\begin{pmatrix}e^{-4\pi t}\cos 2\pi t&-\sin 2\pi t\\\ e^{-4\pi t}\sin
2\pi t&\cos 2\pi t\end{pmatrix},$ (51)

the eigenvalues of the monodromy matrix $X(1)$ equal $1$ and $e^{-4\pi}$ with
eigenspaces spanned by $f(\gamma(0))=(0,1)^{\mathsf{T}}$ and
$(1,0)^{\mathsf{T}}$, respectively. Then, $w^{\mathsf{T}}=(0,1)$ and, by
definition,

$\lambda(t)=X^{-1,\mathsf{T}}(t)w=\begin{pmatrix}-\sin 2\pi t\\\ \cos 2\pi
t\end{pmatrix}.$ (52)

and

$Q(t)=I_{2}-f(\gamma(t))\lambda^{\mathsf{T}}(t)=\begin{pmatrix}\cos^{2}2\pi
t&\cos 2\pi t\sin 2\pi t\\\ \cos 2\pi t\sin 2\pi t&\sin^{2}2\pi
t\end{pmatrix}.$ (53)

It follows that

$TQ(0)\left(\int_{0}^{t}G(s)G^{\mathsf{T}}(s)\mathrm{d}s\right)Q^{\mathsf{T}}(0)=\begin{pmatrix}\mathcal{I}(t)&0\\\
0&0\end{pmatrix},$ (54)

where $\mathcal{I}(0)=0$. The result of Proposition 2.13 then implies that
every periodic solution of (34) is of the form

$\displaystyle C(t)$ $\displaystyle=e^{-8\pi
t}\left(\mathcal{I}(t)+\frac{\mathcal{I}(1)}{e^{8\pi}-1}\right)\begin{pmatrix}\cos^{2}2\pi
t&\cos 2\pi t\sin 2\pi t\\\ \cos 2\pi t\sin 2\pi t&\sin^{2}2\pi
t\end{pmatrix}$ $\displaystyle\quad+b\begin{pmatrix}\sin^{2}2\pi t&-\sin 2\pi
t\cos 2\pi t\\\ -\sin 2\pi t\cos 2\pi t&\cos^{2}2\pi t\end{pmatrix}$ (55)

for some $b$, such that $\lambda^{\mathsf{T}}C\lambda\equiv b$. We obtain a
periodic rescaled covariance function provided that $b=0$. Its eigenvalues
equal

$e^{-8\pi t}\left(\mathcal{I}(t)+\frac{\mathcal{I}(1)}{e^{8\pi}-1}\right)$
(56)

and $0$ with eigenvectors $(\cos 2\pi t,\sin 2\pi t)^{\mathsf{T}}$ and $(-\sin
2\pi t,\cos 2\pi t)^{\mathsf{T}}$, respectively. With the noise defined in
terms of the matrix

$F(x)=\begin{pmatrix}x_{1}x_{2}\\\ x_{2}^{2}\end{pmatrix}$ (57)

and a single standard Brownian motion, we obtain

$\mathcal{I}(t)=\frac{e^{8\pi t}(5-4\cos 4\pi t-2\sin 4\pi t)-1}{40}$ (58)

and the nonzero eigenvalue of the periodic rescaled covariance function equals

$\frac{5-4\cos 4\pi t-2\sin 4\pi t}{40}.$ (59)

We compare the analytical predictions with the approximate results obtained
from an implementation of the periodic boundary-value problem in the software
package coco using continuous, piecewise-polynomial interpolants for
$\gamma(t)$, $\lambda(t)$, and $C(t)$ with base points on a uniform mesh and
constrained by a discretization of the governing differential equations at
collocation nodes associated with each mesh interval. Figure 1 shows the non-
trivial eigenvalue of the covariance matrix obtained using coco and the
analytical expression (59). The results are in excellent agreement.

Figure 1: Comparison of the non-trivial eigenvalue of the covariance matrix
function for the vector field in Eq. (49) obtained from a discretization of
the governing boundary-value problem in Proposition 2.15 in the software
package coco and from the analytical expression reported in (59),
respectively.

Finally, we validate our leading-order analysis by comparing the predictions
with direct numerical simulation of the original Itô SDE. To generate the
stochastic trajectories, we use an explicit Euler Maruyama integrator in
Matlab with time-step $\mathrm{d}t=10^{-4}$, noise intensity $\sigma=0.1$, and
Brownian increments generated using the randn command. The system is
integrated for a duration of 500 periods of the underlying limit cycle. In
Fig. 2, we show the resultant stochastic trajectories projected along the
normalized radial eigenvector $e(\tau)=(\cos 2\pi\tau,\sin
2\pi\tau)^{\mathsf{T}}$ corresponding to the nonzero eigenvalue of the
covariance matrix $C(\tau)$. Specifically, for each point on the stochastic
trajectory $x(t)$, we use the hyperplane definition to compute the
corresponding value of $\tau(t)$ and project the perturbation
$x(t)-\gamma(\tau(t))$ onto $e(\tau(t))$. These results are compared with the
one and two standard deviation predictions computed from $\sigma$ times the
square root of the eigenvalue given in (59). We observe good qualitative
agreement between the predicted variations and those obtained from direct
integration.

Figure 2: Comparison of leading-order theoretical predictions and a numerical
time history describing the influence of small Brownian noise on the dynamics
near a limit cycle of the deterministic vector field Eq. (49) in cartesian
(upper panel) and polar (lower panel) coordinates. Stochastic trajectories
(red dots) were obtained using an Euler-Maruyama scheme with $\sigma=0.1$ and
$\mathrm{d}t=10^{-4}$ for $500$ periods of the limit cycle. Dashed curves
represent the deterministic limit cycle while solid curves represent predicted
deviations from the limit cycle equal to one and two standard deviations,
respectively, computed using (59). In the lower panel, the vertical axis
equals the radial deviation from the limit cycle, since the normalized radial
eigenvector $e(\tau(t))=\gamma(\tau(t))$ and
$\lambda(\tau(t))^{\mathsf{T}}\left(x(t)-\gamma(\tau(t))\right)=0$. There,
filled circles represent the statistical mean (blue) and standard deviation
(black) of simulated data collected in bins of width $\Delta\tau=1/40$.

In the case of periodically excited deterministic dynamics, the analysis in
the previous section carries over in its entirety by assigning the identity
matrix to $Q$ in (42) and (34). Consider, for example, the periodic solution

$\gamma(t)=\begin{pmatrix}\sin 2\pi t\\\ \cos 2\pi t\end{pmatrix}$ (60)

with period $1$ of the linear $1$-periodic dynamical system

$\dot{x}=2\pi f(t,x),\quad\mbox{ where }f(t,x):=\begin{pmatrix}x_{2}\\\
-2x_{2}-x_{1}+2\cos 2\pi t\end{pmatrix},$ (61)

such that $T=2\pi$ and

$f(t,\gamma(t))=\begin{pmatrix}\cos 2\pi t\\\ -\sin 2\pi t\end{pmatrix}.$ (62)

In this case

$X(t)=e^{-2\pi t}\begin{pmatrix}1+2\pi t&2\pi t\\\ -2\pi t&1-2\pi
t\end{pmatrix}$ (63)

with the (sole) eigenvalue of the monodromy matrix $X(1)$ equal to
$e^{-2\pi}$. Since the vector field is non-autonomous, it is not the case that
$X^{-1}(t)f(t,\gamma(t))\equiv f(0,\gamma(0))$, but this property is also no
longer needed for the construction of $Q$.

Now suppose that $F(x)=(0,x_{1})^{\mathsf{T}}$. Substitution in (42) then
yields

$C(0)=\frac{1}{32}\begin{pmatrix}5&-1\\\ -1&1\end{pmatrix}$ (64)

and from (34),

$C(t)=\frac{1}{32}\begin{pmatrix}4+\cos 4\pi t-\sin 4\pi t&-\cos 4\pi t-\sin
4\pi t\\\ -\cos 4\pi t-\sin 4\pi t&4-3\cos 4\pi t-\sin 4\pi t\end{pmatrix}$
(65)

with pair of positive eigenvalues

$\frac{4-\cos 4\pi t-\sin 4\pi t\pm\sqrt{3+2\cos 8\pi t+\sin 8\pi t}}{32}.$
(66)

The panels in Fig. 3 demonstrate excellent agreement between the density
clouds of points of intersection of stochastic trajectories with several
constant-phase sections and ellipses spanned by the eigenvectors of the
corresponding covariance matrix and with radii given by the square roots of
the eigenvalues in (66).

Figure 3: Comparison of leading-order theoretical predictions and a numerical
time history describing the influence of small Brownian noise on the dynamics
near a limit cycle of the harmonically excited, damped, linear oscillator in
Sec. 2.2. The top and bottom panels show samples (red dots) in sections of
constant excitation phase ($t=0.1$ and $t=0.9$, respectively) obtained using
an Euler-Maruyama scheme with $\sigma=0.1$ and $\mathrm{d}t=10^{-4}$ for
$1,000$ periods of excitation. Solid (black) curves represent deviations from
the intersection of the limit cycle (centered black dot) in terms of one and
two standard deviations predicted using (65) and (66). Dashed (blue) curves
represent the one-standard deviation curve predicted from the statistical
covariance matrix for the simulated data.

##  3 Quasiperiodic Tori

A natural generalization of limit cycle dynamics is quasiperiodic motion on a
transversally stable invariant torus, which we treat in this section. As in
the previous section, we first derive a quasiperiodic covariance boundary-
value problem to describe small fluctuations around a deterministic skeleton
solution to leading order in the noise intensity. Next, we validate our
approach against the results of stochastic numerical integration.

###  3.1 Theoretical derivations

Let $\mathbb{S}\sim[0,1]$. For some positive scalar222As in the previous
section, for an autonomous vector field, $T$ is not typically known a priori
and is, instead, solved for simultaneously with $\gamma$ by appending suitable
phase conditions to the quasiperiodic boundary-value problem in $x$ for fixed
$\rho$, e.g., $h_{i}(\gamma(0,0))=0$ for some smooth scalar-valued functions
$h_{i}$, $i=1,\ldots,m+1$. As before, we omit explicit parameter dependence
for notational simplicity. $T$, suppose that $\gamma(\phi,t)\in\mathbb{R}^{n}$
denotes a family of quasiperiodic solutions of the smooth dynamical system

$\dot{x}=Tf(x),\qquad\mbox{ where }x:=x(t)\in\mathbb{R}^{n},$ (67)

on an invariant torus, i.e., such that (67) is satisfied for
$x(t)=\gamma(\phi,t)$ for all $\phi\in\mathbb{S}^{m}$, and there exists a
function $u:\mathbb{S}^{m+1}\rightarrow\mathbb{R}^{n}$ and vector
$\rho\in\mathbb{R}^{m}$ so that $r^{\mathsf{T}}\rho\neq 1$ for all rational
vectors $r\in\mathbb{R}^{m}$ and $\gamma(\phi,t)=u(\phi+\rho t,t)$. In
particular, $\gamma(\phi,t+1)=\gamma(\phi+\rho,t)$ for all $t$. When $m=0$,
the analysis reduces to the case of a periodic orbit.

Let $X(\phi,t)$ denote the solution to the corresponding variational initial-
value problem

$\partial_{t}X=T\mathrm{D}f(\gamma)X,\quad X(\cdot,0)=I_{n}.$ (68)

By quasiperiodicity of $\gamma$, $X(\phi,t+1)=X(\phi+\rho,t)X(\phi,1)$.
Furthermore, it follows that
$X^{-1}(\cdot,t)\partial_{t}\gamma(\cdot,t)\equiv\partial_{t}\gamma(\cdot,0)$
and
$X^{-1}(\cdot,t)\partial_{\phi_{i}}\gamma(\cdot,t)\equiv\partial_{\phi_{i}}\gamma(\cdot,0)$,
$i=1,\ldots,m$, since equality holds at $t=0$ and differentiation with respect
to $t$ followed by substitution of (67) and (68) shows that the left-hand
sides are constant. As a special case, again by quasiperiodicity of $\gamma$,
it follows that $\partial_{t}\gamma(\cdot,0)$ and
$\partial_{\phi_{i}}\gamma(\cdot,0)$, $i=1,\ldots,m$, are right nullvectors of
the operator

$\Gamma_{0}:\delta(\cdot)\mapsto
X(\cdot-\rho,1)\delta(\cdot-\rho)-\delta(\cdot)$ (69)

and that $\partial_{t}\gamma(\cdot,t)$ and
$\partial_{\phi_{i}}\gamma(\cdot,t)$, $i=1,\ldots,m$, are right nullvectors of
the operator

$\Gamma_{t}:\delta(\cdot)\mapsto
X(\cdot,t)\Gamma_{0}\left[X^{-1}(\cdot,t)\delta(\cdot)\right].$ (70)

Analogously with the periodic case and again following [15], assume next that
$\gamma$ is normally hyperbolic in the sense that there exists a
quasiperiodic, smooth family of projections of the form

$Q(\phi,t):=X(\phi,t)\left(I_{n}-\partial_{t}\gamma(\phi,0)w_{t}^{\mathsf{T}}(\phi)-\sum_{i=1}^{m}\partial_{\phi_{i}}\gamma(\phi,0)w_{\phi_{i}}^{\mathsf{T}}(\phi)\right)X^{-1}(\phi,t)$
(71)

for some periodic functions $w_{t}(\phi)$ and $w_{\phi_{i}}(\phi)$,
$i=1,\ldots,m$, such that $\Gamma_{0}$ is a bijection with bounded inverse on
the space of functions $\phi\mapsto Q(\phi,0)\delta(\phi)$ for arbitrary
smooth functions $\delta$ on $\mathbb{S}$. As before, the family $Q(\phi,t)$
and the complementary family $I_{n}-Q(\phi,t)$ respect the spectral
decomposition of $\Gamma_{t}+\mathrm{Id}$.

In further analogy to the periodic case, we formally derive conditions on the
nullspaces of $\Gamma_{0}$; the calculations are formal since we operate with
transposition and integration instead of duality pairings. In particular, we
let

$\int_{\mathbb{S}^{n}}w^{\mathsf{T}}(\phi)\left[\cdot\right]\,\mathrm{d}\phi$
(72)

denote a linear functional on the space of continuous functions on
$\mathbb{S}$ in terms of some continuous, periodic function $w(\phi)$. Then,
the left nullspace of $\Gamma_{0}$ are all functions $w(\phi)$ such that

$\int_{\mathbb{S}^{n}}w^{\mathsf{T}}(\phi)\Gamma_{0}[\delta(\phi)]\,\mathrm{d}\phi=0$
(73)

for arbitrary $\delta(\phi)$. Claims about uniqueness should be understood in
this context as limited to spaces in which such representations are possible.

####  3.1.1 The adjoint boundary-value problem

We mirror the treatment in the previous section by deriving the adjoint
boundary-value problems in Lemma 3.3 and a projection of the local dynamics
onto directions transversal to the torus in Corollary 3.6. The discussion is
adapted from the treatment of two-tori ($m=1$) in [15] .

######  Lemma 3.1.

The linear functionals

$\int_{\mathbb{S}^{m}}w_{t}^{\mathsf{T}}(\phi)\left[\cdot\right]\,\mathrm{d}\phi\mbox{
and
}\int_{\mathbb{S}^{m}}w_{\phi_{i}}^{\mathsf{T}}(\phi)\left[\cdot\right]\,\mathrm{d}\phi,\,i=1,\ldots,m$
(74)

are the unique elements of the left nullspace of $\Gamma_{0}$ such that

$w_{t}^{\mathsf{T}}(\cdot)\partial_{t}\gamma(\cdot,0)\equiv
1,\,w_{t}^{\mathsf{T}}(\cdot)\partial_{\phi_{i}}\gamma(\cdot,0)\equiv 0\equiv
w_{\phi_{i}}^{\mathsf{T}}(\cdot)\partial_{t}\gamma(\cdot,0)$ (75)

for $i=1,\ldots,m$ and

$w_{\phi_{i}}^{\mathsf{T}}(\cdot)\partial_{\phi_{j}}\gamma(\cdot,0)\equiv\delta_{ij}$
(76)

for $i,j=1,\ldots,m$, where $\delta_{ij}=1$ when $i=j$ and $0$ otherwise.

###### Proof.

Let $\ast$ represents either $t$ or $\phi_{i}$. Then, quasiperiodicity
$Q(\phi,1)=Q(\phi+\rho,0)$ holds if and only if
$w_{\ast}^{\mathsf{T}}(\phi)X^{-1}(\phi,1)=w_{\ast}^{\mathsf{T}}(\phi+\rho)$.
But,

$\displaystyle\int_{\mathbb{S}^{m}}w_{\ast}^{\mathsf{T}}(\phi)\Gamma_{0}\left[\delta(\phi)\right]\,\mathrm{d}\phi=\int_{\mathbb{S}^{m}}w_{\ast}^{\mathsf{T}}(\phi)\left(X(\phi-\rho,1)\delta(\phi-\rho)-\delta(\phi)\right)\,\mathrm{d}\phi$
$\displaystyle\qquad\qquad=\int_{\mathbb{S}^{m}}\left(w_{\ast}^{\mathsf{T}}(\phi+\rho)X(\phi,1)-w_{\ast}^{\mathsf{T}}(\phi)\right)\delta(\phi)\,\mathrm{d}\phi$
(77)

for an arbitrary periodic function $\delta(\phi)$. It follows that
quasiperiodicity is equivalent to
$\int_{\mathbb{S}^{m}}w_{\ast}^{\mathsf{T}}(\phi)\left[\cdot\right]\,\mathrm{d}\phi$
lying in the left nullspace of $\Gamma_{0}$. The conditions (75) and (76)
follow by examining the equality $Q^{2}(\phi,0)=Q(\phi,0)$. Finally, to show
uniqueness, suppose that

$\int_{\mathbb{S}^{m}}\tilde{w}^{\mathsf{T}}(\phi)\left[\cdot\right]\,\mathrm{d}\phi$
(78)

is any left nullvector of $\Gamma_{0}$ and let $v(\phi)=Q(\phi,0)\delta(\phi)$
for an arbitrary continuous, periodic function $\delta(\phi)$. Since
$\Gamma_{0}$ is onto functions of this form, there exists a periodic function
$u(\phi)$ such that $\Gamma_{0}\left[u(\phi)\right]=v(\phi)$ and,
consequently,

$\int_{\mathbb{S}^{m}}\tilde{w}^{\mathsf{T}}(\phi)Q(\phi,0)\delta(\phi)\,\mathrm{d}\phi=\int_{\mathbb{S}^{m}}\tilde{w}^{\mathsf{T}}(\phi)v(\phi)\,\mathrm{d}\phi=\int_{\mathbb{S}^{m}}\tilde{w}^{\mathsf{T}}(\phi)\Gamma_{0}\left[u(\phi)\right]\,\mathrm{d}\phi=0.$
(79)

But since $\delta$ is arbitrary, this implies that
$\tilde{w}^{\mathsf{T}}(\phi)Q(\phi,0)\equiv 0$ and, consequently, that
$\tilde{w}(\phi)\in\mathrm{span}\\{w_{t}(\phi),w_{\phi_{1}}(\phi),\ldots,w_{\phi_{m}}(\phi)\\}$
for each $\phi$. The claim follows with the imposition of (75) and (76). ∎

######  Corollary 3.2.

The right nullspace of $\Gamma_{0}$ consists of functions that are tangential
to the manifold $(\phi,t)\mapsto\gamma(\phi,t)$ at every $\phi$ and $t=0$.

###### Proof.

In the same notation as in Lemma 3.1,

$X(\phi-\rho,1)\partial_{\ast}\gamma(\phi-\rho,0)w_{\ast}^{\mathsf{T}}(\phi-\rho)=\partial_{\ast}\gamma(\phi,0)w_{\ast}^{\mathsf{T}}(\phi)X(\phi-\rho,1)$
(80)

and, consequently,
$\Gamma_{t}\left[Q(\phi,t)\cdot\right]=Q(\phi,t)\Gamma_{t}\left[\cdot\right]$.
Now suppose that $\delta(\phi)$ is a right nullvector of $\Gamma_{0}$. It
follows that
$\Gamma_{0}\left[Q(\phi,0)\delta(\phi)\right]=Q(\phi,0)\Gamma_{0}\left[\delta(\phi)\right]=0$.
Since $\Gamma_{0}$ is one-to-one on the space of functions $\phi\mapsto
Q(\phi,0)\delta(\phi)$, $Q(\phi,0)\delta(\phi)\equiv 0$ and, consequently,

$\delta(\phi)=\mathrm{span}\\{\partial_{t}\gamma(\phi,0),\partial_{\phi_{1}}\gamma(\phi,0),\ldots,\partial_{\phi_{m}}\gamma(\phi,0)\\},$
(81)

for each $\phi$. ∎

######  Lemma 3.3.

The functions

$\lambda_{t}(\phi,t):=\left(X^{-1}\right)^{\mathsf{T}}(\phi,t)w_{t}(\phi)\mbox{
and
}\lambda_{\phi_{i}}(\phi,t):=\left(X^{-1}\right)^{\mathsf{T}}(\phi,t)w_{\phi_{i}}(\phi)$
(82)

for $i=1,\ldots,m$ are the unique solutions to the adjoint boundary-value
problem

$0=-\partial_{t}\mu^{\mathsf{T}}-T\mu^{\mathsf{T}}\mathrm{D}f(\gamma),\,0=\mu(\phi+\rho,0)-\mu(\phi,1)$
(83)

where

$1=\int_{\mathbb{S}^{m+1}}\mu^{\mathsf{T}}(\phi,t)\partial_{t}\gamma(\phi,t)\,\mathrm{d}t\,\mathrm{d}\phi$
(84)

and

$0=\int_{\mathbb{S}^{m}}\mu^{\mathsf{T}}(\phi,1)\partial_{\phi_{j}}\gamma(\phi,1)\,\mathrm{d}\phi,\,j=1,\ldots,m$
(85)

in the case of $\lambda_{t}$, and

$0=\int_{\mathbb{S}^{m+1}}\mu^{\mathsf{T}}(\phi,t)\partial_{t}\gamma(\phi,t)\,\mathrm{d}t\,\mathrm{d}\phi$
(86)

and

$\delta_{ij}=\int_{\mathbb{S}^{m}}\mu^{\mathsf{T}}(\phi,1)\partial_{\phi_{j}}\gamma(\phi,1)\,\mathrm{d}\phi,\,j=1,\ldots,m$
(87)

in the case of $\lambda_{\phi_{i}}$.

###### Proof.

Let $\mu$ denote an arbitrary solution to the boundary-value problem. It
follows that $\mu^{\mathsf{T}}(\phi,t)\partial_{t}\gamma(\phi,t)\equiv
v_{t}(\phi)$, $\mu^{\mathsf{T}}(\phi,t)\partial_{\phi_{i}}\gamma(\phi,t)\equiv
v_{\phi_{i}}(\phi)$, $i=1,\ldots,m$, and
$\mu^{\mathsf{T}}(\phi,t)X(\phi,t)\equiv V(\phi)$ for some periodic functions
$v_{t}(\phi)$, $v_{\phi_{i}}(\phi)$, $i=1,\ldots,m$, and $V(\phi)$, since
differentiation with respect to $t$ followed by substitution of (67), (68),
and (83) shows that the left-hand sides are independent of $t$. By
quasiperiodicity,

$v_{t}(\phi-k\rho)=\mu^{\mathsf{T}}(\phi,0)\partial_{t}\gamma(\phi,0)\mbox{
and
}v_{\phi_{i}}(\phi-k\rho)=\mu^{\mathsf{T}}(\phi,0)\partial_{\phi_{i}}\gamma(\phi,0)$
(88)

for any integer $k$. Since $\rho$ is irrational, $v_{t}$ and $v_{\phi_{i}}$
must be independent of $\phi$. Additionally,

$V^{\mathsf{T}}(\phi+\rho)X(\phi,1)-V^{\mathsf{T}}(\phi)=\mu^{\mathsf{T}}(\phi+\rho,0)X(\phi,1)-\mu^{\mathsf{T}}(\phi,1)X(\phi,1)=0$
(89)

and it follows that

$\int_{\mathbb{S}^{m}}V^{\mathsf{T}}(\phi)\left[\cdot\right]\,\mathrm{d}\phi$
(90)

lies in the left nullspace of $\Gamma_{0}$. ∎

######  Corollary 3.4.

In the notation of the previous lemma,

$Q(\phi,t)=I_{n}-\partial_{t}\gamma(\phi,t)\lambda_{t}^{\mathsf{T}}(\phi,t)-\sum_{i=1}^{m}\partial_{\phi_{i}}\gamma(\phi,t)\lambda_{\phi_{i}}^{\mathsf{T}}(\phi,t)$
(91)

is a projection along
$\mathrm{span}\\{\partial_{t}\gamma(\phi,t),\partial_{\phi_{1}}\gamma(\phi,t),\ldots,\partial_{\phi_{m}}\gamma(\phi,t)\\}$
and onto a hyperplane orthogonal to $\lambda_{t}(\phi,t)$ and
$\lambda_{\phi_{i}}(\phi,t)$ for $i=1,\ldots,m$.

We again investigate the geometry on a neighborhood of $\gamma$ in terms of
the projections $Q$. These provide a local decomposition into state-space
displacements tangential and transversal to the corresponding invariant torus.

######  Lemma 3.5.

For $x$ sufficiently close to $\gamma(\phi,t)$ for some $\phi$ and $t$, there
exist unique $\psi$ (with $\psi\approx\phi$) and $\tau$ (with $\tau\approx
t$), such that

$x=\gamma(\psi,\tau)+x_{\mathrm{tr}},$ (92)

where $\lambda_{t}(\psi,\tau)^{\mathsf{T}}x_{\mathrm{tr}}=0$ and
$\lambda_{\phi_{i}}(\psi,\tau)^{\mathsf{T}}x_{\mathrm{tr}}=0$ and,
consequently, $Q(\psi,\tau)x_{\mathrm{tr}}=x_{\mathrm{tr}}$.

###### Proof.

Let

$H(x,\psi,\tau)=\begin{pmatrix}\lambda_{t}^{\mathsf{T}}(\psi,\tau)\left(x-\gamma(\psi,\tau)\right)\\\
\lambda_{\phi_{1}}^{\mathsf{T}}(\psi,\tau)\left(x-\gamma(\psi,\tau)\right)\\\
\vdots\\\
\lambda_{\phi_{m}}^{\mathsf{T}}(\psi,\tau)\left(x-\gamma(\psi,\tau)\right)\end{pmatrix}$
(93)

It follows that $H(\gamma(\phi,t),\phi,t)=0$ and

$\partial_{\psi,\tau}H(\gamma(\phi,t),\phi,t)=\begin{pmatrix}0&-1\\\
-I_{m}&0\end{pmatrix}.$ (94)

The claim follows from the implicit function theorem. ∎

As in the previous section, we perform formal calculations with differentials
in order to capture leading-order effects in the derivation of a suitable
covariance function and the corresponding boundary-value problem.

######  Corollary 3.6.

For the differential $\mathrm{d}x$, Lemma 3.5 implies the existence of
differentials $\mathrm{d}\tau$ and $\mathrm{d}x_{\mathrm{tr}}$, such that

$x+\mathrm{d}x=\gamma(\psi+\mathrm{d}\psi,\tau+\mathrm{d}\tau)+x_{\mathrm{tr}}+\mathrm{d}x_{\mathrm{tr}},$
(95)

where
$\lambda_{t}^{\mathsf{T}}(\psi+\mathrm{d}\psi,\tau+\mathrm{d}\tau)(x_{\mathrm{tr}}+\mathrm{d}x_{\mathrm{tr}})=0$
and
$\lambda_{\phi_{i}}^{\mathsf{T}}(\psi+\mathrm{d}\psi,\tau+\mathrm{d}\tau)(x_{\mathrm{tr}}+\mathrm{d}x_{\mathrm{tr}})=0$
for $i=1,\ldots,m$, and, consequently,

$Q(\psi+\mathrm{d}\psi,\tau+\mathrm{d}\tau)(x_{\mathrm{tr}}+\mathrm{d}x_{\mathrm{tr}})=x_{\mathrm{tr}}+\mathrm{d}x_{\mathrm{tr}}.$
(96)

It follows that

$\displaystyle\mathrm{d}x_{\mathrm{tr}}$
$\displaystyle-T\mathrm{D}f(\gamma(\psi,\tau))x_{\mathrm{tr}}\mathrm{d}\tau-\sum_{i=1}^{m}\partial_{\phi_{i}}Q(\psi,\tau)x_{\mathrm{tr}}\mathrm{d}\psi_{i}$
$\displaystyle\qquad=Q(\psi,\tau)\left(\mathrm{d}x-T\mathrm{D}f(\gamma(\psi,\tau))x_{\mathrm{tr}}\mathrm{d}\tau\right).$
(97)

###### Proof.

Eq. (95) implies that

$\mathrm{d}x_{\mathrm{tr}}=\mathrm{d}x-\partial_{t}\gamma(\psi,\tau)\mathrm{d}\tau-\sum_{i=1}^{m}\partial_{\phi_{i}}\gamma(\psi,\tau)\mathrm{d}\psi_{i}.$
(98)

Similarly, from (96), we obtain

$\mathrm{d}x_{\mathrm{tr}}=\partial_{t}Q(\psi,\tau)x_{\mathrm{tr}}\mathrm{d}\tau+\sum_{i=1}^{m}\partial_{\phi_{i}}Q(\psi,\tau)x_{\mathrm{tr}}\mathrm{d}\psi_{i}+Q(\psi,\tau)\mathrm{d}x_{\mathrm{tr}}.$
(99)

From the definition of $Q$ it follows that

$\displaystyle\partial_{t}Q(\psi,\tau)$
$\displaystyle=-T\mathrm{D}f(\gamma(\psi,\tau))\left(\partial_{t}\gamma(\psi,\tau)\lambda_{t}^{\mathsf{T}}(\psi,\tau)+\sum_{i=1}^{m}\partial_{\phi_{i}}\gamma(\psi,\tau)\lambda_{\phi_{i}}^{\mathsf{T}}(\psi,\tau)\right)$
$\displaystyle\qquad+T\left(\partial_{t}\gamma(\psi,\tau)\lambda_{t}^{\mathsf{T}}(\psi,\tau)+\sum_{i=1}^{m}\partial_{\phi_{i}}\gamma(\psi,\tau)\lambda_{\phi_{i}}^{\mathsf{T}}(\psi,\tau)\right)\mathrm{D}f(\gamma(\psi,\tau))$
(100)

and, consequently, that

$\partial_{t}Q(\psi,\tau)x_{\mathrm{tr}}=T(I_{n}-Q(\psi,\tau))\mathrm{D}f(\gamma(\psi,\tau))x_{\mathrm{tr}},$
(101)

since $\lambda_{t}^{\mathsf{T}}(\psi,\tau)x_{\mathrm{tr}}=0$ and
$\lambda_{\phi_{i}}^{\mathsf{T}}(\psi,\tau)x_{\mathrm{tr}}=0$ for
$i=1,\ldots,m$. Substitution of (98) and (101) into the right-hand side of
(99) yields (3.6), since $Q(\psi,\tau)\partial_{t}\gamma(\psi,\tau)=0$ and
$Q(\psi,\tau)\partial_{\phi_{i}}\gamma(\psi,\tau)=0$ for $i=1,\ldots,m$. ∎

####  3.1.2 The covariance boundary-value problem

In lieu of the deterministic dynamics in (67), consider again the Itô SDE

$\mathrm{d}x=Tf(x)\mathrm{d}t+\sigma\sqrt{T}F(x)\mathrm{d}W_{t}$ (102)

in terms of the noise intensity $\sigma$, and a vector
$W_{t}\in\mathbb{R}^{m}$ of independent standard Brownian motions. In the
notation of the preceding lemmas, and as before under a small noise assumption
and up to leading-order terms, the following proposition holds.

######  Proposition 3.7.

For $x$ sufficiently close to $\gamma(\phi,t)$ for some $\phi$ and $t$, to
first order in $\|x_{\mathrm{tr}}\|$ and $\sigma$,

$\mathrm{d}x_{\mathrm{tr}}=T\mathrm{D}f(\gamma(\psi,\tau))x_{\mathrm{tr}}\mathrm{d}\tau+\sigma\sqrt{T}Q(\psi,\tau)F(\gamma(\psi,\tau))\mathrm{d}W_{\tau}.$
(103)

###### Proof.

For $x_{\mathrm{tr}}=0$ and $\sigma=0$, we obtain $\tau=t$, $\psi=\phi$, and
$\mathrm{d}x_{\mathrm{tr}}=0$. In this limit, (102) and (98) imply that

$0=\partial_{t}\gamma(\psi,\tau)(\mathrm{d}t-\mathrm{d}\tau)-\sum_{i=1}^{m}\partial_{\phi_{i}}\gamma(\psi,\tau)\mathrm{d}\psi_{i},$
(104)

from which it follows that $\mathrm{d}\tau=\mathrm{d}t$ and
$\mathrm{d}\psi_{i}=0$, $i=1,\ldots,m$, to zeroth order in
$\|x_{\mathrm{tr}}\|$ and $\sigma$. To first order in $\|x_{\mathrm{tr}}\|$
and $\sigma$, (102) then yields

$Q(\psi,\tau)\left(\mathrm{d}x-T\mathrm{D}f(\gamma(\psi,\tau))x_{\mathrm{tr}}\mathrm{d}\tau\right)=\sigma\sqrt{T}Q(\psi,\tau)F(\gamma(\psi,\tau))\mathrm{d}W_{\tau}$
(105)

and substitution in (3.6) yields the claimed result. ∎

The Ornstein-Uhlenbeck process in (103) approximates the local dynamics near
$\gamma$ for orbits that remain near $\gamma$ to leading order but neglects
cumulative effects of higher-order terms, e.g., differences between
$\mathrm{d}t$ and $\mathrm{d}\tau$ and higher-order terms in $\mathrm{d}\psi$.
Given an initial condition $x_{\mathrm{tr}}(\psi,0)$, (103) may be solved
explicitly to yield

$\displaystyle x_{\mathrm{tr}}(\psi,\tau)$
$\displaystyle=X(\psi,\tau)x_{\mathrm{tr}}(\psi,0)+\sigma\sqrt{T}X(\psi,\tau)\int_{0}^{\tau}X^{-1}(\psi,s)Q(\psi,s)F(\gamma(\psi,s)))\mathrm{d}W_{s}$
$\displaystyle=X(\psi,\tau)x_{\mathrm{tr}}(\psi,0)+\sigma\sqrt{T}Q(\psi,\tau)X(\psi,\tau)\int_{0}^{\tau}G(\psi,s)\mathrm{d}W_{s},$
(106)

where $G(\psi,s)=X^{-1}(\psi,s)F(\gamma(\psi,s))$. With

$\Lambda(\psi,\tau):=\begin{pmatrix}\lambda_{t}(\psi,\tau)&\lambda_{\phi_{1}}(\psi,\tau)&\cdots&\lambda_{\phi_{m}}(\psi,\tau)\end{pmatrix}$
(107)

and

$\Omega(\psi):=\begin{pmatrix}w_{t}(\psi)&w_{\phi_{1}}(\psi)&\cdots&w_{\phi_{m}}(\psi)\end{pmatrix}$
(108)

it follows that
$\Lambda^{\mathsf{T}}(\psi,\tau)x_{\mathrm{tr}}(\psi,\tau)\equiv 0$ provided
that $\Omega^{\mathsf{T}}(\psi)x_{\mathrm{tr}}(\psi,0)=0$. The rescaled
covariance matrix

$C:=\frac{1}{\sigma^{2}}\mathbb{E}\left[x_{\mathrm{tr}}x_{\mathrm{tr}}^{\mathrm{T}}\right]$
(109)

is then given by

$\displaystyle C(\psi,\tau)$
$\displaystyle=X(\psi,\tau)C(\psi,0)X^{\mathsf{T}}(\psi,\tau)$
$\displaystyle\qquad+TQ(\psi,\tau)X(\psi,\tau)\left(\int_{0}^{\tau}G(\psi,s)G^{\mathsf{T}}(\psi,s)\mathrm{d}s\right)X^{\mathsf{T}}(\psi,\tau)Q^{\mathsf{T}}(\psi,\tau),$
(110)

It follows that
$\Lambda^{\mathsf{T}}(\psi,\tau)C(\psi,\tau)\Lambda(\psi,\tau)$ is independent
of $\tau$ and equals $0$ if and only if
$\Omega(\psi)^{\mathsf{T}}C(\psi,0)\Omega(\psi)=0$.

######  Lemma 3.8.

Let

$\mathcal{I}(\psi):=\int_{0}^{1}G(\psi,s)G^{\mathsf{T}}(\psi,s)\mathrm{d}s$
(111)

Then, for arbitrary $C(\psi,0)$,

$\displaystyle
C(\psi,k+1)=X(\psi+k\rho,1)C(\psi,k)X^{\mathsf{T}}(\psi+k\rho,1)$
$\displaystyle\qquad+TQ(\psi+k\rho,1)X(\psi+k\rho,1)\mathcal{I}(\psi+k\rho)X^{\mathsf{T}}(\psi+k\rho,1)Q^{\mathsf{T}}(\psi+k\rho,1).$
(112)

###### Proof.

By (3.1.2),

$\displaystyle C(\psi,k)=X(\psi,k)C(\psi,0)X^{\mathsf{T}}(\psi,k)$
$\displaystyle\qquad+TQ(\psi,k)X(\psi,k)\left(\int_{0}^{k}G(\psi,s)G^{\mathsf{T}}(\psi,s)\mathrm{d}s\right)X^{\mathsf{T}}(\psi,k)Q^{\mathsf{T}}(\psi,k).$
(113)

Moreover, by the definition of $G$,

$\mathcal{I}(\psi+k\rho)=X(\psi,k)\left(\int_{k}^{k+1}G(\psi,s)G^{\mathsf{T}}(\psi,s)\mathrm{d}s\right)X^{\mathsf{T}}(\psi,k).$
(114)

The claim follows by substitution. ∎

We proceed to assume that the quasiperiodic orbit is transversally stable,
such that

$\displaystyle X(\psi,k)$
$\displaystyle=\partial_{t}\gamma(\psi+k\rho,0))w_{t}^{\mathsf{T}}(\psi)$
$\displaystyle\qquad+\sum_{i=1}^{m}\partial_{\phi_{i}}\gamma(\psi+k\rho,0))w_{\phi_{i}}^{\mathsf{T}}(\psi)+\mathcal{O}\left(\mathrm{exp}(-k/\tau_{\mathrm{tr}})\right)$
(115)

for some positive constant $\tau_{\mathrm{tr}}$ and large $k$. Then,

######  Proposition 3.9.

Given a constant matrix $B$, there exists an initial condition $C(\psi,0)$
with $\Omega(\psi)^{\mathsf{T}}C(\psi,0)\Omega(\psi)=B$, such that
$C(\psi,k)=C(\psi+k\rho,0)$ for all $k$.

###### Proof.

The previous lemma shows that $C(\psi,1)=C(\psi+\rho,0)\Rightarrow
C(\psi,k)=C(\psi+k\rho,0)$ for all $k$. Suppose that

$\displaystyle C(\psi,0)$
$\displaystyle=T\sum_{n=1}^{\infty}Q(\psi,0)X(\psi-n\rho,n)\mathcal{I}(\psi-n\rho)X^{\mathsf{T}}(\psi-n\rho,n)Q^{\mathsf{T}}(\psi,0)$
$\displaystyle\quad+\nabla(\psi,0)B\nabla^{\mathsf{T}}(\psi,0).$ (116)

where

$\nabla(\psi,\tau):=\begin{pmatrix}\partial_{t}\gamma(\psi,\tau)&\partial_{\phi_{1}}\gamma(\psi,\tau)&\cdots&\partial_{\phi_{m}}\gamma(\psi,\tau)\end{pmatrix}$
(117)

The series converges by the assumption of transversal stability. Moreover,

$\displaystyle X(\psi,1)C(\psi,0)X^{\mathsf{T}}(\psi,1)=C(\psi+\rho,0)$
$\displaystyle\qquad-
TQ(\psi+\rho,0)X(\psi,1)\mathcal{I}(\psi)X^{\mathsf{T}}(\psi,1)Q^{\mathsf{T}}(\psi+\rho,0)$
(118)

and the claim follows by evaluation of (3.1.2) at $\tau=1$. ∎

######  Corollary 3.10.

Given a constant matrix $B$, the function $C(\psi,\tau)$ in (3.1.2) with
$\Omega(\psi)^{\mathsf{T}}C(\psi,0)\Omega(\psi)=\beta$ converges to the
function obtained with $C(\psi,0)$ given by (3.1.2) as
$\tau\rightarrow\infty$.

######  Proposition 3.11.

Let $C_{\mathrm{qper}}$ denote the function in (3.1.2) obtained with
$C(\psi,0)$ given by (3.1.2) when $B=0$. It follows that, when the matrix
$A=0$, $C_{\mathrm{qper}}$ is the unique quasiperiodic solution of the
quasiperiodic Lyapunov equation (cf. [24])

$\displaystyle\partial_{t}C$
$\displaystyle=T\left(\mathrm{D}f(\gamma)C+C\mathrm{D}f^{\mathsf{T}}(\gamma)+QF(\gamma)F^{\mathsf{T}}(\gamma)Q^{\mathsf{T}}+\nabla
A\nabla^{\mathsf{T}}\right),$ (119)

for which $\Omega(\psi)^{\mathsf{T}}C(\psi,0)\Omega(\psi)=0$ and that no
quasiperiodic solution exists otherwise.

###### Proof.

It follows from (83) that

$\frac{d}{dt}\Lambda(\psi,t)^{\mathsf{T}}C(\psi,t)\Lambda(\psi,t)=AT$ (120)

Thus, $C$ is quasiperiodic only if $A=0$. The transformation
$C=X\tilde{C}X^{\mathsf{T}}$ then implies that

$\partial_{t}\tilde{C}(\psi,t)=TQ(\psi,0)G(\psi,t)G^{\mathsf{T}}(\psi,t)Q^{\mathsf{T}}(\psi,0)$
(121)

and the conclusion follows by integration. ∎

We conclude by noting that

$\Omega(0)^{\mathsf{T}}C(0,0)\Omega(0)=0\Rightarrow\Lambda^{\mathsf{T}}(0,k)C(0,k)\Lambda(0,k)=0$
(122)

Assuming quasiperiodicity, it follows that

$\Omega(k\rho)^{\mathsf{T}}C(k\rho,0)\Omega(k\rho)=0$ (123)

for every $k$. Since $\rho$ is irrational and by continuity, it follows that
$\Omega(\phi)^{\mathsf{T}}C(\phi,0)\Omega(\phi)=0$ for all $\phi$.

###  3.2 Numerical examples

As a first example, we take inspiration from the second example in Sec. 2.2
and consider the non-autonomous, two-dimensional system of SDEs

$\displaystyle\mathrm{d}x_{1}$ $\displaystyle=\frac{2\pi}{\omega}\left(-\Omega
x_{2}+x_{1}\left(1+\sqrt{x_{1}^{2}+x_{2}^{2}}(\cos 2\pi
t-1)\right)\right)\mathrm{d}t+\sigma\sqrt{\frac{2\pi}{\omega}}x_{1}x_{2}\mathrm{d}W_{t}$
(124) $\displaystyle\mathrm{d}x_{2}$
$\displaystyle=\frac{2\pi}{\omega}\left(\Omega
x_{1}+x_{2}\left(1+\sqrt{x_{1}^{2}+x_{2}^{2}}(\cos 2\pi
t-1)\right)\right)\mathrm{d}t+\sigma\sqrt{\frac{2\pi}{\omega}}x_{2}^{2}\mathrm{d}W_{t}$
(125)

of the form (102) with $T=2\pi/\omega$, $f$ given by the $1$-periodic vector
field

$f(t,x)=\begin{pmatrix}-\Omega
x_{2}+x_{1}\left(1+\sqrt{x_{1}^{2}+x_{2}^{2}}(\cos 2\pi t-1)\right)\\\ \Omega
x_{1}+x_{2}\left(1+\sqrt{x_{1}^{2}+x_{2}^{2}}(\cos 2\pi
t-1)\right)\end{pmatrix},$ (126)

and $F(x)=(x_{1}x_{2},x_{2}^{2})^{\mathsf{T}}$. Rather than consider the
suspended autonomous state space obtained by appending a third component to
the vector field equal to $1$, we retain the two-dimensional description but
omit $w_{t}$ and $\lambda_{t}$ in the definition of $Q$.

A quasiperiodic solution of the deterministic dynamics is given by

$\gamma(\phi,t)=\frac{1+\omega^{2}}{1+\omega^{2}-\cos 2\pi t-\omega\sin 2\pi
t}\begin{pmatrix}\cos 2\pi(\phi+\rho t)\\\ \sin 2\pi(\phi+\rho
t)\end{pmatrix}$ (127)

provided that $\rho=\Omega/\omega$ is irrational. With the notation
$b:=1+\omega^{2}-\cos 2\pi t-\omega\sin 2\pi t$, $c_{\phi,t}:=\cos
2\pi(\phi+\rho t)$, and $s_{\phi,t}:=\sin 2\pi(\phi+\rho t)$, we obtain

$X(\phi,t)=\frac{\omega^{2}}{b^{2}}\begin{pmatrix}e^{-2\pi
t/\omega}\omega^{2}c_{\phi,0}c_{\phi,t}+bs_{\phi,0}s_{\phi,t}&e^{-2\pi
t/\omega}\omega^{2}s_{\phi,0}c_{\phi,t}-bc_{\phi,0}s_{\phi,t}\\\ e^{-2\pi
t/\omega}\omega^{2}c_{\phi,0}s_{\phi,t}-bs_{\phi,0}c_{\phi,t}&e^{-2\pi
t/\omega}\omega^{2}s_{\phi,0}s_{\phi,t}+bc_{\phi,0}c_{\phi,t}\end{pmatrix}$
(128)

such that, in particular,
$\partial_{\phi}\gamma(\phi,t)=X(\phi,t)\partial_{\phi}\gamma(\phi,0)$ while
it no longer holds that
$\partial_{t}\gamma(\phi,t)=X(\phi,t)\partial_{t}\gamma(\phi,0)$, since $f$ is
non-autonomous. It is now straightforward to verify the equalities
$w_{\phi}^{\mathsf{T}}(\phi)X^{-1}(\phi,1)=w_{\phi}^{\mathsf{T}}(\phi+\rho)$
and $w_{\phi}^{\mathsf{T}}(\phi)\partial_{\phi}\gamma(\phi,0)\equiv 1$ for

$w_{\phi}(\phi)=\frac{\omega^{2}}{2\pi(1+\omega^{2})}\begin{pmatrix}-\sin
2\pi\phi\\\ \cos 2\pi\phi\end{pmatrix}$ (129)

and we obtain the family of projections

$Q(\phi,t)=X(\phi,t)\left(I_{2}-\partial_{\phi}\gamma(\phi,0)w_{\phi}^{\mathsf{T}}(\phi)\right)X^{-1}(\phi,t)=\begin{pmatrix}c_{\phi,t}^{2}&s_{\phi,t}c_{\phi,t}\\\
s_{\phi,t}c_{\phi,t}&s_{\phi,t}^{2}\end{pmatrix}$ (130)

onto the radial direction at each point of the corresponding torus.
Substitution in (3.1.2) (with $B=0$), followed by extensive algebraic
manipulation using Mathematica, then yields the family

$C(\phi,0)=\frac{\left(1+\omega^{2}\right)^{4}\left(1+\Omega^{2}-c_{2\phi,0}-\Omega
s_{2\phi,0}\right)}{4\omega^{8}\left(1+\Omega^{2}\right)}Q(\phi,0).$ (131)

In Fig. 4, we compare this prediction to density clouds of points along
stochastic trajectories sampled at $t\in\mathbb{Z}$ and observe excellent
agreement.

Figure 4: Comparison of leading-order theoretical predictions and a numerical
time history for the first example in Sec. 3.2 with $\Omega=\pi$ and
$\omega=1$ in cartesian (upper panel) and polar (lower panel) coordinates.
Stochastic trajectories sampled at integer values of the excitation phase (red
dots) were obtained using an Euler-Maruyama scheme with $\sigma=0.1$ and
$\mathrm{d}t=10^{-4}$ for $10,000$ periods of the excitation. Dashed curves
represent the intersection of the deterministic quasiperiodic invariant torus
while solid curves represent predicted deviations from this curve of
intersection equal to one and two standard deviations, respectively, computed
using (131). In the lower panel, the vertical axis equals the radial deviation
from the (circular) curve of intersection, since the normalized radial
eigenvector $e(\psi,0)\parallel\gamma(\psi,0)$ and
$w_{\phi}(\psi)^{\mathsf{T}}\left(x-\gamma(\psi,0)\right)=0$. There, filled
circles represent the statistical mean (blue) and standard deviation (black)
of simulated data collected in bins of width $\Delta\psi=1/40$.

As a second example, this time not amenable to closed-form analysis, we
consider the four-dimensional SDE obtained by adding noise to two coupled Van
der Pol oscillators:

$\displaystyle\mathrm{d}y=\left(\begin{array}[]{c}y_{2}\\\
-\epsilon\left(y_{1}^{2}-1\right)y_{2}-y_{1}+\beta\left(y_{3}-y_{1}\right)\\\
y_{4}\\\
-\epsilon\left(y_{3}^{2}-1\right)y_{4}-\left(1+\delta\right)y_{3}+\beta\left(y_{1}-y_{3}\right)\end{array}\right)\mathrm{d\eta}+\sigma\begin{pmatrix}0\\\
\mathrm{d}W_{\eta,1}\\\ 0\\\ \mathrm{d}W_{\eta,2}\end{pmatrix}.$ (136)

Here $\epsilon,\delta,\beta$ are the problem parameters and
$W_{\eta}\in\mathbb{R}^{2}$ is a vector of independent standard Brownian
motions. Quasiperiodic invariant tori of the deterministic drift vector field
have been investigated by Schilder et al. [40].

In the limit of small $\epsilon$ and $\beta=0$, a quasiperiodic solution of
the rescaled deterministic dynamics with $T=2\pi$ is approximately given by

$\gamma(\phi,t)=\begin{pmatrix}2\sin 2\pi t\\\ 2\cos 2\pi t\\\ 2\sin
2\pi(\phi+\rho t)\\\ 2\rho\cos 2\pi(\phi+\rho t)\end{pmatrix}$ (137)

for irrational $\rho=\sqrt{1+\delta}$. We may use this expression and
corresponding value of $T$ to construct an initial solution guess for
quasiperiodic solutions of the rescaled deterministic dynamics for nonzero
$\epsilon$ with $\beta=0$. Parameter continuation may then be used to track a
family of such solutions for fixed $\rho$ and simultaneous variations of
$\beta$ and $\delta$. An example of such a family for $\rho=140/62\sqrt{2}$
and $\epsilon=0.5$ (computed with coco using the discretization algorithm
described in the appendix with $N=14$ and individual trajectory segments
approximated by continuous piecewise-polynomial functions on 20 mesh
intervals) together with a projection of the invariant torus obtained for
$\beta=0.5$ are shown in Fig. 5.

Figure 5: (Upper panel) Family of quasiperiodic invariant tori for the
deterministic limit of the SDE (136) under simultaneous variations in $\beta$
and $\delta$ computed using coco. Here, $\epsilon=0.5$ and the rotation number
$\rho$ is fixed at $140/62\sqrt{2}$. Each torus is represented by a finite
collection of equal-duration trajectory segments coupled through all-to-all
boundary conditions in terms of the rotation number. (Lower panel) The
quasiperiodic invariant torus obtained for $\beta=0.5$ represented in terms of
29 trajectory segments, each of which is approximated by a continuous,
piecewise-polynomial function on 20 mesh intervals.

In the absence of closed-form analysis, we use built-in support in coco for
the corresponding discretization of the adjoint boundary-value problem (see
the Appendix) to compute numerical approximations of $\lambda_{t}(\phi,t)$,
$\lambda_{\phi}(\phi,t)$, and $Q(\phi,t)$ for the torus obtained for
$\beta=0.5$. Finally, we derive a corresponding discretization of the linear
boundary-value problem in Proposition 3.11 and solve for the covariance matrix
function $C(\phi,t)$. We use the norm of the computed matrix $A$ to verify the
accuracy of the computation, as this is predicted to equal $0$ in the original
boundary-value problem. The results of this analysis are visualized in Fig. 6.
In this case, the 2-norm of $A$ equals $9\times 10^{-9}$, giving us confidence
in the validity of the analysis.

Figure 6: Surfaces of the two nonzero eigenvalues for the covariance matrix
function $C(\phi,t)$ corresponding to the quasiperiodic invariant torus
obtained for the deterministic limit of the SDE (136) for $\epsilon=0.5$,
$\beta=0.5$, $\delta=1.9422$ and the rotation number $\rho=140/62\sqrt{2}$.

As in the previous examples, we use points along numerically integrated
stochastic trajectories to validate the theoretical predictions. For points on
such trajectories, we first use the constraints in Lemma 3.5 to identify the
corresponding values of $\psi$ and $\tau$. We then proceed to locate
“intersections” with $\tau=0.5$ by finding consecutive trajectory points on
opposite sides of this section and use linear interpolation through these
points to capture a point on $\tau=0.5$ and compute the corresponding value of
$\psi$ (see Fig. 7). Finally, we project the resultant deviations
$x_{\mathrm{tr}}$ onto the predicted eigenvectors of the corresponding
covariance matrix associated with the two nonzero eigenvalues (see Fig. 8).
The agreement between the statistics for the resultant point clouds and the
predicted standard deviations is excellent.

Figure 7: A point cloud of $5,800$ “intersections” (blue) along a numerically
integrated stochastic trajectory of the SDE (136) with the section $\tau=0.5$
overlaid on the quasiperiodic invariant torus for the corresponding
deterministic dynamics obtained for $\beta=0.5$ (see the lower panel of Fig.
5) and the corresponding points $\gamma(\psi,0.5)$ (red). Stochastic
trajectories were obtained using an Euler-Maruyama scheme with $\sigma=0.1$
and $\mathrm{d}t=10^{-4}$ and associated with points on the torus per the
method described in the text.

Figure 8: Comparison of leading-order theoretical predictions and the results
of analysis of the point cloud in Fig. 7 obtained for the SDE (136). For each
point $x$ in this cloud, the deviation $x-\gamma(\psi(x),0.5)$ is projected
onto the normalized eigenvectors $e_{1}(\psi(x),0.5)$ and $e_{2}(\psi(x),0.5)$
corresponding to the two nonzero eigenvalues of the predicted covariance
matrix $C(\psi(x),0.5)$. Dashed curves represent points on the deterministic
quasiperiodic invariant torus while solid curves represent predicted
deviations from this curve of intersection equal to one and two standard
deviations, respectively, obtained by multiplying the square roots of the
eigenvalues by $\sigma=0.1$. Filled circles represent the statistical mean
(blue) and standard deviation (black) of the projected simulated data
collected in bins of width $\Delta\psi=1/29$.

##  4 Conclusions

The results derived here provide a computationally attractive alternative to
and generalization of the earlier work by Bashkirtseva and Ryashko [7], Guo
_et al._ [23], and Zhao _et al._ [43]. Where these references rely on
projections onto hyperplanes orthogonal to the vector field, transversal to
the vector field but otherwise arbitrary, and stroboscopic sections in a case
of periodically excitations on top of an existing limit cycle motion, the
proposed boundary-value problem formulation is closely integrated with a
foliation of transversal hyperplanes that may be defined invariantly through a
defining boundary-value problem, without the need for an explicit coordinate
basis or additional approximation. From the invariance of this foliation under
the linearized deterministic flow, there further follows a closed-form
solution for the covariance function in the form of a rapidly convergent
series for autonomous as well as periodically excited systems.

A further advancement is the explicit recognition of the existence of
$(m+1)^{2}$ conserved quantities for the forward dynamics of the covariance
differential equation with unique periodic or quasiperiodic solutions in each
corresponding level surface. Following a recipe analogous to that developed by
Muñoz-Almaraz _et al._ [36], we removed this degeneracy by imposing a finite
set of transversality constraints and including symmetry-breaking parameters
collected in a matrix $A$ (scalar $a$ for the periodic orbit case) which must
equal $0$ on the sought solution. As described in conjunction with an
implementation of a suitable discretization of the covariance boundary-value
problem in the software package coco, the norm of $A$ provides a useful
measure of the truncation error. Notably, this is not available in our
alternative implementation that relies on the Moore-Penrose pseudo-inverse to
solve the overconstrained boundary-value problem with $A$ set to $0$ _a
priori_.

The availability of three examples for which analytical expressions may be
found for the adjoint variables and covariance matrix functions is another
useful contribution of this paper. As the authors can attest, such analytical
expressions are invaluable for debugging computational software that seeks to
make predictions about problems with uncertainty, especially for eliminating
various sources of systematic bias. In the present study, they highlighted an
early erroneous assumption about the symmetry of the matrices $B$ and $A$, and
led to the final imposition of the full set of $(m+1)^{2}$ transversality
conditions $\Omega(0)^{\mathsf{T}}C(0,0)\Omega(0)=0$.

A numerical solver for the covariance boundary-value problem was strictly only
needed for the second example in Section 3.2. Nevertheless, our coco
implementation relies on a general-purpose routine for one- (limit cycles) and
two-dimensional tori that is compatible with and may be appended to the
corresponding nonlinear continuation problem for the torus together with the
associated adjoint boundary-value problems. Although the final implementation
solves for the covariance function only after the torus and the adjoint
functions have been obtained, it is possible to embed the covariance boundary-
value problem with the full continuation problem. This comes at the expense of
a larger problem dimension but with the benefit that conditions may be imposed
on its solution. For example, it might be of interest to track families of
local maxima of the eigenvalues of the covariance function under variations in
problem parameters or to hold such local maxima fixed during continuation.

Finally, although the covariance boundary-value problem has been investigated
here as a means to characterize noise-induced statistics near transversally
stable periodic orbits of quasiperiodic invariant tori, we expect that it may
serve a useful, but different purpose also in the unstable, but normally
hyperbolic case. For such a periodic orbit or quasiperiodic torus, a solution
to the corresponding covariance boundary-value problem no longer captures the
statistics of trajectories that remain near the deterministic limit cycle over
long times, since none do. With suitably chosen noise, however, we expect that
such a solution could be used as an alternative means of quantifying the
stability of the deterministic object. We are particularly interested in
exploring this in the quasiperiodic case, as a substitute for analysis of the
solution to the variational problem.

## Acknowledgments

This material is based upon work performed while the second author (HD) was
supported by and served at the National Science Foundation. Any opinion,
findings, and conclusions or recommendations expressed in this material are
those of the authors and do not necessarily reflect the views of the National
Science Foundation. CK also wants to thank the VolkswagenStiftung for support
via a Lichtenberg Professorship.

##  Appendix A Adjoint Conditions and Problem Discretization

It is, of course, the rare exception when a closed form analysis is both
possible and expeditious. In all other cases, the boundary-value problems in
the main body of this paper must be implemented numerically in order to
compute $\gamma$, $\lambda$ or $\Lambda$, and $C$ and explore variations in
$C$ along continuous families of such solutions.

Of the three boundary-value problems, the one for the adjoint $\lambda$ and
$\Lambda$ stand out in that the differential equation and periodic or
quasiperiodic boundary conditions are supplemented by normalization conditions
involving integrals. We show here that the adjoint boundary-value problems,
including the integral conditions, may be obtained in stages from a calculus-
of-variations approach applied to a suitable constraint Lagrangian (as
discussed in [32]). We rely on this when using the built-in problem
constructors in coco to derive regular zero problems for discretized
approximations of $\gamma$ and $\lambda$ or $\Lambda$.

The following two propositions and subsequent corollary review results from
[15]. Original to this paper is Proposition A.6 which translates directly to
the coco implementation used to compute $\gamma$, $\Lambda$, and $C$ for the
second example in Sec. 3.2.

######  Proposition A.1.

Suppose that $h(\gamma(0))=0$ and $\partial_{x}h(\gamma(0))f(\gamma(0))\neq 0$
for some function $h$. Then, the boundary-value problem

$0=-\dot{\mu}^{\mathsf{T}}-T\mu^{\mathsf{T}}\mathrm{D}f(\gamma),\,0=\mu(0)-\mu(1),\,1=\int_{0}^{1}\mu^{\mathsf{T}}f(\gamma)\,\mathrm{d}t$
(138)

is equivalent to the adjoint conditions obtained by imposing vanishing
variations of

$L=S+\int_{0}^{1}\mu^{\mathsf{T}}\left(\dot{x}-Sf(x)\right)\,\mathrm{d}t+w^{\mathsf{T}}(x(0)-x(1))+\kappa
h(x(0))$ (139)

under variations in $x$ and $S$, when evaluated at $x=\gamma$ and $S=T$.

###### Proof.

The integral condition results directly from vanishing variations with respect
to $S$. Vanishing variations with respect to $x$ yield the differential
equation and the additional conditions

$0=\mu(1)^{\mathsf{T}}-w,\,0=-\mu(0)+w+\kappa\partial_{x}h(x(0))$ (140)

Multiplication of both equations by $f(\gamma(0))$ and use of periodicity and
the constancy of $\mu^{\mathsf{T}}f(\gamma)$ then yields $\kappa=0$ and the
claim follows. ∎

Recall the notation

$\nabla(\phi,t):=\begin{pmatrix}\partial_{t}\gamma(\phi,t)&\partial_{\phi_{1}}\gamma(\phi,t)&\cdots&\partial_{\phi_{m}}\gamma(\phi,t)\end{pmatrix}$
(141)

from the section on quasiperiodic orbits.

######  Proposition A.2.

Suppose that $h_{i}(\gamma(0,0))=0$, $i=1,\ldots,m+1$ and that the matrix

$\begin{pmatrix}\partial_{x}h_{1}(\gamma(0,0))\partial_{t}\gamma(0,0)&\cdots&\partial_{x}h_{m+1}(\gamma(0,0))\partial_{t}\gamma(0,0)\\\
\vdots&\ddots&\vdots\\\
\partial_{x}h_{1}(\gamma(0,0))\partial_{\phi_{m}}\gamma(0,0)&\cdots&\partial_{x}h_{m+1}(\gamma(0,0))\partial_{\phi_{m}}\gamma(0,0)\end{pmatrix}$
(142)

is invertible. Then, the boundary-value problem

$\displaystyle 0=-\dot{\mu}^{\mathsf{T}}$
$\displaystyle-T\mu^{\mathsf{T}}\mathrm{D}f(\gamma(\phi,t)),0=\mu(\phi+\rho,0)-\mu(\phi,1),$
$\displaystyle\begin{pmatrix}1&0&\cdots&0\end{pmatrix}=\int_{\mathbb{S}^{m}}\mu^{\mathsf{T}}\nabla\,\mathrm{d}\phi$
(143)

is equivalent to the adjoint conditions obtained by imposing vanishing
variations of

$\displaystyle L$ $\displaystyle=\ln
S+\int_{\mathbb{S}^{m}}\int_{0}^{1}\mu^{\mathsf{T}}\left(\dot{x}-Sf(x)\right)\,\mathrm{d}t\,\mathrm{d}\phi$
$\displaystyle\qquad+\int_{\mathbb{S}^{m}}w^{\mathsf{T}}(\phi)(x(\phi+\varrho,0)-x(\phi,1))\,\mathrm{d}\phi+\sum_{i=1}^{m+1}\kappa_{i}h_{i}(x(0,0))$
(144)

under variations in $x$, $S$, and $\varrho$, when evaluated at $x=\gamma$,
$S=T$, and $\varrho=\rho$.

###### Proof.

Vanishing variations with respect to $x$ yield the differential equation and
the additional conditions

$0=\mu^{\mathsf{T}}(\phi,1)-w^{\mathsf{T}}(\phi)$ (145)

and

$0=-\mu^{\mathsf{T}}(\phi,0)+w^{\mathsf{T}}(\phi-\rho)+\delta_{\mathrm{D}}(\phi)\sum_{i=1}^{m+1}\kappa_{i}\partial_{x}h_{i}(\gamma(0,0))$
(146)

in terms of the Dirac delta function $\delta_{\mathrm{D}}$. Multiplication of
the first condition from the right by $\nabla(\phi,1)$, of the second
condition from the right by $\nabla(\phi,0)$, addition of the results using
the constancy of $\mu^{\mathsf{T}}\nabla$ with respect to $t$ to eliminate
terms including $\mu$, and integration over $\mathbb{S}^{m}$ using
quasiperiodicity of $\nabla$ then yields $\kappa_{i}=0$, $i=1,\ldots,m+1$,
thus establishing the quasiperiodicity of $\mu$. The integral conditions now
follow from variations with respect to $S$ and $\varrho$. ∎

######  Corollary A.3.

By replacing the term $\ln S$ in the expression for $L$ by $-\varrho_{i}$, one
obtains the integral conditions associated with $\lambda_{\phi_{i}}$.

One notes, in the quasiperiodic case, the absence of partial derivatives with
respect to components of $\phi$ in either of the differential equations for
$\gamma$, $\Lambda$, or $C$. As described in [14], this suggest a
discretization of all $\phi$-dependent unknowns in terms of truncated Fourier
series such that, e.g.,

$\gamma(\phi,t)=a_{0}(\phi_{-1},t)+\sum_{k=1}^{N}a_{k}(\phi_{-1},t)\cos 2\pi
k\phi_{1}+b_{k}(\phi_{-1},t)\sin 2\pi k\phi_{1}$ (147)

for some integer $N$ and where $\phi_{-k}$ denotes the argument sequence
$\phi_{k+1},\ldots,\phi_{m}$.

######  Lemma A.4.

Let $\phi_{1,j}:=(j-1)/(2N+1)$ and define the notation

$B(\phi_{1}):=\begin{pmatrix}1&\cos 2\pi\phi_{1}&\sin 2\pi\phi_{1}&\cdots&\cos
2N\pi\phi_{1}&\sin 2N\pi\phi_{1}\end{pmatrix}$ (148)

and

$G(\phi_{-1},t):=\begin{pmatrix}\gamma_{1}(\phi_{-1},t)\\\ \vdots\\\
\gamma_{2N+1}(\phi_{-1},t)\end{pmatrix},$ (149)

where $\gamma_{j}(\phi_{-1},t):=\gamma(\phi_{1,j},\phi_{-1},t)$. Then, there
exists an invertible $(2N+1)\times(2N+1)$ matrix $\mathcal{F}$ such that

$\gamma(\phi,t)=\left(B(\phi_{1})\mathcal{F}\otimes
I_{n}\right)G(\phi_{-1},t)$ (150)

and

$\mathcal{F}^{\mathsf{T}}\begin{pmatrix}1&&&\\\ &1/2&&\\\ &&\ddots&\\\
&&&1/2\end{pmatrix}\mathcal{F}=\frac{1}{2N+1}I_{2N+1}.$ (151)

###### Proof.

By inspection, $\mathcal{F}^{-1}=$

${\small\begin{pmatrix}1&\cos 2\pi\phi_{1,1}&\sin 2\pi\phi_{1,1}&\cdots&\cos
2\pi N\phi_{1,1}&\sin 2\pi N\phi_{1,1}\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\ 1&\cos 2\pi\phi_{1,2N+1}&\sin
2\pi\phi_{1,2N+1}&\cdots&\cos 2\pi N\phi_{1,2N+1}&\sin 2\pi
N\phi_{1,2N+1}\end{pmatrix}}$ (152)

and the claim follows by properties of the trigonometric functions. ∎

######  Corollary A.5.

Let

$\mathcal{R}(\rho_{1}):=\begin{pmatrix}1&&&&&\\\ &\cos 2\pi\rho_{1}&\sin
2\pi\rho_{1}&&&\\\ &-\sin 2\pi\rho_{1}&\cos 2\pi\rho_{1}&&&\\\ &&&\ddots&&\\\
&&&&\cos 2\pi N\rho_{1}&\sin 2\pi N\rho_{1}\\\ &&&&-\sin 2\pi N\rho_{1}&\cos
2\pi N\rho_{1}\end{pmatrix}$ (153)

and

$\mathcal{D}:=2\pi\begin{pmatrix}0&&&&&\\\ &0&1&&&\\\ &-1&0&&&\\\
&&&\ddots&&\\\ &&&&0&N\\\ &&&&-N&0\end{pmatrix}$ (154)

Then, $B^{\prime}(\phi_{1})=B(\phi_{1})\mathcal{D}$,
$\mathcal{R}^{\prime}(\rho_{1})=\mathcal{R}(\rho_{1})\mathcal{D}$, and

$\gamma(\phi+\rho,t)=\left(B(\phi_{1})\mathcal{R}(\rho_{1})\mathcal{F}\otimes
I_{n}\right)G(\phi_{-1}+\rho_{-1},t),$ (155)

######  Proposition A.6.

Let

$M(\phi_{-1},t):=\begin{pmatrix}\mu(\phi_{1,1},\phi_{-1},t)\\\ \vdots\\\
\mu(\phi_{1,2N+1},\phi_{-1},t)\end{pmatrix}.$ (156)

The discretized form of the adjoint boundary-value problem in Proposition A.2
consists of the differential equations

$0=-\partial_{t}M^{\mathsf{T}}-TM^{\mathsf{T}}\begin{pmatrix}\mathrm{D}f(\gamma_{1})&&\\\
&\ddots&\\\ &&\mathrm{D}f(\gamma_{2N+1})\end{pmatrix}$ (157)

the boundary conditions

$0=\left(\mathcal{F}\otimes
I_{n}\right)M(\phi_{-1},1)-\left(\mathcal{R}(\rho_{1})\mathcal{F}\otimes
I_{n}\right)M(\phi_{-1}+\rho_{-1},0),$ (158)

and the integral conditions

$1=\frac{1}{2N+1}\int_{\mathbb{S}^{m-1}}M^{\mathsf{T}}\partial_{t}G\,\mathrm{d}\phi_{-1},$
(159)
$0=\frac{1}{2N+1}\int_{\mathbb{S}^{m-1}}M^{\mathsf{T}}\partial_{\phi_{i}}G\,\mathrm{d}\phi_{-1},$
(160)

for $i=2,\ldots,m$, and

$0=\int_{\mathbb{S}^{m-1}}M^{\mathsf{T}}\left(\mathcal{F}^{\mathsf{T}}\mathcal{D}^{\mathsf{T}}\begin{pmatrix}1&&&\\\
&1/2&&\\\ &&\ddots&\\\ &&&1/2\end{pmatrix}\mathcal{F}\otimes
I_{n}\right)G\,\mathrm{d}\phi_{-1},$ (161)

These are equivalent to the adjoint conditions obtained by imposing vanishing
variations of

$\displaystyle L$ $\displaystyle=(2N+1)\ln
S+\int_{\mathbb{S}^{m-1}}\int_{0}^{1}M^{\mathsf{T}}\left(\partial_{t}\Xi-
SF\right)\,\mathrm{d}t\,\mathrm{d}\phi_{-1}$
$\displaystyle+\int_{\mathbb{S}^{m-1}}w^{\mathsf{T}}(\phi_{-1})\bigg{(}\left(\mathcal{R}(\varrho_{1})\mathcal{F}\otimes
I_{n}\right)\Xi(\phi_{-1}+\varrho_{-1},0)-\left(\mathcal{F}\otimes
I_{n}\right)\Xi(\phi_{-1},1)\bigg{)}\,\mathrm{d}\phi_{-1}$
$\displaystyle\qquad\qquad\qquad+\sum_{i=1}^{m+1}\kappa_{i}h_{i}(x(0,0))$

where

$\Xi(\phi_{-1},t):=\begin{pmatrix}x(\phi_{1,1},\phi_{-1},t)\\\ \vdots\\\
x(\phi_{1,2N+1},\phi_{-1},t)\end{pmatrix},\,F(\phi_{-1},t):=\begin{pmatrix}f(x(\phi_{1,1},\phi_{-1},t))\\\
\vdots\\\ f(x(\phi_{1,2N+1},\phi_{-1},t))\end{pmatrix}$ (162)

under variations in $\Xi$, $S$, and $\varrho$, when evaluated at $\Xi=G$,
$S=T$, and $\varrho=\rho$.

###### Proof.

The form of the discretization is an immediate consequence of the results of
the preceding lemma and corollary. Variations of $L$ with respect to $\Xi$
yield the differential equations as well as the conditions

$\displaystyle 0$
$\displaystyle=M^{\mathsf{T}}(\phi_{-1},1)-w^{\mathsf{T}}(\phi_{-1})\left(\mathcal{F}\otimes
I_{n}\right)$ (163) $\displaystyle 0$
$\displaystyle=-M^{\mathsf{T}}(\phi_{-1},0)+w^{\mathsf{T}}(\phi_{-1}-\rho_{-1})\left(\mathcal{R}(\varrho_{1})\mathcal{F}\otimes
I_{n}\right)$
$\displaystyle\qquad+\delta_{\mathrm{D}}(\phi_{-1})\begin{pmatrix}1&0&\cdots&0\end{pmatrix}\otimes\sum_{i=1}^{m+1}\kappa_{i}\partial_{x}h_{i}(\gamma_{1}(0,0))$
(164)

in terms of the Dirac delta function. Invertibility of the matrix in (142)
again implies that $\kappa_{i}=0$ for $i=1,\ldots,m$, these yield the
discretized boundary conditions since
$\mathcal{F}\mathcal{F}^{\mathsf{T}}=\mathcal{R}(\rho_{1})\mathcal{F}\mathcal{F}^{\mathsf{T}}\mathcal{R}^{\mathsf{T}}(\rho_{1})$.
∎

######  Corollary A.7.

The discretized form of the integral conditions obtained from Corollary A.3
follows by replacing $\ln S$ in the expression for $L$ by $-\varrho_{i}$.

## References

  * [1] Zaid Ahsan, Christian Kuehn, and Harry Dankowicz.  A boundary-value problem for covariance analysis of stochastically perturbed limit cycles.  In Proceedings of the ASME 2022 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. Volume 9: 18th International Conference on Multibody Systems, Nonlinear Dynamics, and Control (MSNDC), 2022.
  * [2] S. Baars, J.P. Viebahn, T.E. Mulder, C. Kuehn, F.W. Wubs, and H.A. Dijkstra.  Continuation of probability density functions using a generalized Lyapunov approach.  Journal of Computational Physics, 336:627–643, 2017.
  * [3] I. Bashkirtseva.  Stochastic sensitivity analysis: Theory and numerical algorithms.  IOP Conference Series: Materials Science and Engineering, 192(1), 2017.
  * [4] I.A. Bashkirtseva and L.B. Ryashko.  Stochastic sensitivity of 3D-cycles.  Mathematics and Computers in Simulation, 66(1):55–67, 2004.
  * [5] Irina Bashkirtseva, Alexander B. Neiman, and Lev Ryashko.  Stochastic sensitivity analysis of noise-induced suppression of firing and giant variability of spiking in a Hodgkin-Huxley neuron model.  Phys. Rev. E, 91:052920, 2015.
  * [6] Irina Bashkirtseva and Lev Ryashko.  Stochastic sensitivity of the closed invariant curves for discrete-time systems.  Physica A: Statistical Mechanics and its Applications, 410:236–243, 2014.
  * [7] Irina Bashkirtseva and Lev Ryashko.  Sensitivity analysis of stochastically forced quasiperiodic self-oscillations.  Electronic Journal of Differential Equations, 2016(240):1–12, 2016\.
  * [8] Irina Bashkirtseva and Lev Ryashko.  Stochastic bifurcations, chaos and phantom attractors in the Langford system with tori.  International Journal of Bifurcation and Chaos, 30(16), 2020.
  * [9] Nils Berglund and Barbara Gentz.  Noise-induced phenomena in slow-fast dynamical systems: a sample-paths approach.  Springer Science & Business Media, 2006.
  * [10] Saranya Biswas, Aasifa Rounak, Przemysław Perlikowski, and Sayan Gupta.  Characterising stochastic fixed points and limit cycles for dynamical systems with additive noise.  Communications in Nonlinear Science and Numerical Simulation, 101:105870, 2021.
  * [11] Sergio Bittanti and Patrizio Colaneri.  Stochastic Periodic Systems, pages 291–310.  Springer London, London, 2009.
  * [12] Paolo Bolzern and Patrizio Colaneri.  The periodic Lyapunov equation.  SIAM Journal on Matrix Analysis and Applications, 9(4):499–512, 1988.
  * [13] Yu-Chen Cheng and Hong Qian.  Stochastic limit-cycle oscillations of a nonlinear system under random perturbations.  Journal of Statistical Physics, 182(3), 2021.
  * [14] Harry Dankowicz and Frank Schilder.  Recipes for Continuation.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2013\.
  * [15] Harry Dankowicz and Jan Sieber.  Sensitivity analysis for periodic orbits and quasiperiodic invariant tori using the adjoint method.  Journal of Computational Dynamics, 9(3):329–369, 2022.
  * [16] Eusebius J Doedel, Alan R Champneys, Fabio Dercole, Thomas F Fairgrieve, Yu A Kuznetsov, B Oldeman, RC Paffenroth, B Sandstede, XJ Wang, and CH Zhang.  AUTO-07p: Continuation and bifurcation software for ordinary differential equations.  https://github.com/auto-07p/auto-07p.  Accessed: 2023-08-08.
  * [17] Hua Fan, Xiuming Shan, Jian Yuan, and Yong Ren.  Covariances of linear stochastic differential equations for analyzing computer networks.  Tsinghua Science and Technology, 16(3):264–271, 2011.
  * [18] J. Galan-Vioque, F.J.M. Almaraz, and E.F. Macías.  Continuation of periodic orbits in symmetric Hamiltonian and conservative systems.  European Physical Journal: Special Topics, 223(13):2705–2722, 2014\.
  * [19] Jianzhong Gao, Juping Ji, Yanping Liu, Feng Zhang, Ruiwu Wang, and Hao Wang.  A confidence ellipse analysis for stochastic dynamics model of Alzheimer’s disease.  Nonlinear Dynamics, 111(14):13513–13533, 2023.
  * [20] Koushik Garain and Partha Sarathi Mandal.  Stochastic sensitivity analysis and early warning signals of critical transitions in a tri-stable prey-predator system with noise.  Chaos, 32(3), 2022.
  * [21] Giambattista Giacomin, Christophe Poquet, and Assaf Shapira.  Small noise and long time phase diffusion in stochastic limit cycle oscillators.  Journal of Differential Equations, 264(2):1019–1049, 2018.
  * [22] Kongming Guo, Jun Jiang, and Zigang Li.  Diffusion and persistence of rotor/stator synchronous full annular rub response under weak random perturbations.  Journal of Vibration Engineering and Technologies, 8(4):599–611, 2020.
  * [23] Kongming Guo, Jun Jiang, and Yalan Xu.  Approximation of stochastic quasi-periodic responses of limit cycles in non-equilibrium systems under periodic excitations and weak fluctuations.  Entropy, 19(6), 2017.
  * [24] A. Halanay, C. Tudor, T. Morozan, and A. Halanay.  Tracking almost periodic signals under white noise perturbations.  Stochastics, 21(4):287–301, 1987.
  * [25] Desmond Higham and Peter Kloeden.  An Introduction to the Numerical Simulation of Stochastic Differential Equations.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2021\.
  * [26] Mengya Huang, Anji Yang, Sanling Yuan, and Tonghua Zhang.  Stochastic sensitivity analysis and feedback control of noise-induced transitions in a predator-prey model with anti-predator behavior.  Mathematical Biosciences and Engineering, 20(2):4219–4242, 2023\.
  * [27] Jochen Jungeilges, Trygve Kastberg Nilssen, Tatyana Perevalova, and Alexander Satov.  Transitions between metastable long-run consumption behaviors in a stochastic peer-driven consumer network.  Discrete and Continuous Dynamical Systems - Series B, 26(11):5849–5871, 2021.
  * [28] Christian Kuehn.  Deterministic continuation of stochastic metastable equilibria via Lyapunov equations and ellipsoids.  SIAM Journal on Scientific Computing, 34(3):A1635–A1658, 2012.
  * [29] Christian Kuehn.  Efficient gluing of numerical continuation and a multiple solution method for elliptic PDEs.  Applied Mathematics and Computation, 266:656–674, 2015.
  * [30] Christian Kuehn.  Numerical continuation and SPDE stability for the 2D cubic-quintic Allen–Cahn equation.  SIAM/ASA Journal on Uncertainty Quantification, 3(1):762–789, 2015\.
  * [31] Christian Kurrer and Klaus Schulten.  Effect of noise and perturbations on limit cycle systems.  Physica D: Nonlinear Phenomena, 50(3):311–320, 1991.
  * [32] Mingwu Li and Harry Dankowicz.  Staged construction of adjoints for constrained optimization of integro-differential boundary-value problems.  SIAM Journal on Applied Dynamical Systems, 17(2):1117–1151, 2017\.
  * [33] Stilianos Louca.  Stable limit cycles perturbed by noise, 2018.
  * [34] Ian R. Manchester.  Transverse dynamics and regions of stability for nonlinear hybrid limit cycles.  IFAC Proceedings Volumes, 44(1):6285–6290, 2011.
  * [35] H. Mena and L. Pfurtscheller.  An efficient SPDE approach for El Niño.  Applied Mathematics and Computation, 352:146–156, 2019.
  * [36] F.J. Muñoz-Almaraz, E. Freire, J. Galán, E. Doedel, and A. Vanderbauwhede.  Continuation of periodic orbits in conservative and Hamiltonian systems.  Physica D: Nonlinear Phenomena, 181(1):1–38, 2003.
  * [37] L. Ryashko, I. Bashkirtseva, A. Gubkin, and P. Stikhin.  Confidence tori in the analysis of stochastic 3D-cycles.  Mathematics and Computers in Simulation, 80(2):256–269, 2009.
  * [38] Leonid A. Safonov and Yoshiharu Yamamoto.  Noise-driven switching between limit cycles and adaptability in a small-dimensional excitable network with balanced coupling.  Phys. Rev. E, 73:031914, 2006.
  * [39] S.M. Salman, A.M. Yousef, and A.A. Elsadany.  Dynamic behavior and bifurcation analysis of a deterministic and stochastic coupled logistic map system.  International Journal of Dynamics and Control, 10(1):69–85, 2022\.
  * [40] Frank Schilder, Hinke M. Osinga, and Werner Vogt.  Continuation of quasi-periodic invariant tori.  SIAM Journal on Applied Dynamical Systems, 4(3):459–488, 2005.
  * [41] Yukihiro Tadokoro, Hiroya Tanaka, and MI Dykman.  Noise-induced switching from a symmetry-protected shallow metastable state.  Scientific Reports, 10(1):1–10, 2020.
  * [42] Pascal Wang, Daniele Castellana, and Henk A. Dijkstra.  Improvements to the use of the trajectory-adaptive multilevel sampling algorithm for the study of rare events.  Nonlinear Processes in Geophysics, 28(1):135–151, 2021.
  * [43] Feng Zhao, Yang Li, and Xian-Bin Liu.  Noise induced escape from stable invariant tori.  Journal of Statistical Mechanics: Theory and Experiment, 2022(10), 2022.