Computational singular perturbation method for nonstandard slow-fast systems

TL;DR

This paper extends the computational singular perturbation (CSP) method to nonstandard slow-fast systems with normally hyperbolic attracting critical manifolds, providing new formulas and demonstrating its applicability on complex examples.

Contribution

It adapts the CSP method to nonstandard systems, offering new formulas and the first concrete demonstrations on genuinely nonstandard examples.

Findings

Successfully extended CSP to nonstandard systems

Derived new formulas for the adapted CSP method

Demonstrated applicability on complex nonstandard examples

Abstract

The computational singular perturbation (CSP) method is an algorithm which iteratively approximates slow manifolds and fast fibers in multiple-timescale dynamical systems. Since its inception due to Lam and Goussis, the convergence of the CSP method has been explored in depth; however, rigorous applications have been confined to the standard framework, where the separation between `slow' and `fast' variables is made explicit in the dynamical system. This paper adapts the CSP method to {\it nonstandard} slow-fast systems having a normally hyperbolic attracting critical manifold. We give new formulas for the CSP method in this more general context, and provide the first concrete demonstrations of the method on genuinely nonstandard examples.

Equations

$$z^{\prime }$$

$$S=\{z\in {\mathbb{R}}^n:NW(z)=0\Rightarrow W(z)=0\}.$$

$$\left(\begin{array}{c}{p}^{\prime }\\ d^{\prime }\\ r^{\prime }\\ b^{\prime }\end{array}\right)$$

$$S=\{(p,d,r,b)\in {\mathbb{R}}^4:(r,b)=\psi (p,d)\}$$

$$\psi (p,d)=\left(\frac{k_{-5}(-k_{4}d{p}^2+k_{-4}(c-d))}{k_4(k_5{p}^2+k_{-5})},\frac{k_5{p}^2(-k_{4}d{p}^2+k_{-4}(c-d))}{k_4(k_5{p}^2+k_{-5})}\right).$$

$$z^{\prime }=H(z,\epsilon )=N(z)f(z)+\epsilon G(z,\epsilon )$$

$$z^{\prime }=\frac{dz}{dt}=H(z,\epsilon )=N(z)f(z)+\epsilon G(z,\epsilon ),$$

$$S=\{z\in {\mathbb{R}}^n:f(z)=0\},$$

$$\dot{z}=\frac{dz}{d\tau }=\frac{1}{\epsilon }H(z,\epsilon )=\frac{1}{\epsilon }N(z)f(z)+G(z,\epsilon ),$$

$$z^{\prime }=H(z,0)=N(z)f(z).$$

$$T_z{\mathbb{R}}^n$$

$${\Pi }^S$$

$${\Pi }^N$$

$$DH(z,0)N(z){|}_S=(N(z)Df(z))N(z){|}_S=N(z)(Df(z)N(z)){|}_S.$$

$$\dot{z}$$

$${\Pi }^{\mathcal{V}}$$

$${\Pi }^{\mathcal{W}}$$

$${\Pi }^S$$

$$N$$

$$x^{\prime }$$

$$y^{\prime }$$

$${\varphi }_t({\mathcal{F}}^s(p))$$

$$||{\varphi }_t(p)-{\varphi }_t(q)||$$

$$y$$

$$z^{\prime }$$

$$H^{\prime }$$

$$\left(\begin{array}{c}{H}_x^{\prime }\\ H_y^{\prime }\end{array}\right)$$

$$g(z)$$

$$[X,H{]}_i$$

$$g^{\prime }$$

$$\Lambda (A,B,H)$$

$$\Lambda$$

Full text

1 Introduction

Most systems in nature consist of processes that evolve on disparate timescales and the observed dynamics in such systems reflect these multiple timescale features as well. Mathematical models of such multiple timescale systems are considered singular perturbation problems with slow-fast (or two timescale) problems as the most common. Models of homogeneously mixed biochemical reactions such as substrate-enzyme or ligand-receptor kinetics are prime examples.

An interesting and pervasive feature of these biochemical reaction systems is the observed transition from transient fast kinetics to long-term slow kinetics, wherein the system settles onto a so-called quasi-steady state (QSS). Geometrically, this QSS is perceived as a lower dimensional slow manifold. Identifying such an attracting lower dimensional slow manifold provides the means to reduce the dimension of biochemical reaction systems. Such QSS reduction techniques are frequently employed in the biochemical literature with Michaelis-Menten-type laws as prime examples; see e.g. [12].

The mathematical foundation to justify such a QSS reduction is given by Tikhonov’s [27] respectively Fenichel’s [4] work on normally hyperbolic attracting slow manifolds in singular perturbation problems. Goeke, Noethen and Walcher [22, 6] provide a comprehensive discussion on the general setup of Fenichel’s geometric singular perturbation theory (GSPT) with an emphasis on explaining when a QSS reduction is justified or when it leads to erroneous results.

Schauer and Heinrich [25] explored a homotopy approach to find QSS approximations in biochemical reaction networks, by formulating a perturbation problem in terms of a small parameter $\varepsilon>0$ amplifying the ratio of magnitudes of slow and fast reaction rates.111We note that a small perturbation parameter $\varepsilon>0$ could be properly identified via dimensional analysis. After using stoichiometry to derive the dynamical system from the network, the fast reactions are grouped into a vector $W$, and slow into a vector $V$:

[TABLE]

with state vector $z\in\mathbb{R}^{n}$, $N$ and $R$ are constant stoichiometric matrices with full column rank $m_{N},m_{R}<n$, and $\varepsilon>0$ is the homotopy parameter. With this splitting, a $k$-dimensional manifold of stationary states (QSS), $1\leq k<n$, can be identified in the singular limit $\varepsilon\to 0$,

[TABLE]

We emphasize here that the slow-fast splitting (1) is nonstandard from the point of view of GSPT in the sense that slow-fast reactions are distinguished rather than slow-fast variables. Nevertheless, the theory of Fenichel [4] still applies: under the appropriate geometric condition (see Section 2 for details), an invariant slow manifold $S_{\varepsilon}$ perturbs from $S$ for sufficiently small $\varepsilon>0$.

As a motivating example, consider the following four-dimensional slime mold cell communication model [5]:

[TABLE]

Here, $p$ refers to the concentration of cAMP; $d$ and $r$ represent transmembrane receptors; and $b$ refers to the bound state of $r$. The parameters $k_{i}\geq 0$ are constant reaction rates, the parameter $c\geq 0$ represents a conserved quantity that arises from the model reduction of the original five-dimensional model due to Stiefenhofer [26], and the parameter $S\geq 0$ denotes the concentration of ATP which is assumed to be constant. This system contains a two-dimensional manifold of equilibria given as a graph

