Gait modeling and optimization for the perturbed Stokes regime

Abstract

Many forms of locomotion, both natural and artificial, are dominated by viscous friction in the sense that without power expenditure they quickly come to a standstill. From geometric mechanics, it is known that for swimming at the “Stokesian” (viscous; zero Reynolds number) limit, the motion is governed by a reduced-order “connection” model that describes how body shape change produces motion for the body frame with respect to the world. In the “perturbed Stokes regime” where inertial forces are still dominated by viscosity, but are not negligible (low Reynolds number), we show that motion is still governed by a functional relationship between shape velocity and body velocity, but this function is no longer linear in shape change rate. We derive this model using results from singular perturbation theory and the theory of noncompact normally hyperbolic invariant manifolds. Using the theoretical properties of this reduced-order model, we develop an algorithm that estimates an approximation to the dynamics near a cyclic body shape change (a “gait”) directly from observational data of shape and body motion. This extends our previous work which assumed kinematic “connection” models. To compare the old and new algorithms, we analyze simulated swimmers over a range of inertia-to-damping ratios. Our new class of models performs well on the Stokesian regime and over several orders of magnitude outside it into the perturbed Stokes regime, where it gives significantly improved prediction accuracy compared to previous work. In addition to algorithmic improvements, we thereby present a new class of models that is of independent interest. Their application to data-driven modeling improves our ability to study the optimality of animal gaits and our ability to use hardware-in-the-loop optimization to produce gaits for robots.

Appendices

Appendix A: Derivation of the equations of motion

In this and the following section, we consider systems more general than those considered earlier and in so doing assume that the reader is familiar with some basic concepts in geometric mechanics and differential geometry: Lie groups, group actions, and principal bundles. We refer the reader to Kobayashi and Nomizu [26], Marsden and Ratiu [31], Lee [27], and Bloch [3] for the relevant standard definitions related to Lie groups and group actions, and we refer the reader to Kobayashi and Nomizu [26], Marsden et al. [28], Marsden [29], and Bloch [3] for material on bundles.

We consider a mechanical system on a configuration space Q whose Lagrangian is of the form kinetic minus potential energy. We will also consider this system to be subjected to external viscous forcing arising from a Rayleigh dissipation function and also subjected to an external force exerted by the locomoting body. We are interested in the situation that we have a smooth action $\theta :G \times Q \rightarrow Q$ of a Lie group G on Q, such that the Lagrangian, viscous forces, and external force are all symmetric under the action. In this case, we say that G is a symmetry group.

In “Appendix A.1,” we will define some geometric quantities on Q which encode information about the symmetry and the dynamics. Working in coordinates induced by a local trivialization, in “Appendix A.2” we derive the equations of motion in terms of these quantities. In “Appendix A.3,” we recall how the equations become governed by the so-called viscous connection in the Stokesian limit [7, 25], which will set the stage for our derivation in “Appendix B” of a corrected reduced-order model for the perturbed Stokes regime.

1.1 The mechanical and viscous connections

In this section, we define the mechanical and viscous (or Stokes) connections, roughly following Kelly and Murray [25]. We consider a Lagrangian $L:\mathsf {T}Q \rightarrow \mathbb {R}$ which is invariant under the lifted action $\mathsf {D}\theta _g$ of G on $\mathsf {T}Q$ (here $\mathsf {D}$ denotes the derivative or pushforward). We assume the Lagrangian to be of the form kinetic minus potential energy, where kinetic energy is given by $\frac{m}{2}k$, where $m > 0$ is a dimensionless mass parameter, k is a smooth symmetric bilinear form, and mk is the kinetic energy metric. In what follows, we assume that k is positive definite when restricted to tangent spaces to G orbits, but not necessarily that k is positive definite on all tangent vectors.¹¹ Denoting by $\mathfrak {g}$ the Lie algebra of G and $\mathfrak {g}^*$ its dual, we define the (Lagrangian) momentum map$J:\mathsf {T}Q \rightarrow \mathfrak {g}^*$ via

$$\begin{aligned} \langle J(v_q), \xi \rangle = \langle \mathbb {F}L(v_q), \xi _Q(q) \rangle = mk_q(v_q, \xi _Q(q)),\nonumber \\ \end{aligned}$$

(14)

where $v \in \mathsf {T}_q Q$ and $\xi \in \mathfrak {g}$. Here $\mathbb {F}L:\mathsf {T}Q\rightarrow \mathsf {T}^* Q$ is the fiber derivative of L given by $\mathbb {F}L(v_q)(w_q):=\frac{\partial }{\partial s}|_{s=0} L(v_q +$s$w_q)$, and the smooth vector field $\xi _Q$ on Q is the infinitesimal generator defined by $\xi _Q(q):=\frac{\partial }{\partial s}|_{s=0}\theta _{\exp (s\xi )}(q)$. We define the mechanical connection$\Gamma _{\text {mech}}:\mathsf {T}Q \rightarrow \mathfrak {g}$ via $\Gamma _{\text {mech}}(v_q):=\mathbb {I}^{-1}(q)J(v_q)$, where $\mathbb {I}(q):\mathfrak {g}\rightarrow \mathfrak {g}^*$ is the locked inertia tensor defined via

$$\begin{aligned} \langle \mathbb {I}(q) \xi , \eta \rangle:= & {} \langle \mathbb {F}L(\xi _Q(q)),\eta _Q(q)\rangle \nonumber \\= & {} mk_q(\xi _Q(q),\eta _Q(q)), \end{aligned}$$

(15)

where $\xi ,\eta \in \mathfrak {g}$.

We now follow an analogous procedure to define the viscous connection $\Gamma _{\text {visc}}:\mathsf {T}Q \rightarrow \mathbb {R}$. We consider a Rayleigh dissipation function $R :\mathsf {T}Q \rightarrow \mathbb {R}$ defined in terms of a G-invariant smooth symmetric bilinear form $\nu$ on Q: $R(v_q):=\frac{c}{2}\nu _q(v_q,v_q)$, where $c>0$ is a dimensionless parameter representing the amount of damping or dissipation in the system due to viscous forces. As with k, we assume that $\nu$ is positive definite when restricted to tangent spaces to G orbits, but not necessarily that $\nu$ is positive definite on all tangent vectors.¹² The corresponding force field $F_R:\mathsf {T}Q \rightarrow \mathsf {T}^* Q$ is given by minus the fiber derivative of R, $F_R:=\mathbb {F}(-R)$. We define a map $K:\mathsf {T}Q \rightarrow \mathfrak {g}^*$, analogous to the momentum map J, via

$$\begin{aligned} \langle K(v_q), \xi \rangle = \langle F_R(v_q), \xi _Q(q) \rangle = -c\nu _q(v_q, \xi _Q(q)),\nonumber \\ \end{aligned}$$

(16)

where $v \in \mathsf {T}_q Q$ and $\xi \in \mathfrak {g}$. We define the viscous connection or Stokes connection$\Gamma _{\text {visc}}:\mathsf {T}Q \rightarrow \mathfrak {g}$ via $\Gamma _{\text {visc}}(v_q):=\mathbb {V}^{-1}(q)K(v_q)$, where $\mathbb {V}(q):\mathfrak {g}\rightarrow \mathfrak {g}^*$ is defined via

$$\begin{aligned} \langle \mathbb {V}(q) \xi , \eta \rangle:= & {} \langle F_R(\xi _Q(q)),\eta _Q(q)\rangle \nonumber \\= & {} -c\nu _q(\xi _Q(q),\eta _Q(q)), \end{aligned}$$

(17)

where $\xi ,\eta \in \mathfrak {g}$.

Using the G-invariance of L and $\nu$, a calculation shows that $\Gamma _{\text {mech}}$ and $\Gamma _{\text {visc}}$ are equivariant with respect to the adjoint action of G on $\mathfrak {g}$:

