{'fraction_non_alphanumeric': 0.04152159525796334, 'fraction_numerical': 0.0066181336863004635, 'mean_word_length': 4.301906941266209, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Differential Geometric Formulations of Slow Invariant Manifold Computation: Geodesic Stretching and Flow Curvature

The theory of slow invariant manifolds (SIMs) is the foundation of various model-order reduction techniques for dissipative dynamical systems with multiple time-scales, e.g. in chemical kinetic models. The construction of SIMs and many approximation methods exploit the restrictive requirement of an explicit time-scale separation parameter. Most of those methods are also not formulated covariantly, i.e. in terms of tensorial constructions. We propose an intrinsically coordinate-free differential geometric approximation criterion approximating normally attracting invariant manifolds (NAIMs). We translate some ideas behind existing approximation approaches, the stretching based diagnostics (SBD) and the flow curvature method (FCM) to tensors of Riemannian geometry, specifically to spacetime curvature in extended phase space. For that purpose we derive from flow-generating smooth vector fields a metric tensor such that the original dynamical system is a geodesic flow on a Riemannian manifold. We apply the resulting method to test models.


Introduction
A wide range of natural processes are modeled by high dimensional dynamical systems with multiple time-scales, for example in chemical kinetics. Their numerical treatment is challenging due to high dimension and stiffness resulting from spectral gaps. The existence of different time-scales usually correlates with a bundling behavior of solution trajectories near a lower-dimensional invariant (FCM) [14,15] can used as guidance, when a suitable framework is used. The SBD refers to normal hyperbolicity. Translating the SBD approach into such a tensorial setting is a major focus of this work.
The paper is organized as follows: The foundations of normal hyperbolicity and GSPT as well as the benefit of tensorial formulations are discussed in Section 1.1, Section 1.2 and Section 1.3 respectively. A suitable geometric setting is motivated and developed in Section 2. The resulting setting provides the possibility to make use of various notions of intrinsic curvature to analyze the bundling behavior of trajectories near the NAIMs. In Section 3, we derive which of these notions are an evident choice by exploiting the geometrical foundation of the SBD [1]. In Section 4 we introduce new viewpoints on the FCM and illustrate, how it can reformulated in our Riemannian geometry setting.

Normal Hyperbolicity
The concepts of normal hyperbolicity is central for the tensorial SIM approximation proposed in Section 2 and Section 3. NHIMs have the defining property to admit a hyperbolic splitting into whitney sum of flow-invariant subbundles of the tangent bundle: Here, stable and unstable subbundles are indicated by E s and E u respectively. Since we are primarily interested in exponentially attracting manifolds. Hence, we restrict ourself to the case without unstable directions (though the following approach can be modified to suit the more general case). The sum (1)  The latter property is formalized by considering so-called generalized Lyapunovtype numbers (see [41]), which is discussed in Section 3.2 in more detail.

Geometric Singular Perturbation Theory
The established analytical foundation of the theory of SIMs is introduced in [9][10][11][12] and can be applied to explicit slow-fast systems. This is a class of dynamical systems which can be written in the form where 0 < ε 1. In the former setting, x s and x f are called slow and fast variables respectively, such that x = [x s , x f ] and n s + n f = n. Fenichel showed (see [11]) a persistence property under certain conditions: If there exists a NHIM for the vector field associated with ε = 0, then there is also a NHIM for small ε > 0. If a NHIM is constructed in this specific context in this particular way, we call it a SIM. In case of exponentially attracting invariant manifolds, a SIM is also a NAIM. In this case, a SIM is represented by the mapping where h ε can be expressed by a power series (asymptotic expansion) in ε: The functions h k (x f ) in eq. (2) are iteratively calculated by use of the so-called invariance equation and matching of the coefficient with respect to ε-powers (matched asymptotic expansion).
Many of the computational methods referred to in the Introduction aim at an approximation of the ε-Taylor series to some order in some coordinate system. Our aim is to find a local tensorial, coordinate independent, intrinsically geometric quantity approximating normal hyperbolicity itself, without restriction to singular perturbed slow-fast systems. We also obtain a SIM approximation for slow-fast systems as a mere consequence.

Advantages of a covariant formulation
We use tensors fields as the fundamental concept in this work to propose a suitable covariant approximation method. Tensors are multilinear mappings from a cartesian product of vector spaces and dual spaces to R. In our case, the vector space is the tangent space T p M of the solution manifold M of the dynamical systemẋ = f (x) in extended phase space (including a time coordinate). We base our approach on purely geometric concepts and choose our tensors accordingly: (a) In Section 2.2 we introduce a specific metric (rank two tensor).
(b) The metric is utilized to define the Riemann curvature tensor (rank four tensor) in Section 3.1.
(c) By plugging-in tangent vectors (rank one tensors) we then receive the deviation tensor (rank two) and a scalar value corresponding to a certain curvature (rank zero tensor) in Section 3.3.
For a given dynamical systemẋ = f (x), the set of variables X := {x 1 , . . . , x n } together with a time-variable τ forms a canonical coordinate-system for the space-time manifold M . All calculations in the above steps (a) -(c) are based on the bundling behavior of the solution trajectories ofẋ = f (x) regarding the so-called 'parent coordinates' X and a vector field f expressed in these coordinates and thus defining the solution manifold geometrically. The parent coordinates induce a basis B X of each tangent space (and cotangent space). Each tensor T can then be expressed by its coefficients in the X-coordinate frame. Conversely, once a basis and all coefficients regarding that basis are specified, T is well-defined as a coordinate-free object in the following manner: Its action on every applicable combination of vectors and covectors does not depend on the choice of the basis. We can take another basis B Y (induced by a chart Y ) and the coefficients of T transform in a certain manner, as do the coefficients of vectors and covectors we are inserting. This is what covariant means in our context: The definite way tensor coefficients change when transforming the basis or equivalently the coordinate system.
The scalar curvature value in (see (c)) is used to determine the approximate NAIM location. Crucially, this curvature value is just the evaluation of a tensor and does not depend on the used basis. Hence, we can choose our own (socalled 'utilized') coordinate chart Y by applying some transformation y = φ(x) and calculate steps (a) -(c) regarding Y . We receive the same tensors -still containing the same information of the original systemẋ = f (x) -but now with different coefficients regarding Y .
This approach is applicable for any system of the formẋ = f (x) (as long as f is smooth enough) and without having a presumed division into slow and fast variables. We can choose slow and fast variables or transform the coordinates entirely, according to our needs.
Remark: This approach requires us to originally choose 'parent coordinates' X in which the dynamic systemẋ = f (x) is expressed. All calculations -regardless of the choice of 'utilized coordinates' Y -still contain the geometry regarding X. These computations are not invariant with respect to transformations, implying that our approximation is not invariant, too. This non-invariance is a necessity, since especially SIMs are not invariant to transformations as well. Consider the systems with 0 < ε < 1. System (I) can be transformed into (II) by System (I) has a SIM at {x 1 = 0} and (II) at {y 2 = 0} (not invariant regarding φ). Our method -and basically all other methods -can identify both SIMs separately. Covariance means that we can use φ as a local map and {y 1 , y 2 } as local coordinates in order to calculate the SIM of system (I). {x 1 , x 2 } are 'parent coordinates' and {y 1 , y 2 } are the utilized ones. In doing so, we receive a SIM at y 2 = 0 which is incorrect with regard to system (II), but correct in terms of (I) by identifying x 2 = φ −1 (y) = 0 ⇔ y 2 = 0. Conversely, if we choose {y 1 , y 2 } as parent coordinates, we receive different tensors and a SIM at y 1 = 0 (correct for system (II)). In Section 2.3, we demonstrate how tensors and their coefficients change when choosing a different 'parent system' and how a covariant change of coefficients is calculated.

Solution Trajectories as Geodesics in Spacetime
This section briefly introduces the main differential geometric setting of this work. It also discusses the motivation of choosing this specific setting based on geometric observations and physical analogies.

Geometric Motivation
In dissipative multiple time-scale systems solution initially fast trajectories for arbitrary initial values often converge towards an invariant submanifold while slowing down.
In extended phase space, by the introducing time τ as an additional axis, this behaviour is still observed.
is an open set. We call the dynamical system d dt section 2.1 illustrates how bundling behaviors of trajectories the original system relates to bundling of those of the extended one for the two-dimensional linear system d dt The goal of this section is to formulate a setting capable of translating established SIM approximation methods into Riemannian geometry. Riemannian manifolds (M, g) -consisting of a smooth manifold M and a metric tensor field g -are the foundations of that framework. We suggest a Riemannian manifold (M = M f , g = g f ) -which we call f -manifold -depending on a given vector field f and show that the metric g f constructed from f qualifies for an appropriate geometric setting.
A recent work of Heiter and Lebiedz [17] reformulates flow-invariance by vanishing specific time-sectional curvatures of submanifolds of the extended phasespace R n+1 . Our work adopts from [17] the idea of considering the extended phase-space (x, τ ) ∈ R n ×R in order to construct the metric g f with the desirable properties.

Utilization of differential geometry
The field of Riemann geometry offers a wide variety of geometric quantities defined intrinsically, i.e. without reference to an embeeding space. We define a metric g on the open set E × R (which we call M) giving rise to a Riemannian manifold (M, g). A detailed overview of Riemann geometry and curvature can be found e.g. in [31]. Our metric g is chosen in a way that turns every solution trajectory into a geodesic -a shortest connection path with respect to the metric g. All differential geometric quantities used in this work depend on this specific metric g which itself depends on the given dynamical system and is computed from the generating vector field.
We integrate the former ideas into a mathematical formalism and introduce the basic notions of differential geometry. In the following definitions, we always assume that f ∈ C ∞ (E, R n ) for some open set E ⊂ R n and n ∈ N is fixed. We call the first n coordinates of this chart x 1 (p), . . . , x n (p) = p 1 , . . . , p n the state-components. In contrast, the last coordinate τ (p) = p n+1 is the so-called time-component.
Let T p M and T p M denote the tangent space and cotangent space respectively for each point p. The set of derivatives in the direction of each coordinate forms a basis of T p M. These tangent vectors are denoted by The corresponding dual basis consisting of covectors is denoted by dx 1,p , ..., dx n,p and dτ p . For k, ∈ N, T k (M) indicates the set of all k-times covariant andtimes contravariant tensor fields on M. T 1 0 (M) = T M and T 0 1 M represent the tangent bundle and cotangent bundle respectively. We denote the base vector fields ∂ i by The base covector fields dx i , dτ are defined in the same manner. We define the specific metric g used in this work as a tensor field and state its basic properties: defines a smooth tensor field g on M, where ⊗ indicates the tensor product.
For every fixed p ∈ M, g p is represented by its components g ij = g p,ij := The metric g p (as a symmetric bilinear form on tangent space) applied to the tuple (v p , w p ) then can be calculated by The components of g p can be deduced from definition 3 and read where Id n indicates the n × n identity matrix.
Proposition 1. Let M be defined as before and p = [x p , τ p ] ∈ M an arbitrary point. The tensor g p is a metric for every fixed p, independent of the values of Proof. The tensor g p is a bilinear form at every point p by the definition of the tensor product. The symmetry of the matrix (g p,ij ) 1≤i,j≤n+1 in equation eq. (5) implies the pointwise symmetry of g p . It suffices to show that the matrix (g p,ij ) 1≤i,j≤n+1 is positive definite for every value of f (x(p)). Since the identity Id n is positive definite implying that all its minors are positive, we only have to show that det(g p,ij ) > 0. Adding f i (x p )-times the i-th column to the (n + 1)st Hence, the matrix g p,ij is positive definite and g p is a metric tensor for each p ∈ M.
Corollary 1. For any given smooth function f , the tuple (M, g) is a Riemannian manifold.
We call the tuple (M, g) f-manifold. The right hand side of the extended system in definition 1 defines a smooth vector field T : M → T M on M. Its coordinate representation is given by The extended system is a dynamic system on M. The core property of the metric g is formalized in the following theorem: be the corresponding f-manifold and ∇ = ∇ g be the Levi-Civita connection which preserves g. Let γ : (−ε, ε) → M be a solution curve of the extended system of f on M. Then, γ is a geodesic with respect to ∇. In particular, γ satisfies the geodesic equation With regard to the coordinates (x, τ ), the former equality reads with Γ k ij (γ(t)) being the Christoffel symbols of the Levi-Civita connection ∇ evaluated at the point γ(t).
Proof. By calculation of the Christoffel symbols Γ k ij , see Appendix.
Every solution trajectory of the extended system has equal velocity, since In this sense, the metric is a normalizer of the time parametrization of solution trajectories. The proposed setting shares similarities with general relativity where the trajectories of free-falling particles are geodesics with regard to a metric representing a gravitational field. This interpretation motivates the use of concepts from general relativity -such as geodesic deviation -to represent attractiveness. This is new in the SIM context and the idea is implemented in Section 3.

Covariant Transformation of the metric
We briefly demonstrate how the coefficients of the metric tensor g from definition 3 change under simple re-scaling of parent coordinates. Letẋ = f (x) be the parent system with coordinates X := {x 1 , . . . , x n , τ }. The metric tensor from the last section regarding this choice of parent coordinates is denoted by [g] (1) .
Its coefficients regarding X are denoted by g In Y −coordinates, the phase-space system then readsẏ k = (a k f (x(y))) k=1,...,n . We can express [g] (1) by the means of Y and calculate Inserting the transformation 8 and original coefficients from definition 3, we receive Alternatively, when declaringẏ k = (a k f (x(y))) k as parent system, we receive a different metric [g] (2) . The coefficients have to be chosen according to definition 3. Hence, the resulting coefficients with respect to Y are The metrics are different, thus all derived tensors as well as our calculated manifold.

Geodesic Stretching Approach
In this section, we introduce a new tensorial method to approximate normal attractiveness by exploiting the previously introduced setting in the following way: We translate the notion of geodesic deviation to the concept of the ration of stretching rates from SBD (see [1,2]). The results are what we call geodesic stretching rates which turn out to be specific sectional curvatures, in coordinates corresponding to some curvature tensor entries. In 3.4 we apply the resulting method to non-linear test-models. In 3.4 we derive a SIM approximation method from this approach and apply the resulting method to non-linear test-models.

Deviation
In general relativity geodesic deviation is used to describe relative behavior of neighboring geodesics corresponding to the relative acceleration of nearby particles in free-fall. It is defined by plugging in a tangent vector y p -representing the instantaneous velocity of the geodesic -into the first and third argument of the Riemann curvature tensor which is denoted by The result is a tensor field depending on the tangent vectors y p . The input vector v p of this reduced tensor represents a small displacement between the neighboring geodesics, while the output stands for the deviation. On the fmanifolds (M, g) from the previous Section 2, there is one set of geodesics of special interest: The solution trajectories of the extended system bundling near a NAIM. Hence, an evident choice is y p = T p for all p ∈ M and receive tensorfield depending on T p , leading to the following definition: Let R = R f be the corresponding Riemann curvature tensor. We call the tensor field Remark 1. The christoffel symbols, the curvature tensor and the f-deviation do not depend on the explicit time τ .
Proof. The components of g are independent of the time τ , implying that the components g ij of the inverse metric tensor are also time-independent. Since we use the Levi-Civita connection, the Christoffel symbols are calculated by derivatives of time-independent quantities and the statement holds for Γ k ij . Using the same argumentation, we conclude this property for the curvature tensor R and the f-deviation S.
Based on its properties and geometric interpretations, the f-deviation appears to be well-suited to be turned into a geometric criterion to approximate normal attractiveness. We now aim to deduce a scalar, curvature-based quantity from the f-deviation that intuitively represents the bundling behavior. In order to do so, we are guided by an existing, geometric approach to characterize normal attractiveness: The SBD, introduced in [1], [2] by Adrover et al. .

Original Stretching approach
SBD in dissipative and chaotic system is a local, geometric reduction approach to multiple time scale dynamics. Let M be an embedded submanifold of R n with x being the euclidean coordinates and p ∈ M an arbitrary point. Normal attractiveness requires vector dynamics in tangential direction to be dominated by the dynamic in normal direction. This property is encoded in the generalized Lyapunov-type number (see [21] and [40]) Here, v t = v(x(t)) and w t = w(x(t)) represent the vector dynamic along a solution of the ODEẋ = f (x) for tangential and normal initial vectors respectively. Combined with demanding v −t to approach infinity at a certain rate, normal attractiveness requires σ(p) < 1. Hence, we compare the convergence rates of the limits lim t→∞ v −t and lim t→∞ w −t .
According to [2], an evident local quantity to approximate contraction rates of vector dynamics in time, is to take the first derivative of the square: where the brackets ·, · represent the euclidean inner product. Normalizing the previous term by division by 2 v 2 -implying that the resulting terms do not depend on the vector's length -leads to the definition of so-called stretching rates ω p (T ) and ω p (N ) for tangential and normal vectors respectively, where Instead of considering the limit of the quotient ( w −t / v −t ), Adrover et. al. take the quotient of the stretching rates instead. These are so-called stretching ratios r(p) = ω p (N )/ω p (T ) and replace the generalized Lyapunov-type numbers σ(p). We receive local geometric quantities representing NAIM properties: Normal stretching ω p (N ) represents bundling behavior near M , tangential stretching ω p (T ) indicates acceleration/deceleration alongside M , while their ratio r(p) locally approximates normal attractiveness. An evident way of introducing a lower-dimensional manifold for model-order reduction is to maximize/minimize the stretchting rates and/or the corresponding ratio.

Geodesic stretching
The stretching rates incorporate a geometric interpretation which is adopted to be transferred in the differential geometric setting from Section 2. We can interpret (R n , ·, · ) as a Riemann manifold, equipped with euclidean metric.
The term J f v represents the so-called vector dynamics of the flow φ(t) related to the differential equationẋ = f (x). The flow differential Dφ t : T x(0) R n → T x(t) R n propagates perturbation vectors v x(0) along φ t . Hence the map Dφ t indicates how solution trajectories ofẋ = f (x) diverge or converge. Aiming to extract a local rate of deviation, a direct calculation yields Hence, the mapping J : T x R n → T x R n , v → J f v assigns a perturbation vector v to a local rate of deviation. This interpretation shares major similarities with the f-deviation introduced in definition 4. An evident transfer of the ω x (v) into the coordinate-free setting from Section 2 is to replace (R n , ·, ·, ) by (M, g) and apply the f-deviation instead of J. The result is the so-called geodesic stretching rate: Definition 5. Let (M, g) be defined as in Section 2, S the f-deviation, p ∈ M an arbitrary point and v p ∈ T p M. The mapping is called geodesic stretching. The image ϑ p (v p ) is denoted as geodesic stretching rate of v p . By definition, we can calculate geodesic stretching rates for every tangent vector v p ∈ T p M, especially those with non-vanishing time component dτ p (v p ) = 0. We exclude the latter ones and only consider those of the following subspace: Definition 6. Let p ∈ M be arbitrary. We call the subspace This restriction appears natural, since the time axis τ is included artificially and we have no valid interpretation for "perturbation in explicit time". A visualization of the geodesic stretching rates can be found in section 3.3.
There exists a curvature-based correspondent ϑ p (v p ) that is well-defined for v p ∈ T p,s , formalized in the following theorem: Theorem 2. Let (M, g), the stretching rates ϑ p , the subspace T p,s defined as above, v p ∈ T p,s arbitrary. Then, ϑ p (v p ) equals the sectional curvature of the subspace σ v p spanned by the vectors T p and v p .
Proof. Let p ∈ M and v p ∈ T p,s be arbitrary. The set {∂ 1,p , . . . ∂ n,p , T p } is an orthonormal basis of T p M for all p ∈ M with respect to g p , implying g p (T p , v p ) = 0. Using the fact that g p (T p , T p ) = 1 we can calculate where σ v p := span(T p , v p ) and K p (σ v p ) is the sectional curvature of σ v p .
Corollary 2. The geodesic stretching rate ϑ(v p ) is a covariant intrinsically geometric quantity for every v p ∈ T p,s and p ∈ M.
The aim is to exploit the rates ϑ(v p ) in an analogous way as the stretching characterize NAIMs by decomposing the tangent space into tangential and normal directions of a submanifold U ⊂ R n . Since M is a space-time manifold, we select the pure-state space T p,x and split it in the same manner as in the previous subsection. The result is the following definition: Definition 7. Let x ∈ R n and U ⊂ R n be an embedded submanifold with x ∈ U . Assume dim(T x U ) = k and vectors u 1,x , . . . u n,x ∈ T x R n satisfying Let [u 1j , . . . , u nj ] be the euclidean coordinates of u x,j for all j ∈ {1, . . . , n}, p = [x, τ ] ∈ M for an arbitrary τ and x from above. We then define the projected tangent space T tan p,s and projected normal space T orth p,s by Definition 8. Let U be a submanifold of R n , T tan p,s (U ) and T orth p,s (U ) defined as above. We define the tangential and orthogonal stretching rate as respectively.
Geodesic stretching rates are now used similarly to the original stretching rates in [2]: Tangential-and normal geodesic stretching rates measure deceleration and contraction respectively, while their ratio approximates the property of normal attractiveness.
A consequence of the geodesic stretching rate construction is the fact that it is capable of providing an approximation for a specifically constructed NAIM: the SIM in the context of slow-fast systems which is subject of the following section. We call the resulting method Geodesic Stretching method (GSM).

Testing the Geodesic Stretching method
We conclude this section by applying the above method non-linear test models. We go into more detail on how this construction is deployed to approximate one-dimensional SIMs in these two-dimensional models. Finally, we discuss how the approach can then be applied to more general cases. Consider the non-linear Davis-Skodje system (see [7]) where the parameter η > 1 measures time-scale separation. This system has a one-dimensional SIM with graph representation  Hence, both the tangential and normal space of each trajectory is onedimensional as well. Following definition 7, we get the subspaces Since the subspaces are one-dimensional and the geodesic stretching rate does not depend on the length of each vector, we get Θ tan section 3.4 depicts tangential and respectively orthogonal stretching rates in the vicinity of the SIM. Very close to the SIM, the tangential stretching rate gets small, while the orthogonal one gets particularly large. This observation matches the description of the SBD, introduced in the beginning of this section. By fixing one variable (e.g. x 1 ) and maximizing/minimizing orthogonal/tangential geodesic stretching rates with respect to the other variable, we get at least an adequate approximation of the SIM. Figure 3.4 shows that the former criterion is not exact. There, we fix x 1 = 1 and consider both tangential and orthogonal geodesic stretching as a function of x 2 . Both resulting one-dimensional graphs have an extremum at around x 2 = 0.4985, while the SIM point is at 0.5. We directly conclude that the ratio between both rates is also extremal at around x 2 = 0.4985.

Comparing the Geodesic Stretching Method
In this subsection, we refer to the chemical reaction mechanism taken from [6] aiming to compare our method to different well-estabished methods to approximate the SIM. In this model, we have four species with respective concentrations c 1 , c 2 , c 3 and c 4 . Apart from the ILDM, the following SIM ap-  Figure 5: Illustration of different SIM approximation methods applied to the test system defined by (11) and (12) respectively, similar to Figure 4(a) in [6]. A SIM is given by the vertical line c 3 = 0.1.
proximation methods are applied in [6]: • Quasi-Equilibrium-Manifold (QEM) • Spectral-Quasi-Equilibrium-Manifold (SQEM) • Symmetric-Entropic-Intrinsic-Low-Dimensional-Manifold (SEILDM) Considering the conservation law and choosing a specific set of constants for that mechanism (see [6] for more details), we receive the two-dimensional system We apply our Geodesic Stretching Method (GSM) to this system. In this case, we have a one-dimensional SIM -according to the GSPT-definition. The SIM is the vertical line c 3 = 0.1 with equilibrium point (c 1 = 0.5, c 3 = 0.1). In Figure  3.4 you can see the corresponding approximations plotted in the (c 3 − c 1 )-plane. Applying the GSM, c 1 is our reaction-progress variable.
We can see that the approximation error of the QEM and the SQEM is comparatively large. Both SEILDM and GSM provide significantly better approximations with similar deviations from the SIM. One of two candidates provided by the solution of the ILDM equation coincides exactly with the SIM which is unsurprising, since the SIM is linear and the approximation error of the ILDM is proportional to the curvature of the SIM [19].

General applicability of Geodesic Stretching
The previous two test models are comparatively simple in the following manner: Both are (a) two-dimensional with a one-dimensional SIM, (b) slow-fast systems, (c) formulated on the (linear) phase-space R 2 .
We briefly point out how this approach can be applied to systems that are not limited to the properties (a) -(c) in our test examples, it can be used for systems of every (finite) dimension n with a SIM of an arbitrary dimension 1 ≤ m < n (addressing point (a)). Definitions 7 and 8 already cover how this more general case is tackled: Consider the subspace T tan p,s (U ) of the whole tangent space that is tangent to a potential SIM U and calculate the maximum stretching rate for each tangent vector of this subspace (denoted as Θ tan p (U )). This generalization is still covariant since every tangent vector is a tensor with a coordinate-independent meaning.
In case of a general systemẋ = f (x), the object we want to approximate is a NAIM, as the notion of a SIM is tied to slow-fast systems. The geodesic stretching approach is formulated for a system of the formẋ = f (x) and approximates normal attractiveness. Hence, there is no limitation to slow-fast systems and Restriction (b) is addressed by construction.
Concerning point (c), imagine a model containing adiabatic constraints g(x) = 0 where g : R n → R with 1 ≤ < n. Instead of operating on an open subset of R n , the dynamic is now formulated on a manifold. Conveniently, all exploited notions and tensors used to define this geodesic stretching approach are formulated on Riemann manifolds. We only need to choose a parent coordinate system and define the metric on this manifold according to definition 3. This approach can be applied directly, once this is done.

4
Differential Geometric Interpretation of the Flow Curvature Method Sections 2 and 3 cover a new covariant approach on manifolds to approximate NAIM's core property by utilizing a differential-geometry setting. In this section, we aim to investigate a potential covariant reformulation of an existing geometry-based approach to SIMs: The flow-curvature method (FCM) by Ginoux [14]. We begin by briefly stating its main properties:

Flow Curvature Method in a nutshell
The foundation of this ansatz are higher curvatures of trajectories in the phase space R n . In case of an n-dimensional dynamical systemẋ = f (x) with (n − 1)dimensional SIM, its FCM-approximation is defined by the union of all points p with vanishing n-th curvature of the trajectory. It can be shown that the FCM is capable of approximating a SIM to order n [15]. The former criterion is satisfied if and only if Its solution is called flow curvature manifold which is not flow-invariant, as long as d dt J f (x(t)) does not vanish. This can directly be seen by calculating d dt (Ψ(x(t))) .
On the other hand, this manifold is also non-invariant regarding coordinate transformations. Away from each fixed point, we can locally transform the system into a constant system, e.g.ẏ = g(y) ≡ c, with y = Φ (c) (x), for a given systemẋ = f (x). In these new coordinates, every point satisfies the flow curvature criterion. This argument also proves that the FCM is not covariant. The value of Ψ -and crucially -whether or not Ψ vanishes depends on the coordinate chart, as Φ (c) shows.
An evident way of finding a covariant reformulation of the FCM is to take steps similar to the approach used to develop the geodesic stretching method.
Step one: Translate the existing method into a manifold-based setting.
Step two: Modify the method in a sensible way, such that the scalar value (here: Ψ) is the evaluation of a tensor and can be expressed in any appropriate coordinate chart. We implement the first step in the following subsections.

Flow derivatives as covariant derivatives in euclidean space
Consider the manifold (M, g e ) = (R n , ·, · ) where ·, · = g e represents the euclidean inner product at each point p ∈ M. Let ∂ j,p ∈ T p M indicates the tangent vector in the direction of the j-th coordinate. Let dx j,p represent the dual basis on each point p, we receive g (e) p = n j=1 dx j,p ⊗ dx j,p ∀p ∈ R n .
The christoffel symbols of the Levi-Civita connection all vanish.
By definition, both the gramian matrix and determinant are coordinate independent for every metric g. In case of the euclidean metric g e , a direct calculation shows Using the multiplicativity of the determinant, we conclude that the definition criterion for the FCM can be written in the following coordinate independent way: Φ(p) = 0 ⇔ D p (∇ The last transformation finishes the embedding of the FCM into the Riemann geometry framework. All translated quantities yield are well-defined within this field. Unfortunately, the utilized notions are not tensors -and as mentioned in the section 4.1 -a modification is now needed. One might be guided by the geometric interpretation of the FCM: Inspect the highest curvature of a solution trajectory. Find a way to translate this interpretation into a differential geometry framework to define the appropriate tensors.

Summary, Conclusion and Outlook
The aim of this work was to present a route to formulate tensorial representation of normal attractiveness in multiple time-scale systems, enabling tensorial approximations for slow invariant manifolds for slow-fast systems. For this purpose we introduce a general differential geometric setting in Section 2. We exploit the notion of intrinsic curvature to reformulate the stretching-based analysis. This covariant formulation on manifolds makes this work a novelty in this context. We exemplarily apply the resulting approach to the Davis-Skodje system and the Michaels-Menten model in Section 3.4. In Section 4 we also reformulate the flow curvature method by expressing its utilized flow derivatives by covariant ones. Our ideas might be useful as a general guideline towards finding tensorial reformulations of established SIM methods. The authors share the opinion that the field of differential geometry is an appropriate frame as to express essential SIM quantities intrinsically as e.g. shown in [17].
where the entries marked by a * are determined by the symmetry Γ k ij = Γ k ji . In case k = n + 1, the components we receive are: Thus, the covariant derivatives of the base vector fields ∇ ∂i,p ∂ j vanish. Let h ∈ T R n be a smooth vector field with components h k (x). We receive h k ∂x i ∂ k,p ∀i = 1, . . . , n.
By linearity we conclude that The component vector of the flow derivatives is calculated by dh(x(t)) dt ∂h k ∂x i which proves the Lemma.