[TABLE]

with

[TABLE]

When $\varepsilon>0$ is sufficiently small, numerical observations reveal the decay of fast transient motion as trajectories are attracted to a low-dimensional invariant set $S_{\varepsilon}$ near $S$, as shown in Figure 1. Our objective is to approximate such invariant sets arbitrarily well.

Several algorithmic & computational techniques have been developed to identify these QSS (i.e. a lower dimensional slow manifold), and to identify the simplified system underlying long-term stability. Lam and Goussis [17, 18] devised the computational singular perturbation (CSP) method for chemical kinetics modelling which is an iterative procedure that generates a sequence of CSP manifolds and CSP fibers which approximate the slow manifolds and fast fibers of the dynamical system; see Section 3 for details. The convergence of these CSP objects has been analyzed by Kaper, Kaper, and Zagaris in a series of papers [30, 31, 32, 11]. They rigorously prove this convergence in the case of slow-fast vector fields satisfying a spectral gap condition for the eigenvalues of the Jacobian of the vector field along an attracting invariant manifold, but focus on examples in the standard form. Indeed, the practical applicability of the CSP method has remained unclear in the case where we cannot easily identify the slow or fast components of the vector field a priori. Lam and Goussis [18] give an example of CSP applied to a linearized problem by populating the initial conditions with (constant) eigenvectors, but pointed out the necessity of a modification for more difficult nonlinear problems where the Jacobian of the vector field is point-dependent.

Our approach is to avoid this defect for nonstandard slow-fast systems of the form

[TABLE]

which includes system (3). The idea is that the singular limit vector field factorization of (6), i.e. the vector field $H(z,0)=N(z)f(z)$, encodes a wealth of geometric information. We will make explicit use of this factorization to provide formulas for a nonstandard version of the CSP method.

Although this paper is concerned with the CSP method, there are in fact a variety of iterative schemes to approximate invariant manifolds; the zero-derivative principle (ZDP), intrinsic low-dimensional manifolds (ILDM), and center manifold normal-form reductions [24] are among a few of these (chapter 11 in [15] give an overview of the first three methods). Many of these schemes are interconnected; for instance, the CSP and ZDP iterations differ only by a multiplicative factor [32].

The paper is organised as follows: In Sec. 2, we describe a general framework for nonstandard slow-fast systems. In Sec. 3 we define the two-step CSP update. Our main results are provided in Secs. 4 and 5: we give specific formulas for initializing the CSP method in the nonstandard context, study the output of the first update carefully, and then give examples demonstrating these formulas. We conclude in Sec. 6 by highlighting some fruitful new connections between the CSP method and the factorization given in the second section. We also provide an appendix (App. A), which expands on the framework underlying the CSP update step.

2 Nonstandard slow-fast dynamical systems

We begin by giving an abbreviated treatment of a general framework for nonstandard slow-fast systems. Much of this material in fact appears in Fenichel’s seminal work on GSPT [4]; see also [22, 6]. This approach has been further developed by Wechselberger [28] and extends the framework to loss of normal hyperbolicity.

We are interested in two-timescale (or slow-fast) dynamical systems of the form (6), which we restate here for convenience:

[TABLE]

with state variable $z\in\mathbb{R}^{n}$, $n\times(n-k)$ matrix $N(z)$ formed by column vectors $N^{i}(z)=(N^{i}_{1}(z),\ldots,N^{i}_{n}(z))^{\top}$ with sufficiently smooth functions $N^{i}:\mathbb{R}^{n}\to\mathbb{R}$, $i=1,\ldots,n-k$, $f(z)=(f_{1}(z),\ldots,f_{n-k}(z))^{\top}$ a column vector of sufficiently smooth functions $f_{i}:\mathbb{R}^{n}\to\mathbb{R}$, $i=1,\ldots,n-k$, $G(z,\varepsilon)=(G_{1}(z,\varepsilon),\ldots,G_{n}(z,\varepsilon))^{\top}$ a column vector of sufficiently smooth functions $G_{i}:\mathbb{R}^{n}\to\mathbb{R}$, $i=1,\ldots,n$, and $\varepsilon\ll 1$ characterizes the ratio of timescales in the system.

Definition 2.1

*Let $S_{0}$ denote the set of equilibria of system (7) in the singular limit $\varepsilon\to 0$. If there exists a subset $S\subseteq S_{0}$ which forms a $k$-dimensional differentiable manifold of equilibria with $1\leq k<n$, then system (7) defines a singular perturbation problem.

Assumption 2.1

System (7) is a singular perturbation problem with a single subset $S\subseteq S_{0}$,

[TABLE]

*which forms a $k$-dimensional differentiable manifold of equilibria, $1\leq k<n$, called the critical manifold.

Assumption 2.2

*In system (7), the matrix $N(z)$ has full (column) rank for all $z\in S$.

Next consider system (7), rescaled from the fast timescale $t$ to the slow timescale $\tau=\varepsilon t$:

[TABLE]

Systems (7) and (9) are equivalent when $\varepsilon>0$ but their singular limits $\varepsilon\to 0$ are not. In fact, they carry complementary, lower dimensional information, and it is a cornerstone of GSPT to concatenate the information from these two limiting problems to deduce the dynamics of the full system (7) respectively (9).

2.1 The layer problem

We focus first on system (7) evolving on the fast time scale $t$.

Definition 2.2

The layer problem of system (7) is given by the formal limit $\varepsilon\to 0$:

[TABLE]

Under Assumption 2.1, the set $S$ forms a $k$-dimensional manifold of equilibria of the layer problem. Hence, the Jacobian $Dh|_{S}$ along $S$ has $k$ trivial zero eigenvalues and $(n-k)$ nontrivial eigenvalues.

Definition 2.3

*A $k$-dimensional critical manifold $S$ is called normally hyperbolic if the $(n-k)$ nontrivial eigenvalues of the Jacobian $Dh|_{S}=(ND\!f)|_{S}$ are bounded away from the imaginary axis.

The existence of a normally hyperbolic manifold implies the following splitting:

[TABLE]