$$\begin{aligned}&\forall g \in G:\Gamma _{\text {mech}}\circ \mathsf {D}\theta _g = \text {Ad}_g \circ \Gamma _{\text {mech}}, \quad \Gamma _{\text {visc}}\circ \mathsf {D}\theta _g\nonumber \\&\quad = \text {Ad}_g \circ \Gamma _{\text {visc}}\end{aligned}$$

(18)

Hence, if the natural projection $\pi _Q:Q \rightarrow Q/G$ from Q to the space of orbits Q / G of points in Q is a principal G-bundle, then the mechanical and viscous connections $\Gamma _{\text {mech}}$ and $\Gamma _{\text {visc}}$ are indeed principal connections; this justifies their titles.

Now in order for our system to move itself through space, we also allow there to be a G-equivariant external force $F_E:\mathbb {R}\times \mathsf {T}Q\rightarrow \mathsf {T}^*Q$ exerted by the locomoting body, subject to the requirement that $F_E$ takes values in the annihilator of $\ker \mathsf {D}\pi _Q$, the distribution tangent to group orbits. This requirement reflects the physically reasonable assumption that the locomoting body can exert only “internal forces” which directly affect only its shape $r\in Q/G$ (c.f. Eldering and Jacobs [7, Sec. 3.3] and Bloch et al. [2, Sec. 4.2]). For future use, we now prove the following

Proposition 1

The derivative of J along trajectories of the G-symmetric mechanical system is given by

$$\begin{aligned} \dot{J} = K, \end{aligned}$$

(19)

making the canonical identifications $\mathsf {T}_J \mathfrak {g}\cong \mathfrak {g}$.

Proof

We compute in a local trivialization on $\mathsf {T}Q$ induced by a chart for Q, so that we may write a trajectory as $(q,\dot{q})$. Note that in such local coordinates, $\mathbb {F}L(q,\dot{q})(v_q) = \frac{\partial L(q,\dot{q})}{\partial \dot{q}}v_q$. Hence

$$\begin{aligned} \begin{aligned} \langle \dot{J}(q,\dot{q}), \xi \rangle&= \frac{\mathrm{d}}{\mathrm{d}t}\left( \frac{\partial L(q(t),\dot{q}(t))}{\partial \dot{q}}\xi _Q(q(t))\right) \\&= \left( \frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial L}{\partial \dot{q}}\right) \xi _Q(q) + \frac{\partial L}{\partial \dot{q}}\mathsf {D}\xi _Q(q)\dot{q}\\&= \left( \frac{\partial L}{\partial q} + F_R + F_E \right) \xi _Q(q)\\&\quad + \frac{\partial L}{\partial \dot{q}}\mathsf {D}\xi _Q(q)\dot{q}, \end{aligned} \end{aligned}$$

(20)

where we obtained the last line using $\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial L}{\partial \dot{q}} - \frac{\partial L}{\partial q}=F_R + F_E,$ which follows from the Lagrange–d’Alembert principle [3, p. 8]. Since $F_E$ annihilates tangent vectors to group orbits, $\langle F_E, \xi _Q(q)\rangle = 0$. Hence, rearranging and letting $\Phi _\xi ^s$ denote the flow of $\xi _Q$, we find

$$\begin{aligned} \begin{aligned} \langle \dot{J}(q,\dot{q}), \xi \rangle =&\frac{\partial }{\partial s} L\left( \Phi _\xi ^s(q(t)), \mathsf {D}\Phi _\xi ^s(q(t))\dot{q}(t)\right) \\&+ \langle F_R(q,\dot{q}), \xi _Q(q)\rangle \\ =&\frac{\partial }{\partial s} L\left( \Phi _\xi ^s(q(t)), \mathsf {D}\Phi _\xi ^s(q(t))\dot{q}(t)\right) \\&+ \langle K(q,\dot{q}), \xi \rangle . \end{aligned} \end{aligned}$$

The derivative term is zero due to the invariance of L under the action of G, so from the arbitrariness of $\xi \in \mathfrak {g}$ we obtain the desired result. $\square$

As a corollary, we obtain a slight generalization of the classical Noether’s theorem.

Corollary 1

(Noether’s theorem) Consider a mechanical system given by a G-invariant Lagrangian of the form kinetic minus potential energy. Assume that the only external forces take values in the annihilator of the distribution tangent to the G orbits. Then the derivative of the momentum map J along trajectories satisfies

$$\begin{aligned} \dot{J}=0. \end{aligned}$$

Proof

Set $K = 0$ in Proposition 1. $\square$

1.2 Local form of the equations of motion

Assuming that the action of G on Q is free and proper [27, Ch. 21] so that $\pi _Q:Q\rightarrow Q/G$ is a principal G-bundle, we now derive the equations in a local trivialization, following [25]. In a local trivialization $U \times G$, $\pi _Q$ simply becomes projection onto the first factor and the G action is given by left multiplication on the second factor. We define $S :=Q/G$ to be the shape space representing all possible shapes of a locomoting body, and we write a point in the local trivialization as $(r,g)\in U\times G$ where $U\subset S$. We assume that U is the domain of a chart for S, so that we have induced coordinates $(r,\dot{r})$ for $\mathsf {T}U$.

Defining the body velocity¹³$\mathring{g} :=\mathsf {D}\mathrm {L}_{g^{-1}}\dot{g}$, the equivariance property (18) of the connection forms $\Gamma _{\text {mech}}, \Gamma _{\text {visc}}$ implies that they may be written in the trivialization as

$$\begin{aligned} \begin{aligned} \Gamma _{\text {mech}}(r,g)\cdot (\dot{r},\dot{g})&= \text {Ad}_g\left( \mathring{g}+A_{\text {mech}}(r)\cdot \dot{r}\right) \\ \Gamma _{\text {visc}}(r,g)\cdot (\dot{r},\dot{g})&= \text {Ad}_g\left( \mathring{g}+A_{\text {visc}}(r)\cdot \dot{r}\right) , \end{aligned} \end{aligned}$$

(21)

where $A_{\text {mech}}:\mathsf {T}U \rightarrow \mathfrak {g}$ and $A_{\text {visc}}:\mathsf {T}U \rightarrow \mathfrak {g}$ are, respectively, the local mechanical connection and local viscous connection. We define a diffeomorphism $(r,\dot{r},g,\dot{g})\mapsto (r,\dot{r},g,p)$, with p the body momentum defined by

$$\begin{aligned} p:=\text {Ad}_g^*{J} \in \mathfrak {g}^*. \end{aligned}$$

(22)

Here $\text {Ad}_g^*$ is the dual of the adjoint action $\text {Ad}_g$ of G on $\mathfrak {g}$. We additionally define

$$\begin{aligned} \begin{aligned} \mathbb {I}_{\text {loc}}&:=\text {Ad}_g^* \mathbb {I}\text {Ad}_g :\mathfrak {g}\rightarrow \mathfrak {g}^*\\ \mathbb {V}_{\text {loc}}&:=\text {Ad}_g^* \mathbb {V}\text {Ad}_g :\mathfrak {g}\rightarrow \mathfrak {g}^* \end{aligned} \end{aligned}$$

(23)

to be the local forms of $\mathbb {I}$ and $\mathbb {V}$. We note that the invariance of the Lagrangian L and Rayleigh dissipation function R under G, together with the general identity $\mathsf {D}\theta _g \xi _{Q}(q) = (\text {Ad}_g\xi )_{Q}(\theta _{g}(q))$, imply that $\mathbb {I}_{\text {loc}}(r),\mathbb {V}_{\text {loc}}(r)$ depend on the shape variable r only.

Rearranging (21), using the expressions (22), (23), and using Proposition 1, we obtain the equations of motion

$$\begin{aligned} \begin{aligned} \mathring{g}&= -A_{\text {mech}}\cdot \dot{r} + \mathbb {I}_{\text {loc}}^{-1}p\\ \dot{p}&= \mathbb {V}_{\text {loc}}(A_{\text {visc}}- A_{\text {mech}})\cdot \dot{r} + \mathbb {V}_{\text {loc}}\mathbb {I}_{\text {loc}}^{-1} p\\&\quad + \text {ad}^*_{\mathbb {I}_{\text {loc}}^{-1}p}p-\text {ad}^*_{A_{\text {mech}}\cdot \dot{r}}p, \end{aligned} \end{aligned}$$

(24)

where we have suppressed the r-dependence of $A_{\text {mech}},A_{\text {visc}},\mathbb {I}_{\text {loc}},\mathbb {V}_{\text {loc}}$ for readability. Notice that the $\dot{p}$ equation is completely decoupled from g.

In this paper, we are interested in the effect of shape changes on body motion, and not on the generation of shape changes themselves. Hence, we have suppressed the equations for $\dot{r},\ddot{r}$ from (24), simply viewing $r, \dot{r}$ as inputs in those equations, but see Bloch et al. [2] for more details on the specific form of the equations.

We merely note that, if the kinetic energy metric is positive definite, then the Lagrangian is hyperregular and our assumption of G-equivariance of the exerted force $F_E$ implies that

$$\begin{aligned} \ddot{r} = f(t,r,\dot{r},\mathbb {I}_{\text {loc}}^{-1} p) \end{aligned}$$

(25)

for some function f which depends on the local trivialization. If the kinetic energy metric is not positive definite (for use in toy examples like those in Sect. 4; see the precise assumptions in “Appendix A.1,” and the footnote there), then we assume that $\ddot{r}$ is given by (25).

1.3 Reduction in the Stokesian limit

From the definitions (15), (17) of $\mathbb {I}_{\text {loc}}, \mathbb {V}_{\text {loc}}$, we see that we may define $\bar{\mathbb {I}}_{\text {loc}}, \bar{\mathbb {V}}_{\text {loc}}$ by

$$\begin{aligned} \mathbb {I}_{\text {loc}}(r) =: m \bar{\mathbb {I}}_{\text {loc}}(r) ~~~~ \mathbb {V}_{\text {loc}}(r) =: c \bar{\mathbb {V}}_{\text {loc}}(r). \end{aligned}$$

Defining the dimensionless parameter $\epsilon :=\frac{m}{c}$ and multiplying both sides of (24) by $\mathbb {I}_{\text {loc}}\mathbb {V}_{\text {loc}}^{-1}$, we obtain the rewritten equations of motion

$$\begin{aligned}&\mathring{g} = -A_{\text {mech}}\cdot \dot{r} + \frac{1}{m}\bar{\mathbb {I}}_{\text {loc}}^{-1}p\nonumber \\&\epsilon \bar{\mathbb {I}}_{\text {loc}}\bar{\mathbb {V}}_{\text {loc}}^{-1}\dot{p} = m\bar{\mathbb {I}}_{\text {loc}}(A_{\text {visc}}- A_{\text {mech}})\cdot \dot{r} + p\\&\qquad \quad + \epsilon \bar{\mathbb {I}}_{\text {loc}}\bar{\mathbb {V}}_{\text {loc}}^{-1} \text {ad}^*_{\mathbb {I}_{\text {loc}}^{-1}p}p\nonumber \\&\quad \qquad - \epsilon \bar{\mathbb {I}}_{\text {loc}}\bar{\mathbb {V}}_{\text {loc}}^{-1}\text {ad}^*_{A_{\text {mech}}\cdot \dot{r}}p.\nonumber \end{aligned}$$

(26)

In considering the limit in which viscous forces dominate the inertia of the locomoting body, Kelly and Murray [25] formally set $\epsilon =0$ in (26) to obtain $p = m\bar{\mathbb {I}}_{\text {loc}}(A_{\text {mech}}-A_{\text {visc}})\cdot \dot{r}$ from the second equation. Substituting this into the first equation of (26), they derive the following form of the equations of motion:

$$\begin{aligned} \mathring{g}=-A_{\text {visc}}\cdot \dot{r}. \end{aligned}$$

(27)

In the language of differential geometry, (27) states that in the Stokesian limit trajectories are horizontal with respect to the viscous connection. We will see in the next section that this reduction can be extended away from the $\epsilon \rightarrow 0$ limit.

Appendix B: Reduction in the perturbed Stokes regime

In Eldering and Jacobs [7], the argument of Kelly and Murray [25] was explained in more detail using the theory of normally hyperbolic invariant manifolds (NHIMs) in the context of geometric singular perturbation theory [14, 20, 22]. The idea is to show that for $\epsilon > 0$ sufficiently small, the dynamics (26) possess an exponentially attractive invariant slow manifold$M_\epsilon$, such that the dynamics restricted to $M_\epsilon$ approach (27) as $\epsilon \rightarrow 0$. We give an alternative argument which yields a result differing from that of Eldering and Jacobs [7] in two ways.

• 1. Eldering and Jacobs [7] give an argument for general mechanical systems without symmetry under the assumption that the configuration space Q is compact, although they do indicate that compactness can be replaced with uniformity conditions using noncompact NHIM theory [9]. Our argument assumes symmetry but allows G to be noncompact, though we do require that $S:=Q/G$ be compact. This enables application of our result to locomotion systems with noncompact symmetry groups, such as the Euclidean group of planar rigid motions $\mathsf {SE}(2)$ as in the systems of Sect. 4.

• 2. Eldering and Jacobs [7] consider the limit $m\rightarrow 0$ while holding c and the force exerted by the locomoting body fixed. This makes sense, because if the exerted force were held fixed while taking $c \rightarrow \infty$, then trivial dynamics would result in the singular limit: The system would not move at all. Rather than holding the exerted force fixed, we will consider the differential equation prescribing the dynamics of the shape variable to be fixed.¹⁴ Under this assumption, we show that the dynamics depend only on the ratio$\epsilon = \frac{m}{c}$, and in particular, the dynamics obtained in the two singular limits $m \rightarrow 0$ and $c \rightarrow \infty$ are the same.

Before stating Theorem 2, we need the following definition.

Definition 1

($C^k_b$time-dependent vector fields) Let M be a compact manifold with boundary, and let $f:\mathbb {R}\times M \rightarrow \mathsf {T}M$ a $C^{k \ge 0}$ time-dependent vector field. Let $(U_i)_{i=1}^n$ be a finite open cover of M and $(V_i, \psi _i)_{i=1}^n$ be a finite atlas for M such that $\bar{U}_i \subset V_i$ for all i, and for each i define $f_i :=(\mathsf {D}\psi _i \circ f\circ (\text {id}_\mathbb {R}\times \psi _i^{-1}))$. We define an associated $C^k$ norm ${||}f{||}_k$ of f via

$$\begin{aligned} {||}f{||}_{ k}:=\max _{1\le i\le n}\max _{\begin{array}{c} 0 \le j \le k\\ x\in \psi _i(\bar{U}_i) \end{array}}{||}\mathsf {D}^j f_i(x){||}, \end{aligned}$$

(28)

where ${||}\mathsf {D}^j f_i(x){||}$ denotes the norm of a j-linear map; here $\mathsf {D}^j f$ includes partial derivatives with respect to time as well as the spatial variables. If ${||}f{||}_k< \infty$, we say that f is $C^k$-bounded and write $f\in C^k_b$. The norm ${||}\cdot {||}_k$ makes the $C^k_b$ time-dependent vector fields into a Banach space. The norms induced by any two such finite covers of M are equivalent and thereby induce a canonical $C^k_b$topology on the space of $C^k_b$ time-dependent vector fields.

Remark 4

Definition 1 defines the $C^k_b$ topology on the space of $C^k_b$ time-dependent vector fields on a compact manifold. As discussed in Eldering [9, Sec. 1.7], this $C^k_b$ topology is finer than the $C^k$ weak Whitney topology and coarser than the $C^k$ strong Whitney topology [18, Ch. 2], but all of these topologies induce the same topology on the subspace of time-independent vector fields due to compactness. Definition 1 is a special case of the definition in Eldering [9, Ch. 2] for the $C^k_b$ topology on $C^k_b$ vector fields on Riemannian manifolds of bounded geometry and on $C^k_b$ maps between such manifolds.