where $T_{z}S$ denotes the tangent space of the critical manifold $S$ at $z$, coinciding with the kernel of the linear map $DH(z,0)$, and $\mathcal{N}_{z}$ is the unique complement of the splitting identified with the quotient space $T_{z}\mathbb{R}^{n}/T_{z}S$. We call $\mathcal{N}_{z}$ the linear fast fiber with basepoint $z\in S$. Repeating this construction across all points of $S$, we obtain the tangent bundle $TS=\cup_{z\in S}T_{z}S$ of $S$, and the transverse linear fast fiber bundle $\mathcal{N}=\cup_{z\in S}\mathcal{N}_{z}$. This splitting induces unique projection operators onto the tangent bundle $TS$ and the linear fast fibre bundle $\mathcal{N}$ as follows:

[TABLE]

Lemma 1

*At any point $z\in S$ of a normally hyperbolic manifold $S$, the rows of $Df(z)$ form a basis of $T_{z}S^{\bot}$ and the columns of $N(z)$ form a basis of $\mathcal{N}_{z}$.

Proof. The manifold $S$ is the zero level set of $f(z)$ which forms a $k$-dimensional manifold, i.e. $D\!f(z)|_{S}$ has full row rank $(n-k)$ for any $z\in S$. Thus the rows of $Df(z)$ form a basis of the orthogonal complement of $T_{z}S$.

Since $S$ is normally hyperbolic, we have $DH(z,0)$ has rank $(n-k)$ which implies $\ker DH(z,0)=T_{z}S$ at points $z\in S$. Under Assumption 2.2, $N(z)$ has full column rank $(n-k)$. Thus the column spaces of $DH(z,0)$ and $N(z)$ coincide at points $z\in S$ and, hence, the columns of $N(z)$ form a basis of $\mathcal{N}_{z}$. $\Box$

Lemma 2

*Let $S$ be a normally hyperbolic manifold. Then the $(n-k)\times(n-k)$ square matrix $D\!fN|_{S}$ is regular, and its eigenvalues are equal to the set of nontrivial eigenvalues of $Dh|_{S}$.

Proof. Let the linear operator $DH(z,0)|_{S}$ act on the basis of the invariant subset $\mathcal{N}$, i.e.

[TABLE]

By Lemma 1, the $(n-k)\times(n-k)$ square matrix $D\!fN|_{S}$ is necessarily a regular matrix with $(n-k)$ nonzero eigenvalues which coincide with the nontrivial eigenvalues of $DH(z,0)|_{S}$. Note, if $p(z)\in\mathbb{R}^{n-k}$ is an eigenvector of $D\!fN|_{S}$ with eigenvalue $\mu(z)$ then $N(z)p(z)\in\mathbb{R}^{n}$ is the corresponding eigenvector of $DH(z,0)$ with the same eigenvalue. $\Box$

Assumption 2.3

*The $k$-dimensional critical manifold $S$ of system (7) is normally hyperbolic and attracting, i.e. all $(n-k)$ nontrivial eigenvalues have negative real part.

Remark 1

System (3) satisfies Assumptions 2.1–2.3. A rich source of examples comes from chemical reaction networks; Schauer and Heinrich [25] provide a derivation of such models with the goal of identifying quasi-steady state approximations in biochemical reaction networks.

*Nevertheless, it is unclear whether general techniques exist to compute such factorizations in a typical slow-fast system. In the case of polynomial vector fields, Goeke and Walcher [6] have used division algorithms for algebraic varieties to factor $H(z,0)$.

Remark 2

*In this paper we do not concern ourselves with singularities of $H(z,0)$ outside of the critical manifold $S$; we only require the local geometric structure provided by the differentiable manifold and transverse fibers in the proceeding arguments. A simple example of a system having an isolated singularity as well as a critical manifold is given in Sec. 5.1.

2.2 The reduced problem

Now consider system (9) evolving on the slow timescale $\tau=\varepsilon t$. The singular limit $\varepsilon\to 0$ of system (9) requires more care, i.e. this limit will only be well-defined provided that we restrict the phase space to $S$ and that $G(z,0)$ is restricted to the tangent bundle $TS$. The projection operator $\Pi^{S}$ (12) allows us to formulate this limit.

Definition 2.4

The reduced problem of system (9) is

[TABLE]

Definition 2.5

Consider the splitting $\mathbb{R}^{n}=\mathcal{V}\oplus\mathcal{W}$, where $\mathcal{V}$ has dimension $n-k$ and $\mathcal{W}$ has dimension $k$. The oblique projection of a vector $p=x+y$, $x\in\mathcal{V}$ and $y\in\mathcal{W}$, onto the subspace $\mathcal{V}$ parallel to $\mathcal{W}$ is a linear map $\Pi^{\mathcal{V}}=(\Pi^{\mathcal{V}})^{2}$ satisfying $\Pi^{\mathcal{V}}(p)=x$.

Suppose $V$ (resp. $U$) is an $n\times(n-k)$ matrix whose column vectors span $\mathcal{V}$ (resp. $\mathcal{W}^{\bot}$). Then it can be shown that $\Pi^{\mathcal{V}}$ has the following matrix representation:

[TABLE]

The (complementary) oblique projection onto the subspace $\mathcal{W}$ parallel to $\mathcal{V}$ is given by

[TABLE]

In the present context, the splitting $T\mathbb{R}^{n}=\mathcal{N}\oplus TS$ induces an oblique projection onto the tangent bundle of the critical manifold, parallel to the fast fibers (see Fig. 2). Lemma 1 gives us a matrix representation using the matrices $N$ and $Df$:

[TABLE]

By Lemma 2, the matrix $D\!fN$ is regular and, hence, the inverse is well defined. Equivalent projection formulas appear in [4] and [6].

Remark 3

It is easily shown that standard singular perturbation problems are a special case of a nonstandard problem (7). Indeed, select

[TABLE]

and describe the variables and vector field in terms of components $z=(x,y)\in\mathbb{R}^{k}\times\mathbb{R}^{(n-k)}$ and $G(x,y,\varepsilon)=(g(x,y,\varepsilon),\tilde{f}(x,y,\varepsilon))$. Then $z^{\prime}=N(z)f(z)+\varepsilon G(z,\varepsilon)$ gives

[TABLE]

which is a standard singular perturbation problem. From (19) we observe that $x$ lists the $k$ slow variables and $y$ lists the $(n-k)$ fast variables. These systems lie in contrast to the nonstandard case, where one cannot expect such a trivial factorization (18) to exist globally.

2.3 Fenichel Theory

In the case of normally hyperbolic critical manifolds, QSS reductions onto a slow invariant manifold are justified by the following theorem.

Theorem 1** (Fenichel’s Theorem [4, 10, 15])**

Given the system (7) with a $C^{r}$-smooth vector field and a compact normally hyperbolic critical manifold, the following hold for $\varepsilon>0$ sufficiently small:

• •