The following theorem concerns a G-symmetric dynamical system on $\mathsf {T}Q$ whose equations of motion are consistent with our assumptions so far, i.e., they are given in local trivializations by (26) and an equation of the form (25).

Theorem 2

Assume that $S=Q/G$ is compact. Let $2 \le k < \infty$, and let $X^\epsilon$ be a $C^k$ family of G-symmetric time-dependent vector fields on $\mathsf {T}Q$ with the following properties:

• 1. For every compact neighborhood with $C^k$ boundary $K_0 \subset \mathsf {T}Q$ and $\epsilon > 0$, $X^\epsilon |_{\mathbb {R}\times K_0}\in C^k_b$ (Definition 1).

• 2. There exists a compact connected neighborhood $K\subset \mathsf {T}S$ of the zero section of $\mathsf {T}S$ with $C^k$ boundary, such that $N:=\mathsf {D}\pi _Q^{-1}(K) \subset \mathsf {T}Q$ is positively invariant for $X^\epsilon$, for all sufficiently small $\epsilon > 0$.

• 3. $X^\epsilon$ is given in each local trivialization $\mathsf {T}(U\times G)$, where U is a chart for S, by (25) and (26):

  $$\begin{aligned} \begin{aligned} \ddot{r}&= f\left( t,r,\dot{r},\frac{1}{m}\bar{\mathbb {I}}_{\text {loc}}^{-1}p\right) \\ \epsilon \bar{\mathbb {I}}_{\text {loc}}\bar{\mathbb {V}}_{\text {loc}}^{-1}\dot{p} =&m\bar{\mathbb {I}}_{\text {loc}}(A_{\text {visc}}- A_{\text {mech}})\cdot \dot{r} + p\\&\quad + \epsilon \bar{\mathbb {I}}_{\text {loc}}\bar{\mathbb {V}}_{\text {loc}}^{-1} \text {ad}^*_{\mathbb {I}_{\text {loc}}^{-1}p}p\\&\quad - \epsilon \bar{\mathbb {I}}_{\text {loc}}\bar{\mathbb {V}}_{\text {loc}}^{-1}\text {ad}^*_{A_{\text {mech}}\cdot \dot{r}}p\\ \mathring{g}&= -A_{\text {mech}}\cdot \dot{r} + \frac{1}{m}\bar{\mathbb {I}}_{\text {loc}}^{-1}p \end{aligned} \end{aligned}$$

  (29)

  for some function f which depends on the local trivialization but is independent of $\epsilon$.

Then for all sufficiently small $\epsilon > 0$, there exists a $C^k$ noncompact normally hyperbolic invariant manifold with boundary $M_\epsilon \subset \mathbb {R}\times N \subset \mathbb {R}\times \mathsf {T}Q$ for the extended dynamics given by the extended vector field $(1,X_\epsilon )$ on $\mathbb {R}\times \mathsf {T}Q$. Additionally, $M_\epsilon$ is uniformly (in time and space) globally asymptotically stable and uniformly locally exponentially stable (with respect to the distance induced by any complete G-invariant Riemannian metric on $\mathsf {T}Q$) for the extended dynamics restricted to $\mathbb {R}\times N$. Finally, there exists $\epsilon _0 > 0$ such that, for each local trivialization $U\times G$, there exists a $C^k$ map $h_\epsilon :\mathbb {R}\times (\mathsf {T}U \cap K) \times (0,\epsilon _0)\rightarrow \mathfrak {g}^*$ such that $M_\epsilon \cap \mathsf {D}\pi _Q^{-1}(\mathsf {T}U \cap K)$ corresponds to

$$\begin{aligned} \{(t,r,\dot{r},p,g): p = h_\epsilon (t,r,\dot{r},\epsilon )\}, \end{aligned}$$

(30)

$$\begin{aligned} h_\epsilon (t,r,\dot{r},\epsilon ) = \mathbb {I}_{\text {loc}}\left[ (A_{\text {mech}}(r)-A_{\text {visc}}(r))\cdot \dot{r} + \mathcal {O}(\epsilon )\right] \end{aligned}$$

(with p defined by (22)), and $h_\epsilon$ together with its partial derivatives of order k or less is bounded uniformly in time. If $f(t,r,\dot{r},\mathbb {I}_{\text {loc}}^{-1} p)$ is independent of t, then $h_\epsilon$ and $M_\epsilon$ are independent of t, and $M_\epsilon$ can be interpreted as a compact NHIM for the (non-extended) dynamics restricted to N.

Remark 5

Note that even if we assume $f\in C^\infty$, we can generally only obtain $C^k$ NHIMs $M_\epsilon$ for k finite. This is because we obtain $M_\epsilon$ as a perturbation of a NHIM $M_0$, and perturbations of $C^\infty$ NHIMs are generally only finitely smooth because the maximum perturbation size $\epsilon$ required to obtain degree of smoothness k for $M_\epsilon$ generally depends on k in such a way that $\epsilon \rightarrow 0$ as $k\rightarrow \infty$. See Eldering [9, Rem. 1.12] and Strien [40] for more discussion.

Remark 6

By replacing compactness of Q / G with uniformity conditions, it should be possible to generalize Theorem 2 to the situation of Q noncompact where either Q / G is noncompact, or where there is no symmetry at all. This was pointed out in Eldering and Jacobs [7, App. 1]. This observation seems important for the consideration of dissipative mechanical systems which are only approximately symmetric under a group G, which seems to be a more realistic assumption.

Remark 7

By taking $\epsilon \rightarrow 0$ in Theorem 2, we find that $p = \mathbb {I}_{\text {loc}}(A_{\text {mech}}-A_{\text {visc}})\cdot \dot{r}$ in the limit. Substituting this into the first equation of (32), we obtain Eq. (24) as in Kelly and Murray [25].

Proof

Preparation of the equations of motion. Throughout the proof, we consider the dynamics in local trivializations of the form $U\times G$ for Q, where U is the domain of a chart for S, so that we have induced coordinates $(r,\dot{r})$ for $\mathsf {T}U$. In such a local trivialization, we would like to use (29) to analyze the dynamics, but there are two (related) problems with this. First, the definition of p depends on m, and this will cause difficulties in verifying Definition 1 to check that certain vector fields are close in the $C^k_b$ topology. Second, we would like to analyze (29) in a singular perturbation framework, but this is difficult to do directly because m explicitly appears, and the size of m may or may not be commensurate with the size of $\epsilon$. To remedy this situation, we change variables via the diffeomorphism $(r,\dot{r},p,g)\mapsto (r,\dot{r},\Omega ,g)$ of $\mathsf {T}U \times \mathfrak {g}^* \times G \rightarrow \mathsf {T}U \times \mathfrak {g}\times G$ where $\Omega \in \mathfrak {g}$ is defined by

$$\begin{aligned} \Omega :=\mathbb {I}_{\text {loc}}^{-1}p = \text {Ad}_{g^{-1}}\Gamma _{\text {mech}}(\dot{g},\dot{r})=\mathring{g} + A_{\text {mech}}\cdot \dot{r}.\nonumber \\ \end{aligned}$$

(31)

Sometimes $\Omega$ is referred to as the (body) locked angular velocity [2, p. 61]. Differentiating $\mathbb {I}_{\text {loc}}\Omega = p$, using (29), and rearranging yield

$$\begin{aligned} \begin{aligned}&\dot{t} = 1\\&\dot{r} = v\\&\dot{v} = f(t,r,v,\Omega )\\&\epsilon \dot{\Omega } = -\,\epsilon \bar{\mathbb {I}}_{\text {loc}}^{-1} \left( \frac{\mathrm{d}}{\mathrm{d}t}\bar{\mathbb {I}}_{\text {loc}}\right) \Omega \\&\quad \quad \quad + \bar{\mathbb {I}}_{\text {loc}}^{-1}\bar{\mathbb {V}}_{\text {loc}}(A_{\text {visc}}- A_{\text {mech}})\cdot v + \bar{\mathbb {I}}_{\text {loc}}^{-1}\bar{\mathbb {V}}_{\text {loc}}\Omega \\&\quad \quad \quad + \epsilon \bar{\mathbb {I}}_{\text {loc}}^{-1}\text {ad}^*_{\mathring{g}}\bar{\mathbb {I}}_{\text {loc}}\Omega , \end{aligned} \end{aligned}$$

(32)

where we have introduced the variable $v:=\dot{r}$. We have written $\text {ad}^*_{\mathring{g}}$ for space reasons, but note that the $\dot{\Omega }$ equation is independent of g since

$$\begin{aligned} \mathring{g} = -A_{\text {mech}}\cdot \dot{r} + \Omega , \end{aligned}$$

(33)

and this implies that $\text {ad}^*_{\mathring{g}} = \text {ad}^*_{\Omega }-\text {ad}^*_{A_{\text {mech}}\cdot \dot{r}}$. We see that (32) is split into slow (t, r, v) and fast $(\Omega )$ variables, which is the appropriate setup for a singular perturbation analysis. The remainder of the proof consists of two parts: (i) proving that the NHIM $M_\epsilon$ exists and (ii) establishing the stability properties of $M_\epsilon$.

Proof that$M_\epsilon$exists. Introducing the “fast time” $\tau :=\frac{1}{\epsilon } t$ and denoting a derivative with respect to $\tau$ by a prime, after the time rescaling we obtain the regularized equations

$$\begin{aligned} t'&= \epsilon \nonumber \\ r'&= \epsilon v\nonumber \\ v'&= \epsilon f(t,r,v,\Omega )\nonumber \\ \Omega '&= -\epsilon \bar{\mathbb {I}}_{\text {loc}}^{-1} \left( \frac{\mathrm{d}}{\mathrm{d}t}\bar{\mathbb {I}}_{\text {loc}}\right) \Omega \nonumber \\&\quad + \bar{\mathbb {I}}_{\text {loc}}^{-1}\bar{\mathbb {V}}_{\text {loc}}(A_{\text {visc}}- A_{\text {mech}})\cdot v\nonumber \\&\quad + \bar{\mathbb {I}}_{\text {loc}}^{-1}\bar{\mathbb {V}}_{\text {loc}}\Omega + \epsilon \bar{\mathbb {I}}_{\text {loc}}^{-1}\text {ad}^*_{\mathring{g}}\bar{\mathbb {I}}_{\text {loc}}\Omega . \end{aligned}$$

(34)

This rescaling of time is equivalent to replacing the vector field $(1,X_\epsilon )$ on $\mathbb {R}\times \mathsf {T}Q$ by $(\epsilon ,\epsilon X_\epsilon )$. We see from (33) and (34) that there is a well-defined $C^k$ time-dependent vector field $\tilde{X}_0$ given by the pointwise limit $\tilde{X}_0:=\lim _{\epsilon \rightarrow 0} \epsilon X_\epsilon$. Given any G-symmetric time-dependent vector field Y on $\mathsf {T}Q$, we let Y / G denote the corresponding reduced vector field on $(\mathsf {T}Q)/G$. Hence, (34) shows that the extended vector field $(1,\tilde{X}_0/G)$ has a smooth embedded submanifold $(M_0/G)$ of critical points whose intersection with a locally trivializable neighborhood is given by

$$\begin{aligned} \{(r,v,\Omega ) \in \mathsf {T}U \times \mathfrak {g}: \Omega = (A_{\text {mech}}-A_{\text {visc}})\cdot v\}, \end{aligned}$$

(35)

and it is readily seen that $M_0/G$ is described globally as the quotient of the Ehresmann connection $M_0 :=\ker \Gamma _{\text {visc}}$ by the lifted action of G on $\mathsf {T}Q$.

Furthermore, $M_0/G$ is a globally exponentially stable NHIM for the $\epsilon = 0$ system. To see this, first note that in any local trivialization t, r, v are constants when $\epsilon = 0$, and hence, $\Omega '$ is of the form $\Omega ' = \bar{\mathbb {I}}_{\text {loc}}^{-1} \bar{\mathbb {V}}_{\text {loc}}\Omega + b$ for a constant b and therefore has a globally exponentially stable equilibrium provided that all eigenvalues of $\bar{\mathbb {I}}_{\text {loc}}^{-1}\bar{\mathbb {V}}_{\text {loc}}$ have negative real part. To see that this is the case, fix a basis of $\mathfrak {g}$ and corresponding dual basis for $\mathfrak {g}^*$, and first consider the product $\mathbb {I}^{-1} \mathbb {V}$. With respect to our chosen basis, $\mathbb {I}, \mathbb {V}$ and their inverses $\mathbb {I}^{-1}, \mathbb {V}^{-1}$ are, respectively, represented by r-dependent matrices $I_{ij}, V_{ij}$ and their inverses $I^{ij}, V^{ij}$. It is immediate from the definitions (15) and (17) that $I_{ij}$ and $V_{ij}$ are, respectively, positive definite and negative definite symmetric matrices (this is why we required the bilinear forms $k, \nu$ to be positive definite when restricted to vectors tangent to G orbits). Since $I_{ij}$ is symmetric positive definite, we may let $(\sqrt{I})_{ij}$ be a matrix square root of $I_{ij}$ and let $(\sqrt{I})^{ij}$ be its inverse. But then the product $I^{ik}V_{kj}$ is similar to the symmetric negative definite matrix $(\sqrt{I})^{ik}V_{k\ell }(\sqrt{I})^{\ell j}$ (Einstein summation implied). Hence, $\mathbb {I}^{-1} \mathbb {V}$ has only eigenvalues with negative real part, and the same is true of $\mathbb {I}_{\text {loc}}^{-1}\mathbb {V}_{\text {loc}}$ because of the similarity $\mathbb {I}_{\text {loc}}^{-1}\mathbb {V}_{\text {loc}}= \text {Ad}_g^{-1} \mathbb {I}^{-1}\mathbb {V}\text {Ad}_g$.

Let $\tilde{\pi }:(\mathsf {T}Q)/G \rightarrow \mathsf {T}S$ denote the projection induced by $\mathsf {D}\pi _Q$. Equation (35) implies that $M_0/G$ is the image of a section $\sigma _0:\mathsf {T}S \rightarrow (\mathsf {T}Q)/G$ of $\tilde{\pi }$. Hence, $(M_0/G)\cap \tilde{\pi }^{-1}(K) = \sigma _0(K)$ is compact, and $M_0/G$ intersects $\tilde{\pi }^{-1}(\partial K)$ transversely. Furthermore, the assumption that $X^\epsilon |_{\mathbb {R}\times K_0} \in C^k_b$ for any compact neighborhood with $C^k$ boundary $K_0 \subset \mathsf {T}Q$ implies that all partial derivatives of f are bounded on compact sets uniformly in time. This makes it clear that for any compact $K_1\subset (\mathsf {T}Q)/G$, $(\epsilon X_\epsilon /G)|_{\mathbb {R}\times K_1}$ can be made arbitrarily close to $(\tilde{X}_0/G)|_{\mathbb {R}\times K_1}$ in the $C^k_b$ topology (Definition 1) by taking $\epsilon > 0$ sufficiently small. Hence, by the noncompact NHIM results of Eldering [9, Sec. 4.1-4.2], it follows that $(M_0/G)\cap \tilde{\pi }^{-1}(K)$ persists in extended state space $\mathbb {R}\times N$ to a nearby attracting NHIM $M_\epsilon /G$ with boundary for $(\epsilon , \epsilon X_\epsilon /G)$.¹⁵ Furthermore, $M_\epsilon /G$ is the image of a section $\sigma _\epsilon :\mathbb {R}\times K \rightarrow (\mathsf {T}Q)/G$ of $\tilde{\pi }$ and is given in each local trivialization of $(\mathsf {T}Q)/G$ by the graph of a function $\Omega = \tilde{h}_\epsilon (t,r,\dot{r},\epsilon )$ which is $C^k$ bounded uniformly in time. By symmetry, the preimage $M_\epsilon = \pi _{\mathsf {T}Q}^{-1}(M_\epsilon /G)$ of $M_\epsilon /G$ via the quotient $\pi _{\mathsf {T}Q}:\mathsf {T}Q \rightarrow (\mathsf {T}Q)/G$ yields a NHIM $M_\epsilon$ for $(\epsilon , \epsilon X_\epsilon )$ (and hence also for $(1,X_\epsilon )$) on the subset $\mathbb {R}\times N$ of $\mathbb {R}\times \mathsf {T}Q$, and $M_\epsilon$ is given in each local trivialization by the graph of the same function $\Omega = \tilde{h}_\epsilon$ as $M_\epsilon /G$ but augmented with trivial dependence on g. The function $h_\epsilon$ from the theorem statement is given by $h_\epsilon = \mathbb {I}_{\text {loc}}\tilde{h}_\epsilon$.