There exists a locally invariant $C^{r}$-smooth, normally hyperbolic slow manifold $S_{\varepsilon}$ that is $C^{r}$ $\mathcal{O}(\varepsilon)$-close to $S$.

• •

The flow on $S_{\varepsilon}$ converges to the reduced flow on $S$ as $\varepsilon\to 0$.

• •

There are $C^{r}$-smooth locally invariant stable and unstable manifolds, $\mathcal{F}^{s}(S_{\varepsilon})$ and $\mathcal{F}^{u}(S_{\varepsilon})$.

• •

These manifolds admit nonlinear, $C^{r-1}$-smooth foliations $\{\mathcal{F}^{s}(p):p\in S_{\varepsilon}\}$ resp. $\{\mathcal{F}^{u}(p):p\in S_{\varepsilon}\}$. Furthermore, these families are positively (resp. negatively) invariant on fibers; i.e. if $\phi_{t}$ is the time-$t\geq 0$* flow map of (7) and if $p,\phi_{t}(p)\in S_{\varepsilon}$, then*

[TABLE]

An analogous statement holds for the fibers $\mathcal{F}^{u}(p)$, with $p\in S_{\varepsilon}$.

• •

There exist constants $C,\lambda>0$ such that if $q\in\mathcal{F}^{s}(p)$, then for $t\geq 0$ we have

[TABLE]

Analogous rate estimates hold for the fibers $\mathcal{F}^{u}(p)$, with $p\in S_{\varepsilon}$.

The invariant slow manifold $S_{\varepsilon}$ and its invariant nonlinear foliation $\mathcal{F}_{\varepsilon}$ (which split into the sub-foliations given in the theorem above) organise the dynamics of the full system: trajectories perturb from concatenated orbits of the layer and reduced problems under the appropriate transversality conditions. In analogy to the $\varepsilon=0$ setting, we denote the linear fast fiber at basepoint $p\in S_{\varepsilon}$ by $\mathcal{N}_{\varepsilon,p}$ and the linear fast fiber bundle by $\mathcal{N}_{\varepsilon}=\cup_{p\in S_{\varepsilon}}\mathcal{N}_{\varepsilon,p}$.

3 The CSP iteration

Suppose that the invariant slow manifold $\mathcal{S}_{\varepsilon}$ of a nonstandard slow-fast system (7) is (locally) given by the graph of a function $y=h_{\varepsilon}(x)$ with $h_{\varepsilon}:\mathbb{R}^{k}\to\mathbb{R}^{n-k}$. We express $h_{\varepsilon}(x)$ as an asymptotic series in $\varepsilon$:

[TABLE]

The terms $h_{j}(x)$ may be obtained reinserting this series into (7) and matching coefficients in orders of $\varepsilon$. The CSP method adopts a different, iterative approach to efficiently compute not only the terms $h_{j}(x)$, but also a similar asymptotic approximation to the linear fast fibers transverse to the slow manifold. The motivating idea is to understand how the variational equation of a dynamical system is affected by changes in basis.

3.1 Geometric framework of the CSP method

Consider a smooth vector field $H(z)$, to which we append the variational equation to produce the dynamical system on the tangent bundle:

[TABLE]

Under some arbitrarily chosen splitting $z=(x,y)$, where $x\in\mathbb{R}^{n-k}$ and $y\in\mathbb{R}^{k}$, we can write the variational equation in block component form:

[TABLE]

Now suppose $A(z)=[A_{f}(z)~{}A_{s}(z)]$ is a smooth, regular $n\times n$ matrix for all $z\in S$ with $B(z)=\begin{pmatrix}B_{s\bot}(z)\\ B_{f\bot}(z)\end{pmatrix}$ its dual. Here, the first and second block columns of $A$, denoted $A_{f}$ and $A_{s}$, are of sizes $n\times(n-k)$ and $n\times k$, respectively; similarly, the first and second block rows $B_{s\bot}$ and $B_{f\bot}$ have sizes $(n-k)\times n$ and $k\times n$, respectively. The vector field $H(z)$ expressed in this new basis is given by

[TABLE]

Definition 3.1

For a smooth $n\times m$ matrix function $X$ and a smooth vector field $H$, $[X,H]$ is denoted the Lie bracket of $X$ and $H$ and defined as the $n\times m$ matrix whose $i$th column is

[TABLE]

*where $i=1,\cdots,m$ and $A_{i}$ refers to the $i$th column of $A$.

It can be shown that the variational equation of the transformed vector field can be written compactly in the form

[TABLE]

where the nonlinear operator $\Lambda(A,B,H)$ has the algebraic structure of a Lie bracket:

[TABLE]

see the appendix (Sec. A.1) for a derivation of this result. The operator $\Lambda$ can be written in block-component form:

[TABLE]

The key insight is that in the presence of a $k$-dimensional invariant manifold $\mathcal{M}$ and corresponding transverse linear fiber bundle $\mathcal{N}$, the vector field $H$ may be expressed in a clever choice of basis $A(z)$ and dual $B(z)$ which block-diagonalizes $\Lambda$. As a consequence, the invariant manifold $\mathcal{M}$ and the linear fiber bundle $\mathcal{N}$ are easily characterized in terms of this basis as follows:

[TABLE]

see the appendix (Sec. A.2) for details. This formalism applies directly to the case of an invariant slow manifold $\mathcal{M}=S_{\varepsilon}$ and its accompanying linear fast fiber bundle $\mathcal{N}=\mathcal{N}_{\varepsilon}$, as defined by Fenichel’s theorem (Sec. 2.3). The characterizations (27)–(28) then state that $S_{\varepsilon}$ is defined by the locus of points $B_{s\bot}(z)H(z)=0$ where the components of the vector field $H(z)$ lying in the direction of the linear fast fibers vanish, and that the linear fast fibers $\mathcal{N}_{p}$ at basepoints $p\in S_{\varepsilon}$ are spanned by the columns of the block component $A_{f}(p)$.

3.2 CSP objects

The CSP iteration acts on a suitably initialised (point-dependent) basis matrix $A^{(0)}(z)$ and dual $B^{(0)}(z)$, producing a sequence $\{(A^{(j)}(z),B^{(j)}(z))\}_{j=0}^{\infty}$ of successively refined bases. We initialise with the pair

[TABLE]

The sequence $\{(A^{(j)}(z),B^{(j)}(z))\}_{j=0}^{\infty}$ in turn defines a sequence of updates $\{\Lambda^{(j)}\}_{j=0}^{\infty}$ to the $\Lambda$ operator:

[TABLE]