Proof of the stability properties of$M_\epsilon$. Fix any complete G-invariant Riemannian metric on¹⁶$\mathsf {T}Q$, so that it descends to a metric on $(\mathsf {T}Q)/G$ making $\pi _{\mathsf {T}Q}:\mathsf {T}Q \rightarrow (\mathsf {T}Q)/G$ into a Riemannian submersion [6, p. 185]. We have distance functions $\tilde{d}$ and d on $\mathsf {T}Q$ and $(\mathsf {T}Q)/G$ induced by these metrics. For $t \in \mathbb {R}$, we let $M_\epsilon (t):=M_\epsilon \cap (\{t\}\times N)$ and $M_\epsilon (t)/G:=\pi _{\mathsf {T}Q}(M_\epsilon (t))$. Given $w\in \mathsf {T}Q$ and its orbit $\pi _{\mathsf {T}Q}(w) \in (\mathsf {T}Q)/G$, it follows that for all $t\in \mathbb {R}$, $\tilde{d}(w, M_\epsilon (t)) = d(\pi _{\mathsf {T}Q}(w), M_\epsilon (t)/G)$.¹⁷ Hence, it suffices to prove that $M_\epsilon /G$ is uniformly globally asymptotically stable and locally exponentially stable for the vector field $(1, X_\epsilon /G)$ on $\mathbb {R}\times \tilde{\pi }^{-1}(K) = \mathbb {R}\times \pi _{\mathsf {T}Q}(N)$, and to do this it suffices to prove the same for $(\epsilon , \epsilon X_\epsilon /G)$.

Fixing an inner product $\langle \,\cdot \,, \,\cdot \,\rangle$ and associated norm ${||}\,\cdot \,{||}$ on $\mathfrak {g}$, we accomplish this in two steps. First, we show that there exists a compact neighborhood $K_0 \subset \pi _{\mathsf {T}Q}(N)$ of $M_\epsilon /G$ such that $K_0$ is positively invariant for the time-dependent flow of $X_\epsilon$ and such that any other compact neighborhood $K_1\subset \pi _{\mathsf {T}Q}(N)$ of $M_\epsilon /G$ flows into $K_0$ after some finite time depending on $K_1$ but independent of the initial time. Second, we show that all trajectories in $K_0$ converge to $M_\epsilon /G$ at a uniform exponential rate. To achieve this second step, we show that in the intersection of each local trivialization with $K_0$, ${||}\Omega -\tilde{h}_\epsilon (t,r,v){||}$ decreases at an exponential rate. Since $(\mathsf {T}Q)/G$ is covered by finitely many local trivializations (by compactness of S) and since all Riemannian metrics are uniformly equivalent on compact sets,¹⁸ this will establish uniform exponential convergence of points in $K_0$ with respect to the distance induced by any Riemannian metric, and in particular the distance d.

Consider a local trivialization $U\times G$ of Q and the associated form (34) of the dynamics restricted to $\tilde{\pi }^{-1}(K \cap \mathsf {T}U)$. Differentiating ${||}\Omega {||}^2$ using the last equation of (34), it is easy to check that $\frac{d}{d\tau }{||}\Omega {||}^2 \rightarrow -\infty$ as ${||}\Omega {||}^2 \rightarrow \infty$, uniformly in $(t,r,v,\epsilon )$ for $\epsilon$ sufficiently small. (This follows from the negative definiteness of $\mathbb {I}_{\text {loc}}^{-1}\mathbb {V}_{\text {loc}}$ and the compactness of K.) Hence, we see that there exists $k_0 > 0$ such that for all $\epsilon$ sufficiently small, $\frac{d}{d\tau }{||}\Omega {||}^2 \le -1$ when ${||}\Omega {||}^2 \ge k_0^2$. Now $k_0$ might depend on the local trivialization, but we can replace $k_0$ with the largest such constant selected from finitely many fixed local trivializations covering Q. Hence, there exists a compact subset $K_0 \subset \pi _{\mathsf {T}Q}(N)$ given by $\{{||}\Omega {||} \le k_0\}$ in each of these fixed local trivializations, such that $K_0$ is positively invariant for the time-dependent flow of $X_\epsilon$ and such that any other compact neighborhood $K_1 \subset \pi _{\mathsf {T}Q}(N)$ of $M_\epsilon /G$ flows into $K_0$ after some finite time independent of the initial time.

It remains only to establish the uniform exponential rate of convergence of trajectories in $K_0$ to $M_\epsilon$. For each local trivialization $U \times G$ of Q, we define the translated variable $\tilde{\Omega }:=\Omega - \tilde{h}_\epsilon (t,r,v,\epsilon )$. Since $M_\epsilon /G$ is invariant, we must have $\tilde{\Omega }' = 0$ whenever $\tilde{\Omega } = 0$. Differentiating $\tilde{\Omega }$ using (34), we therefore find that

$$\begin{aligned} \tilde{\Omega }'&= \left[ -\epsilon \bar{\mathbb {I}}_{\text {loc}}^{-1} \left( \frac{\mathrm{d}}{\mathrm{d}t}\bar{\mathbb {I}}_{\text {loc}}\right) \right. \nonumber \\&\quad \left. +\, \epsilon \bar{\mathbb {I}}_{\text {loc}}^{-1}\text {ad}^*_{\mathring{g}}\bar{\mathbb {I}}_{\text {loc}}+ \epsilon \zeta (t,r,v,\tilde{\Omega }) + \bar{\mathbb {I}}_{\text {loc}}^{-1}\bar{\mathbb {V}}_{\text {loc}}\right] \tilde{\Omega }\nonumber \\&=: \left[ \epsilon A(t,r,v,\tilde{\Omega }) + \bar{\mathbb {I}}_{\text {loc}}^{-1}\bar{\mathbb {V}}_{\text {loc}}(r) \right] \tilde{\Omega }, \end{aligned}$$

(36)

since all of the terms which do not vanish when $\tilde{\Omega } = 0$ must cancel. Here $\zeta$ is defined via Hadamard’s lemma [32, Lemma 2.8]:

$$\begin{aligned}&\zeta (t,r,v,\tilde{\Omega }) \nonumber \\&\quad :=\frac{\partial }{\partial v}\tilde{h}_\epsilon (t,r,v) \int _{0}^{1}\frac{\partial }{\partial \Omega }f(t,r,v,\tilde{h}_\epsilon (t,r,v) \nonumber \\&\quad \quad +\, s \tilde{\Omega })\,ds, \end{aligned}$$

(37)

so that $\zeta (t,r,v,\tilde{\Omega })\tilde{\Omega } = \tilde{h}_\epsilon (t,r,v) f(t,r,v,\tilde{h}_\epsilon + \tilde{\Omega })$. As previously mentioned, the $C^k$ boundedness of $X_\epsilon$ on compact subsets of $\mathsf {T}Q$ implies that $\tilde{h}_\epsilon$, f, and their first k partial derivatives are uniformly bounded on sets of the form $\mathbb {R}\times K_2$ with $K_2$ compact. Hence, whenever ${||}\Omega {||}\le k_0$ and $(r,v) \in U \cap K$, ${||}A(t,r,v,\tilde{\Omega }){||} \le L$ for some constant L depending on the local trivialization; we replace L with the largest such constant chosen from finitely many local trivializations covering Q. Integrating both sides of (36), taking norms using the triangle inequality, and applying Grönwall’s Lemma therefore yield