These approximate operators may be written in block components in analogy to (26). Motivated by the characterizations of the invariant manifold and fast fibers provided in Eqs. (27)–(28), we define the following approximating objects:

Definition 3.2

The CSP manifold of order [math] is the level set

[TABLE]

For integers $j\geq 1$, the CSP manifold of order $j$, denoted $\mathcal{K}^{(j)}$, is

[TABLE]

*where $y=\psi^{(j-1)}(x,\varepsilon)$ is a graph of $\mathcal{K}^{(j-1)}$.  

Definition 3.3

For $j\geq 0$ and $p\in\mathcal{K}^{(j)}$, the CSP fiber of order $j$ is the subspace

[TABLE]

The CSP fiber bundle of order $j$ is the corresponding vector bundle

[TABLE]

Remark 4

*The convergence of the CSP manifolds and fiber bundles to the invariant slow manifold and fast fiber bundle of system (7) depends on the choice of iteration step, as shown in Sec. 3.4.

3.3 One-step and two-step CSP updates

There are two commonly-used variants of the CSP iteration. We will only apply the two-step method in this paper, but it is instructive to introduce the simpler one-step method to clarify the relationship between these iterations and the CSP objects (32)–(33).

Both methods use near-identity transformations in the update step, i.e. multiplication by matrices of the form $I\pm U$ and $I\pm L$, where $I$ is the identity matrix, and $U$ and $L$ are nilpotent matrices of the form

[TABLE]

The block components $\tilde{U}^{(j)}$ and $\tilde{L}^{(j)}$, of respective sizes $(n-k)\times k$ and $k\times(n-k)$, are defined by the constraint that the CPS manifolds and fibers converge asymptotically to the slow manifold and linear fast fibers; see Sec. 3.4. Near-identity update matrices are computationally easy to invert. In particular, we obtain efficient update rules for the dual basis $B^{(j)}$, as shown in the following two methods.

3.3.1 One-step CSP method

The one-step CSP method is given by the iteration rule

[TABLE]

The one-step CSP method only updates ‘half’ of the basis. This update is sufficient if we are interested in computing only the CSP manifolds (32). If we wish to approximate the fast fibers to the invariant manifold in tandem, we must update the remaining blocks of the basis and dual basis matrices. This update is provided by the two-step CSP method.

3.3.2 Two-step CSP method

The two-step CSP method is given by the iteration rule

[TABLE]

3.4 Convergence

Suppose that after $j$ iterates of the CSP two-step method, the CSP manifold of order $j$ defined by (32) is locally expressed as a graph $y=\psi^{(j)}(x,\varepsilon)$ and then expanded as an asymptotic series in $\varepsilon$:

[TABLE]

When the update terms in the near-identity transformations (3.3.2) are defined appropriately, Kaper, Kaper, and Zagaris [30] demonstrated convergence of the CSP manifold to the invariant slow manifold of Fenichel’s theory (as described in Sec. 2.3), in the following sense.

Theorem 2** (Convergence of CSP manifolds [30])**

Suppose the leading-order term of the graph of $\mathcal{K}^{(0)}$, denoted $y=\psi^{(0)}(x,\varepsilon)$, agrees with the graph of the critical manifold $y=h_{0}(x)$ under a local choice of coordinates. Then using the one-step or two-step CSP update rule and for $j$ fixed, we have

[TABLE]

when $\varepsilon>0$ is sufficiently small, where $\psi^{(j)}$ is the graph of $\mathcal{K}^{(j)}$ with asymptotic series expansion (38) and $y=h(x)=\sum_{j=0}^{\infty}h_{j}(x)$ is the asymptotic series expansion of the graph of $S_{\varepsilon}$.

This theorem states if the CSP two-step method is initialised appropriately, then for pairs of indices $0\leq i\leq j$ we have $\psi^{(j)}_{i}(x)=h_{i}(x)$. This justifies the nomenclature ‘CSP manifold of order $j$’: the set $\mathcal{K}^{(j)}$ agrees with $S_{\varepsilon}$ up to order $j$ terms in the asymptotic expansion.

An analogous convergence statement can be given for the CSP fiber bundle (34):

Theorem 3** (Convergence of CSP fibers [31])**

Fix $j\geq 0$ and let $\mathcal{K}^{(j)}$ be given locally by the graph $y=\psi^{(j)}(x)$ for $x\in\mathbb{R}^{k}$ and $y\in\mathbb{R}^{n-k}$. Then given the basepoint $(x,\psi^{(j)}(x))$ on $\mathcal{K}^{(j)}$ (resp. $(x,h_{\varepsilon}(x))$ on $S_{\varepsilon}$), the asymptotic expansions of $\mathcal{L}_{\varepsilon}^{(j)}(x,\psi^{(j)}(x))$ and $\mathcal{N}_{\varepsilon}(x,h(x))$ agree up to and including terms of $\mathcal{O}(\varepsilon^{j})$. Thus, the CSP fast fiber bundle $\mathcal{L}^{(j)}$ is an $\mathcal{O}(\varepsilon^{j})$ approximation to the fast fiber bundle $\mathcal{N}_{\varepsilon}$.

Theorems 2–3 hold when the near-identity updates in (35) are defined as follows:

[TABLE]

With this choice of update step, the CSP two-step method simultaneously block-diagonalizes the CSP operator $\Lambda^{(j)}$ (as defined in (30)) in discrete steps as follows [31]:

Lemma 3

For $j=0,1,\cdots$, the CSP two-step method (3.3.2) provides the asymptotic estimates

[TABLE]

where $\Lambda^{(j)}$ is evaluated on $\mathcal{K}^{(j)}$.

4 Nonstandard CSP updates

We now describe the CSP update step for nonstandard slow-fast systems (7) satisfying Assumptions 2.1–2.3. The key is that bases for the fast and slow subspaces can be ‘read off’ using the factorization $H(z,0)=N(z)f(z)$, providing a natural initial condition for the iteration. By Lemma 1, the columns of $N(z)$ span the linear fast fiber with basepoint $z\in S$, and the columns of $Df(z)^{\top}$ span $T_{z}S^{\bot}$. Thus,

[TABLE]

where the notation $P^{\bot}$ refers to a matrix whose columns form a basis to the subspace orthogonal to the column space of the matrix $P$. We define

[TABLE]

where the matrix prefactors in both block rows (which are well-defined by Lemma 1) normalize the block diagonal components in the product $B^{(0)}A^{(0)}$ to the identity, as required. With this initialization, the leading-order approximation $\mathcal{K}^{(0)}$ of the CSP manifold defined in (31) can be computed:

[TABLE]

This level set implicitly defines the graph $y=\psi^{(0)}(x,\varepsilon)$. On this graph, we expand both sides of this equation in powers of $\varepsilon$. The leading-order $O(\varepsilon^{0})$ coefficient is

[TABLE]

which matches the definition of the critical manifold $S$ as the leading-order part of the asymptotic series of $\mathcal{S}_{\varepsilon}$, $f(x,h_{0}(x),0)=0$; thus, $\psi^{(0)}_{0}(x)=h_{0}(x)$.

The four block components of the operator $\Lambda^{(0)}=B^{(0)}[A^{(0)},H]$ can be computed explicitly in terms of the components $N,f,$ and $G$ in the vector field $H$; compare with (26):

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The vector function $f$ is $\mathcal{O}(\varepsilon)$ on $\mathcal{K}^{(0)}$ by (42); therefore, the latter three block components are of at most $\mathcal{O}(\varepsilon)$ as well. This in turn implies that the update blocks $\tilde{U}^{(0)}$ and $\tilde{L}^{(0)}$ defined in (39) are both $\mathcal{O}(\varepsilon)$ on this set. After one application of the CSP two-step update (3.3.2), the updated bases $A^{(1)}$ and $B^{(1)}$ are

[TABLE]

Without writing out the update in detail, it suffices to note that the block components $\tilde{L}^{(0)}$ and $\tilde{U}^{(0)}$ each introduce $\mathcal{O}(\varepsilon^{1})$ perturbations of the initial block columns of $A^{(0)}$ (respectively block rows of $B^{(0)}$). The first fast fiber update is explored in greater detail using projection maps, in Sec. 6.1. Subsequent CSP iterations will result in higher-order perturbations.

The first update of the CSP manifold $\mathcal{K}^{(1)}$ is defined as follows; see (32):

[TABLE]

In this equation and the following, we use the additional subscript ’[math]’ to refer to those quantities that are computed on points $(x,\psi^{(0)}(x))\in\mathcal{K}^{(0)}$. We compute the order $O(\varepsilon^{1})$ coefficient of the asymptotic series assuming the graph form $y=\psi^{(1)}(x,\varepsilon)$ where the leading order coefficient is defined by $f_{0}=0\,$; see (43):

[TABLE]

which gives

[TABLE]

This result matches a recent result due to Wechselberger [28] for the first-order correction of $\mathcal{S}_{\varepsilon}$ in the nonstandard case.

5 Examples

We now demonstrate the two-step CSP method in four nonstandard examples. Our list of systems increases in complexity: first we consider the case where the slow manifold is equal to the critical manifold and the fast fiber bundle remains unchanged. Second, we consider the case where the slow manifold is equal to the critical manifold but the fast fiber bundle perturbs from the $\varepsilon=0$ case. Third, we consider a system where both the manifold and fiber updates are nontrivial. Finally, we consider the four-dimensional system (3) discussed in the introduction.

5.1 Trivial updates of the slow manifold and fast fibers

Consider the following two-dimensional nonstandard system:

[TABLE]

The critical manifold $S=\{x^{2}+y^{2}=1\}$ is a circle. The Jacobian evaluated along $S$ is

[TABLE]

with one trivial zero eigenvalue and one negative eigenvalue $\lambda=-2$ for all $(x,y)\in S$, implying that the critical manifold $S$ is attracting and normally hyperbolic. Thus, Assumptions 2.1–2.3 are satisfied. Note that the origin is an isolated singularity of $N(z)$, where $f(z)\neq 0$, but we do not consider this singularity in the CSP iteration (see Remark 2).

When $\varepsilon$ is small and positive, there is a clear separation between slow motion along $S$ versus fast motion toward $S$. The dynamics is made obvious if we write the system in polar coordinates and plot a graph of representative trajectories (Fig. 3):

[TABLE]

The coordinate representation (55) decomposes the system into two independent ODEs. The invariant manifold $r=1$, which is independent of $\varepsilon$, can be read off from the first equation of (55). Thus, $S_{\varepsilon}=S$. We can also determine that the (nonlinear) fast fiber bundle $\mathcal{F}_{\varepsilon}$ is foliated by the family of rays extending from the origin, with basepoints on the circle:

[TABLE]

Note that the nonlinear and linear fast fiber bundles coincide in this example: $\mathcal{F}_{\varepsilon}=\mathcal{N}_{\varepsilon}$. The fast fiber bundle is invariant as a family: for $t\in\mathbb{R}$, the fiber $\mathcal{F}_{\varepsilon}(p(\theta_{0}))$ is mapped to $\mathcal{F}_{\varepsilon}(p(\theta_{0}+\varepsilon t))$ under the time $t$ flow map. Furthermore, exponential contraction of the rays toward $S_{\varepsilon}$ is governed by the equation $r^{\prime}=r(1-r^{2})$. Compare these facts about the nonlinear fast fibers with Fenichel’s theorem, Sec. 2.3. This behavior is similarly independent of the value of $\varepsilon>0$.

We initialize the CSP two-step method using (40):

[TABLE]

We begin by computing the initial CSP manifold $\mathcal{K}^{(0)}$ (see (31)). We compute

[TABLE]

and so

[TABLE]

We have

[TABLE]

With this information, the initial $\Lambda^{(0)}$ matrix (30) becomes

[TABLE]

Using (39), we find that $\tilde{U}^{(0)}=0$ and $\tilde{L}^{(0)}=0$.

The triviality of the updates implies that the CSP manifolds and fibers ((32)–(33)) are

[TABLE]

for $j=0,1,2,\cdots$.

5.2 Trivial updates of the slow manifold; nontrivial updates of the fast fibers

We now consider a variant of the previous system, where we modify the first component of $G(x,y)$:

[TABLE]

We still have $S=\{(x,y):x^{2}+y^{2}=1\}$, and it is easy to check that this is also the invariant slow manifold for $\varepsilon>0$: for $p=(x,y)\in S$ we have

[TABLE]

so that $H(p,\varepsilon)\in T_{p}S$. Thus, $S=S_{\varepsilon}$ as before. In this modified system, however, the linear fast fiber $\mathcal{N}_{\varepsilon,p}$ will not be orthogonal to $T_{p}S_{\varepsilon}$ at points $p\in S_{\varepsilon}$.

The CSP step is initialized as usual with (40). The basis matrices are identical to (57)–(58) since we have not modified $N(z)$ or $f(z)$. From (30), the initial operator $\Lambda^{(0)}$ is (compare (60))

[TABLE]