$$\begin{aligned} \begin{aligned} {||}\tilde{\Omega }(\tau ){||}&\le e^{-\lambda (\tau -\tau _0)} e^{\int _{\tau _0}^{\tau }\epsilon {||}A(t(s),r(s),v(s),\tilde{\Omega }(s){||}\,ds }\\&{||}\tilde{\Omega }(\tau _0){||}\\&\le e^{\left[ -\lambda + \epsilon L \right] (\tau -\tau _0)} {||}\tilde{\Omega }(\tau _0){||}. \end{aligned} \end{aligned}$$

(38)

where $-\lambda < 0$ is defined via $-\lambda :=\sup _{r \in S} \max \, \text {spec}(\bar{\mathbb {I}}_{\text {loc}}^{-1}\bar{\mathbb {V}}_{\text {loc}}(r))$ and is strictly negative since S is compact. By the previous discussion, requiring $\epsilon > 0$ to be sufficiently small so that $-\lambda + \epsilon L < 0$ completes the proof. $\square$

Theorem 2 and Remark 7 show that, to zeroth order in $\epsilon$, the dynamics restricted to the slow manifold $M_\epsilon$ are given by the viscous connection model (27). The following theorem shows that the dynamics restricted to $M_\epsilon$ can be explicitly computed to higher order in $\epsilon$. We compute the restricted dynamics to first order in $\epsilon$. Higher-order terms in $\epsilon$ can also be computed recursively, but we choose not to pursue this here.

Theorem 3

Assume the same hypotheses as in Theorem 2. Then the dynamics restricted to the slow manifold $M_\epsilon$ are given in a local trivialization by

$$\begin{aligned} \mathring{g}= & {} -A_{\text {visc}}\cdot \dot{r} + \epsilon \bar{\mathbb {V}}_{\text {loc}}^{-1} \left( \left( \frac{\partial }{\partial _r} \bar{h}_0\right) \dot{r}\right. \nonumber \\&\left. + \left( \frac{\partial }{\partial \dot{r}}\bar{h}_0\right) \ddot{r} - \text {ad}^*_{\mathring{g}}(\bar{h}_0)\right) + \mathcal {O}(\epsilon ^2), \end{aligned}$$

(39)

where

$$\begin{aligned} \bar{h}_0(r,\dot{r}) :=\frac{1}{m}h_0(r,\dot{r}) = \bar{\mathbb {I}}_{\text {loc}}(A_{\text {mech}}(r)-A_{\text {visc}}(r))\cdot \dot{r}, \end{aligned}$$

where we are using the definition $\bar{\mathbb {I}}_{\text {loc}}:=\frac{1}{m}\mathbb {I}_{\text {loc}}$. Alternatively, we may write

$$\begin{aligned} \mathring{g}= & {} -A_{\text {visc}}\cdot \dot{r} + \epsilon \bar{\mathbb {V}}_{\text {loc}}^{-1} \left( \left( \frac{\partial }{\partial r} \bar{h}_0\right) \dot{r}\right. \nonumber \\&\left. +\, \left( \frac{\partial }{\partial \dot{r}}\bar{h}_0\right) f(t,r,\dot{r},\bar{\mathbb {I}}_{\text {loc}}^{-1}\bar{h}_0) - \text {ad}^*_{\mathring{g}}(\bar{h}_0)\right) \nonumber \\&+\, \mathcal {O}(\epsilon ^2), \end{aligned}$$

(40)

for a different $\mathcal {O}(\epsilon ^2)$ term.

Remark 8

Notice the presence, in the second term of (39), of $\bar{h}_0$ rather than $h_0$ of (30). This is important because the expression for $h_0$ contains an $\mathbb {I}_{\text {loc}}= m \bar{\mathbb {I}}_{\text {loc}}$ factor. Because of the possibility that the size of m is commensurate with $\epsilon$, this means that $h_0$ could be $\mathcal {O}(\epsilon )$. However, $\bar{h}_0$ is $\mathcal {O}(1)$, ensuring that the second term is $\mathcal {O}(\epsilon )$ but not $\mathcal {O}(\epsilon ^2)$.

Remark 9

Equations (39) and (40) can be viewed as adding $\mathcal {O}(\epsilon )$ correction terms to the viscous connection model (27), valid in the limit $\epsilon \rightarrow 0$, to account for the more realistic situation that the inertia–damping ratio $\frac{m}{c} = \epsilon$ is small but nonzero.

Proof of Theorem 3

Consider the function

$$\begin{aligned} \tilde{h}_\epsilon (t,r,\dot{r},\epsilon ) :=\mathbb {I}_{\text {loc}}^{-1} h_\epsilon {=} (A_{\text {mech}}(r){-}A_{\text {visc}}(r))\cdot \dot{r} {+} \mathcal {O}(\epsilon ) \end{aligned}$$

from the proof of Theorem 2, and define $\bar{h}_\epsilon :=\bar{\mathbb {I}}_{\text {loc}}\tilde{h}_\epsilon = \frac{1}{m}h_\epsilon$. Since $\bar{h}_\epsilon , \tilde{h}_\epsilon \in C^k$, we may expand them as asymptotic series

$$\begin{aligned} \begin{aligned} \bar{h}_\epsilon&= \bar{h}_0 + \epsilon \bar{h}_1 + \cdots + \epsilon ^{k} \bar{h}_{k} + \mathcal {O}(\epsilon ^{k+1})\\ \tilde{h}_\epsilon&= \tilde{h}_0 + \epsilon \tilde{h}_1 + \cdots + \epsilon ^{k} \tilde{h}_{k} + \mathcal {O}(\epsilon ^{k+1}), \end{aligned} \end{aligned}$$

(41)

where for all i, $\bar{h}_i = \bar{\mathbb {I}}_{\text {loc}}\tilde{h}_i$. We also already know from Theorem 2 that $\tilde{h}_0 =(A_{\text {mech}}-A_{\text {visc}})\cdot \dot{r}$, and therefore, $\tilde{h}_0(t,r,\dot{r}) \equiv \tilde{h}_0(r,\dot{r})$ has no explicit t-dependence. We now compute $\tilde{h}_1$ via a standard technique [20]. Differentiating both sides of the equation $\Omega = \tilde{h}_\epsilon (t,r,\dot{r},\epsilon )$ with respect to time (using (32) to differentiate the left hand side), substituting the second equation of (41) for $\Omega$ in the resulting expression, and retaining terms only up to $\mathcal {O}(\epsilon )$, we obtain

$$\begin{aligned}&-\epsilon \bar{\mathbb {I}}_{\text {loc}}^{-1} \left( \frac{\mathrm{d}}{\mathrm{d}t}\bar{\mathbb {I}}_{\text {loc}}\right) \tilde{h}_0 + \bar{\mathbb {I}}_{\text {loc}}^{-1}\bar{\mathbb {V}}_{\text {loc}}(A_{\text {visc}}- A_{\text {mech}})\cdot \dot{r}\\&\quad + \bar{\mathbb {I}}_{\text {loc}}^{-1}\bar{\mathbb {V}}_{\text {loc}}\left( \tilde{h}_0+\epsilon \tilde{h}_1 \right) \\&\quad + \epsilon \bar{\mathbb {I}}_{\text {loc}}^{-1}\text {ad}^*_{\mathring{g}}\bar{\mathbb {I}}_{\text {loc}}\tilde{h}_0 = \epsilon \dot{\tilde{h}}_0 + \mathcal {O}(\epsilon ^2). \end{aligned}$$

Equating the coefficients of $\epsilon$ yields

$$\begin{aligned} \tilde{h}_1&= \bar{\mathbb {V}}_{\text {loc}}^{-1} \left( \frac{\mathrm{d}}{\mathrm{d}t}\bar{\mathbb {I}}_{\text {loc}}\right) \tilde{h}_0+ \bar{\mathbb {V}}_{\text {loc}}^{-1}\bar{\mathbb {I}}_{\text {loc}}\dot{\tilde{h}}_0 - \bar{\mathbb {V}}_{\text {loc}}^{-1}\text {ad}^*_{\mathring{g}}\bar{\mathbb {I}}_{\text {loc}}\tilde{h}_0\\&= \bar{\mathbb {V}}_{\text {loc}}^{-1}\frac{\mathrm{d}}{\mathrm{d}t}\left( \bar{\mathbb {I}}_{\text {loc}}\tilde{h}_0 \right) - \bar{\mathbb {V}}_{\text {loc}}^{-1}\text {ad}^*_{\mathring{g}}\bar{\mathbb {I}}_{\text {loc}}\tilde{h}_0. \end{aligned}$$

Since $h_1 = \mathbb {I}_{\text {loc}}\tilde{h}_1$ and $\bar{h}_0 = \bar{\mathbb {I}}_{\text {loc}}\tilde{h}_0$, we find

$$\begin{aligned} h_1 = \mathbb {I}_{\text {loc}}\bar{\mathbb {V}}_{\text {loc}}^{-1}\frac{\mathrm{d}}{\mathrm{d}t}\left( \bar{h}_0 \right) - \mathbb {I}_{\text {loc}}\bar{\mathbb {V}}_{\text {loc}}^{-1}\text {ad}^*_{\mathring{g}}\left( \bar{h}_0\right) , \end{aligned}$$

(42)

and therefore (substituting $\ddot{r} = f(t,r,\dot{r},\mathbb {I}_{\text {loc}}^{-1} p) = f(t,r,\dot{r},\tilde{h}_0) + \mathcal {O}(\epsilon )$ and differentiating $\bar{h}_0(r,\dot{r})$ via the chain rule),

$$\begin{aligned} h_\epsilon (t,r,\dot{r},\epsilon )&= \mathbb {I}_{\text {loc}}(A_{\text {mech}}-A_{\text {visc}})\cdot \dot{r} \nonumber \\&\quad + \epsilon \mathbb {I}_{\text {loc}}\bar{\mathbb {V}}_{\text {loc}}^{-1} \left( \left( \frac{\partial }{\partial r} \bar{h}_0\right) \dot{r}\right. \nonumber \\&\quad \left. + \left( \frac{\partial }{\partial \dot{r}}\bar{h}_0\right) f(t,r,\dot{r},\tilde{h}_0) - \text {ad}^*_{\mathring{g}}(\bar{h}_0)\right) \nonumber \\&\quad + \mathbb {I}_{\text {loc}}\mathcal {O}(\epsilon ^2). \end{aligned}$$

(43)

Notice that, since $\tilde{h}_0$ is a function of $r,\dot{r}$ only, the $\mathcal {O}(\epsilon )$ portion of the right-hand side of (43) is a function of $t,r,\dot{r}$ alone and not p. This is required since $h_\epsilon$ is required to be a function of $t,r,\dot{r},\epsilon$ alone, and is the reason that we needed to replace $\ddot{r}$ by $f(t,r,\dot{r},\tilde{h}_0)$ in the $\mathcal {O}(\epsilon )$ term. Substituting (43) into the first equation of (26) yields Eq. (40). Finally, making the substitution $f(t,r,\dot{r},\tilde{h}_0) = \ddot{r} + \mathcal {O}(\epsilon )$ in Eq. (40) yields Eq. (39). $\square$

The following theorem makes clearer the functional form of the dynamics (39), and it removes the $\mathring{g}$ dependence of the right-hand side of (39).

Theorem 1$'$Assume the hypotheses of Theorem 2. For sufficiently small$\epsilon > 0,$then for each local trivialization there exist smooth fields of linear maps B(r) and (1, 2) tensors G(r) such that the dynamics restricted to the slow manifold$M_\epsilon$in the local trivialization satisfy

$$\begin{aligned} \mathring{g} {=} -A_{\text {visc}}(r) \cdot \dot{r} {+} \epsilon B(r)\cdot \ddot{r} +\epsilon G(r)\cdot (\dot{r},\dot{r}) + \mathcal {O}(\epsilon ^2).\nonumber \\ \end{aligned}$$

(44)

Remark 10

The (1, 2) tensors G(r) are not generally symmetric, which is clear from Eq. (46).

Proof

Using the properties of $\text {ad}^*$, we may write $\text {ad}^*_{\mathring{g}}(\bar{h}_0) = (C\cdot \bar{h}_0)\cdot (\mathring{g})$ for an appropriate (r-independent) linear map $C:\mathfrak {g}^* \rightarrow \text {End}(\mathfrak {g})$, and hence, we may rewrite (39) as

$$\begin{aligned}&(\text {id}_\mathfrak {g}+ \epsilon \bar{\mathbb {V}}_{\text {loc}}^{-1} (C\cdot \bar{h}_0) )\cdot (\mathring{g})\\&\quad = -A_{\text {visc}}\cdot \dot{r} +\epsilon \bar{\mathbb {V}}_{\text {loc}}^{-1} \left( \left( \frac{\partial }{\partial r} \bar{h}_0\right) \dot{r} + \left( \frac{\partial }{\partial \dot{r}}\bar{h}_0\right) \ddot{r}\right) \\&\qquad +\, \mathcal {O}(\epsilon ^2). \end{aligned}$$

For sufficiently small $\epsilon$, we may use the identity

$$\begin{aligned}&(\text {id}_\mathfrak {g}+ \epsilon \bar{\mathbb {V}}_{\text {loc}}^{-1} (C\cdot \bar{h}_0) )^{-1}\\&\quad = \text {id}_\mathfrak {g}- \epsilon \bar{\mathbb {V}}_{\text {loc}}^{-1} (C\cdot \bar{h}_0) + \mathcal {O}(\epsilon ^2) \end{aligned}$$

to obtain

$$\begin{aligned} \mathring{g}= & {} -A_{\text {visc}}\cdot \dot{r} + \epsilon \bar{\mathbb {V}}_{\text {loc}}^{-1}(C\cdot \bar{h}_0)\cdot A_{\text {visc}}\cdot \dot{r}\nonumber \\&+\, \epsilon \bar{\mathbb {V}}_{\text {loc}}^{-1} \left( \frac{\partial }{\partial r} \bar{h}_0\right) \dot{r} \nonumber \\&+\, \epsilon \bar{\mathbb {V}}_{\text {loc}}^{-1} \left( \frac{\partial }{\partial \dot{r}}\bar{h}_0\right) \ddot{r} + \mathcal {O}(\epsilon ^2). \end{aligned}$$

(45)

Since $\bar{h}_0(r,\dot{r})= \bar{\mathbb {I}}_{\text {loc}}(r)(A_{\text {mech}}(r)-A_{\text {visc}}(r))\cdot \dot{r}$ is linear in $\dot{r}$, it follows that the second and third terms are bilinear in $\dot{r}$, and the fourth term is linear in $\ddot{r}$. Hence, we may take $B(r):=\bar{\mathbb {V}}_{\text {loc}}^{-1} \left( \frac{\partial }{\partial \dot{r}}\bar{h}_0\right)$ and

$$\begin{aligned} G(r)\cdot (\dot{r},\dot{r}):= & {} \bar{\mathbb {V}}_{\text {loc}}^{-1}(C\cdot \mathbb {I}_{\text {loc}}(A_{\text {mech}}-A_{\text {visc}})\cdot \dot{r})\nonumber \\&\cdot \, A_{\text {visc}}\cdot \dot{r} + \epsilon \bar{\mathbb {V}}_{\text {loc}}^{-1} \frac{\partial }{\partial r}\nonumber \\&\left( \mathbb {I}_{\text {loc}}(A_{\text {mech}}-A_{\text {visc}})\cdot \dot{r}\right) \cdot \dot{r}. \end{aligned}$$

(46)

$\square$