By definition (39), the update term $\tilde{U}^{(0)}$ is trivial since $\Lambda_{fs}^{(0)}=0$; on the other hand, the update quantity $\tilde{L}^{(0)}$ may give a nontrivial first-order correction to the fast fiber. The first column (resp. first row) of the updated basis matrix $A^{(1)}$ (resp. $B^{(1)}$) are

[TABLE]

Applying the definitions (32)–(33), the updated CSP objects are

[TABLE]

Here, the function $y(x)$ refers to a graph of $x^{2}+y^{2}=1$ containing the chosen basepoint $(x,y(x))$. Graphs of the form $x(y)$ can also be used.

5.3 Parabolic critical manifold

Given the system

[TABLE]

The critical manifold of the system is defined as the level set $f=0$, which in this case can be solved globally for the graph $y=1-x^{2}$. This critical manifold is globally attracting and normally hyperbolic—the Jacobian along $S$ has one trivial zero eigenvalue and another eigenvalue $-(3x^{2}+1)<0$.

The CSP method is initialized with the following basis and dual as given in (40):

[TABLE]

The initial CSP manifold (31) written to $\mathcal{O}(\varepsilon^{0})$ order is

[TABLE]

Define the auxiliary function $g(x)=3x^{2}+1$. The operator $\Lambda^{(0)}$ (as defined in (30)) computed to linear order is

[TABLE]

After one application of the two-step CSP method (3.3.2), we obtain the updated basis matrix

[TABLE]

and dual $B^{(1)}|_{\mathcal{K}^{(0)}}$.

This gives us the first update of the CSP objects using the definitions (32)–(33):

[TABLE]

The $\Lambda^{(1)}$ update (30) is given by

[TABLE]

Observe that the off-diagonal elements are now 0 (modulo nonzero $\mathcal{O}(\varepsilon^{2})$ terms) while the diagonal terms are unchanged up to order $\varepsilon$. As we continue to apply the CSP algorithm, the off-diagonal elements can be made to vanish modulo arbitrarily high orders; see Lemma 3. The order we used in the initial basis matrix implies that the (1,1) term governs the fast dynamics and (2,2) term governs the slow dynamics in the decoupled system.

5.4 A four-dimensional system

Recall system (3):

[TABLE]

for parameters $k_{i},S,c>0$, together with the perturbation parameter $\varepsilon>0$. This system has the factorized form $H(z)=N(z)f(z)+\varepsilon G(z,\varepsilon)$. Using (40), the initial bases are given by

[TABLE]

and the dual by $B^{(0)}=(A^{(0)})^{-1}$, whose first two rows are given by the $2\times 4$ matrix $(DfN)^{-1}Df$.

As before, the leading-order term of the graph of $\mathcal{K}^{(0)}$ (defined in (31)) will match the graph of $\{f=0\}$:

[TABLE]

Applying the CSP manifold definition (32), we have

[TABLE]

where $\psi^{(1)}_{0}(p,d)=\psi^{(0)}_{0}(p,d)$ and the first-order correction $\psi^{(1)}_{1}(p,d)$ can be factored as

[TABLE]

for the seven polynomials $\beta_{i},\delta_{i},\gamma$ in $p$, defined as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The fast fiber update, which we do not write down, will be an $\mathcal{O}(\varepsilon)$ correction to the first $4\times 2$ block of $A^{(0)}$, in analogy to the previous examples.

Remark 5

*The tensorial nature of the update step only becomes apparent in systems with dimension higher than two. For example, to compute $\Lambda_{fs}^{(0)}$, we must compute the tensor $DA_{s}^{(0)}$. Such computations quickly become unwieldy for high-dimensional systems after the first or second CSP iterate.

Remark 6

Symbolic computations for these examples were performed using Mathematica Version 11.1.1.0 [29].

6 Conclusion

We have given a detailed treatment of the CSP method for nonstandard slow-fast systems (7), when the vector field admits a geometrically intuitive factorization in the leading-order term, and demonstrated the method explicitly for nontrivial examples. There appear to be several other natural connections between the CSP algorithm and this factorization, which we comment on further.

6.1 CSP as projection

There is a tantalizing connection between the CSP two-step update and the projection operators (12)–(13) which were naturally induced by the normally hyperbolic splitting. Consider the CSP fiber of order 1, $\mathcal{L}^{(1)}=\text{Col}(A_{f}^{(1)})$, defined by (3.3.2). We have

[TABLE]

where the term $A_{s}^{(0)}\tilde{L}^{(0)}$ can be written out explicitly by using the formulas (4) and (4) for the block components of $\Lambda$ in the update formula (39):

[TABLE]

where

[TABLE]

A careful comparison to the matrix representation of oblique projection (Definition 2.5) reveals that in fact $\Pi^{\bot}=\Pi^{S}$ to leading order, when we evaluate this expression on the set $\mathcal{K}^{(0)}$. Since $f=O(\varepsilon)$ on $\mathcal{K}^{(0)}$ as shown in (42), the equation (70) becomes

[TABLE]

where the subscript 0 reminds the reader that these quantities are computed on $\mathcal{K}^{(0)}$, and the vector $\xi_{0}$ is rigorously defined by the two-step CSP rule (3.3.2). The CSP fiber $\mathcal{L}^{(1)}$ is defined on its corresponding CSP manifold $\mathcal{K}^{(1)}$ by Def. (33), but it is known that the fiber approximation is still $\mathcal{O}(\varepsilon^{j})$-accurate if $A_{f}^{(j)}$ is evaluated on $\mathcal{K}^{(j-1)}$ (see Sec. 3.5 in [31]).

In the examples 5.1 and 5.2, $\Pi^{S}$ is an orthogonal projection onto $T_{z}S$, and furthermore the manifold $S=S_{\varepsilon}$ is given by $f=0$ exactly. The computation (71) shows that the fast fiber update is trivial if $[N,G]\in\text{ker~{}}\Pi^{S}$. This is indeed the case in the first example 5.1 but not in the modification 5.2.

These results demonstrate the close relationship between the CSP update step and the projectors which define the fast and slow subsystems of (7). Deeper connections, including a recasting of the CSP method as a projective scheme in which the projections are themselves iteratively updated, remain to be further clarified.

6.2 CSP as a factorization algorithm for nonstandard vector fields

As discussed in Sec. 2, the leading-order factorization gives a compact description of the geometry of the system near the critical manifold $S$ when $\varepsilon=0$, and gives approximate information about the dynamics on the nearby slow manifold $S_{\varepsilon}$ when $\varepsilon>0$. Consider the following conjecture.

Conjecture 1

Vector fields satisfying Assumptions 2.1–2.3 may be re-factorized as follows:

[TABLE]

The new objects $N_{\varepsilon}(z,\varepsilon)$, $f_{\varepsilon}(z,\varepsilon)$, and $\tilde{G}(z,\varepsilon)$ satisfy the following properties:

• •

(Slow manifold as a level set) $S_{\varepsilon}=\{z\in\mathbb{R}^{n}:f_{\varepsilon}(z,\varepsilon)=0\}$.

• •

(Basis of fast fibers) The columns of $N_{\varepsilon}(p,\varepsilon)$ span the linear fast fibers $\mathcal{N}_{\varepsilon}(p)$ at basepoints $p\in S_{\varepsilon}$.

• •

(Invariance) $\tilde{G}(p,\varepsilon)\in T_{p}S_{\varepsilon}$ for $p\in S_{\varepsilon}$.

The assertion is that a vector field factorization with respect to the critical manifold $S$ when $\varepsilon=0$ will induce a new factorization with respect to the slow manifold $S_{\varepsilon}$ when $\varepsilon>0$.

We describe a possible direction to answer this conjecture in the affirmative in the special case where $S_{\varepsilon}=S$ (as in the examples 5.1 and 5.2). The idea is to use the CSP method to obtain iterative refinements of the objects $N_{\varepsilon}(z,\varepsilon)$, $f_{\varepsilon}(z,\varepsilon)$, and $\tilde{G}(z,\varepsilon)$. Writing out the first few terms in the asymptotic series of the CSP fast fiber $\mathcal{L}^{(1)}$ (see (33) and (3.3.2)), we have $A^{(1)}_{f}(z,\varepsilon)=N(z)+\varepsilon A^{(1)}_{f,1}(z)+\mathcal{O}(\varepsilon^{2})$.

The vector field can be re-factorized as follows:

[TABLE]

where

[TABLE]

In comparison to the original factorization $H=Nf+\varepsilon G$, this new factorization has $\varepsilon$-dependence in the ‘leading-order’ term $N^{(1)}f^{(1)}$, and a modified ‘remainder’ term $G^{(1)}$. But observe that $S=S_{\varepsilon}=\{f^{(1)}=0\}$. Furthermore, the column space of the prefactor matrix $N^{(1)}$ is now an $\mathcal{O}(\varepsilon)$-approximation of the fast fibers. Finally, $G^{(1)}(z,\varepsilon)|_{S}=G(z,\varepsilon)\in T_{p}S$ since $S$ is invariant by assumption.

We can extend this idea in the obvious way to refine the factorization to some arbitrary order $j$:

[TABLE]

We remind the reader that the triviality of the level set update $f^{(j)}=f$ is a consequence of the assumption that $S=S_{\varepsilon}$. The general case where the CSP manifolds also update nontrivially is a topic of further study. Here, more care must be taken to write down the factorizations since the fiber basepoints are not fixed. Furthermore, nontrivial cross-terms of the form $N^{(k)}f^{(l)}$, where $k+l>j$, begin to appear at a given iterate $j$.

6.3 CSP for slow manifolds of saddle-type

Normally hyperbolic slow manifolds of saddle-type arise in traveling-wave profiles of the FitzHugh-Nagumo model [1], in the Hodgkin-Huxley model [8], and in models of cardiac pacemakers [14]. They are very challenging to numerically approximate even in the standard case, due to exponential instabilities in both forward and backward time [7]. It is therefore of interest to relax the assumption that the critical manifold be attracting for the CSP method.

Appendix A The formalism underlying the CSP method

The purpose of this section is to collect a few results given across the three very detailed papers [30, 31, 11].

A.1 Change of variables and Lie bracket

In this section, we justify the statement that the variational equation of the transformed vector field has the structure of a Lie bracket as shown in (24)–(25). Differentiating $BA=I$ with respect to $t$, we obtain the identity

[TABLE]

Furthermore, the chain rule gives

[TABLE]

Using identities (72)-(A.1), the variational equation of the transformed vector field $f$ becomes:

[TABLE]

This calculation appears in [11].

A.2 Block-diagonalization of $\Lambda$

In this section, we justify the characterizations of the invariant manifold and transverse fiber bundle in terms of a cleverly chosen basis, as shown in (27)–(28). In the presence of an invariant manifold having a transverse fiber bundle which is invariant as a family (i.e. the flow maps fibers into fibers), $H$ can be expressed in a basis which block-diagonalizes the operator $\Lambda$. This was proven in [30, 31] for the case of a slow manifold having a transverse fast fiber bundle, but the identical argument carries over without distinguishing ‘slow’ versus ‘fast’ components. We demonstrate this for the vanishing of the upper-diagonal block.

Let $\mathcal{M}$ be an invariant manifold with a transverse fiber bundle $\mathcal{N}$; i.e., at a point $z\in\mathcal{M}$, we have the following splitting:

[TABLE]

We write $A=[A_{f}~{}A_{s}]$, where the columns of the $n\times(n-k)$ matrix $A_{F}(z)$ form a basis of $\mathcal{N}_{z}$ and the columns of the $n\times k$ matrix $A_{s}(z)$ form a basis of $T_{z}M$. The corresponding dual basis similarly spans the dual splitting: we have $B=\begin{pmatrix}B_{s\bot}\\ B_{f\bot}\end{pmatrix}$.

Proposition 1

For the choice of basis $A$ above, we have $\Lambda_{fs}=\mathbb{O}_{n-k,k}$.

Proof. We have

[TABLE]

By invariance, $H\in TM$ for points in $\mathcal{M}$, and so $B_{s\bot}H=0$. The directional derivative along $A_{s}$ must therefore be identically 0 on $\mathcal{M}$ as well. This gives us the following identity:

[TABLE]

Here we clarify that $DB_{s\bot}$ is a symmetric bilinear form, and this (matrix) identity is to be understood as taking the standard directional derivative of each column of $A_{s}$ in turn and concatenating the result into a matrix.

Similarly, we have the trivial identity $B_{s\bot}A_{s}=0$ on $\mathcal{M}$, coming from the dual basis criterion. Differentiating with respect to $t$ and using the chain rule, we obtain

[TABLE]

Subtracting identities (76) from (A.2) and using the symmetry of the bilinear form, we obtain $\Lambda_{fs}=0$. $\Box$

The argument for $\Lambda_{sf}=\mathbb{O}_{n-k,k}$ is similar if slightly more involved, using the invariance of the fiber bundle to obtain an identity as above [31].