Sorting-free Hill-based stability analysis of periodic solutions through Koopman analysis

In this paper, we aim to study nonlinear time-periodic systems using the Koopman operator, which provides a way to approximate the dynamics of a nonlinear system by a linear time-invariant system of higher order. We propose for the considered system class a specific choice of Koopman basis functions combining the Taylor and Fourier bases. This basis allows to recover all equations necessary to perform the harmonic balance method as well as the Hill analysis directly from the linear lifted dynamics. The key idea of this paper is using this lifted dynamics to formulate a new method to obtain stability information from the Hill matrix. The error-prone and computationally intense task known by sorting, which means identifying the best subset of approximate Floquet exponents from all available candidates, is circumvented in the proposed method. The Mathieu equation and an n-DOF generalization are used to exemplify these findings.


Introduction
The objective of this paper is to introduce a novel stability method based on the Hill matrix, which differs from the state-of-the-art methods in that a matrix projection is applied before computing an eigenvalue problem. The structure of this projection is obtained by considering the Hill matrix to be a result from the Koopman lift for well-chosen basis functions.
The Koopman framework [1,2] has gained immense popularity in recent years as a versatile tool for various engineering applications, such as system identification [3], model order reduction [4] and feedback control [5]. This is due to an auspicious promise: global representation of a nonlinear system by a linear operator. To this end, in the Koopman framework, the dynamical system is defined through the propagation of functions on the state space, also called observables, over time. Bernard Koopman first described a unitary linear operator which evolves a class of measurable functions along the flow of a conservative system [6], which would later be named the Koopman operator. After some relatively quiet years, the Koopman operator experienced a revival around the turn of the century, when it was shown that its spectral characteristics contain global properties for the underlying dynamical system [7]. This sparked generalizations to nonconservative systems [8,9] and at the same time, numerical and data-driven methods like the Arnoldi method [10] and extended dynamic mode decomposition [2] emerged.
Classically, the Koopman framework is applied to time-autonomous systemsẋ = f(x) and the approximate linear dynamics obtained by the Koopman lift then takes the formż = Az. The incorporation of a time-dependent input v(t) into the dynamics, i.e., x = f(x, v(t)) or simplyẋ = f(x, t), generally poses problems in the Koopman framework as the system can only be approximated by a linear time-invariant (LTI) systemż = Az + Bu(t) if products of state and input are neglected. In this paper we focus on non-autonomous systems for which the input is timeperiodic, i.e., v(t) = v(t + T ). In particular, we propose a specific choice of observable functions which contains observables depending both on state and time, opening the possibility to include products of state and input.
The numerical computation of periodic solutions in time-periodic non-autonomous systems is a task of greatest interest in engineering application. Periodic solutions are of prime importance for, e.g., nonlinear vibration analysis in structural dynamics [11,12], in acoustics [13] and thermo-acoustics [14], nonlinear oscillation analysis in electronic circuits [15] as well as heart rhythm analysis in cardiology [16]. Therefore, it is important to find and characterize periodic solutions, assess their stability properties and also track these quantities along varying system parameters, a process which is called continuation [17].
Naively, attractive periodic solutions can be found by simply simulating numerically over a long time interval, until transient effects are negligible. However, since this method is very dependent on the initial condition, can only find attractive solutions and is computationally expensive, more sophisticated methods have been developed. There are a multitude of periodic solution solvers available, including finite differences [18], shooting [19], multiple shooting [20], collocation and generalized collocation [21] and harmonic balancing [22]. We highlight two methods: • The shooting method [19] uses the Newton method and a numerical ODE solver to find initial conditions x 0 such that x(t 0 + T ) = x(t 0 ) = x 0 , i.e., whose trajectory over a period fulfills the periodicity constraint. Since the monodromy matrix appears in the update step of the Newton method, stability information comes (almost) for free with this method. This monodromy matrix also plays a role in the tangent prediction in an arc-length continuation method. • The harmonic balance method (HBM) [22,23], in contrast to the other mentioned methods, is a purely frequency-based method. The periodic solution is parameterized globally by a finite set of trigonometric functions, and the coefficients are determined using a residual in the frequency domain. The transfer of nonlinear terms into the frequency domain is non-trivial. In practice, this is achieved either using an alternating frequency and time (AFT) method [24] or by providing a recast of the system in quadratic form [25]. One advantage of the HBM is that it automatically provides a filtering effect on the identified periodic solution.
While stability information about the periodic solution branches comes almost for free in the shooting methods as the monodromy matrix is calculated as a necessary continuation step, stability information and hence information about the bifurcations that occur in a system are not so trivial in the HBM. The Hill matrix provides a frequency-based method to assess the stability of periodic solutions [26,27]. It can be found and constructed easily in the numerical asymptotic method (ANM) [25,28], being a continuation method based on the HBM. Hence, computation of stability through its eigenvalues, which approximate the Floquet exponents, is a viable option.
However, two critical problems make the Hill method often unattractive in practice. On the one hand, from a numerical viewpoint, computing the eigenvalues of the large Hill matrix is computationally expensive and potentially inaccurate [29]. On the other hand, for correct assertion of stability, only a non-trivial subset of these eigenvalues must be considered. This process is known in the literature as sorting of Floquet exponent candidates. The determination of this subset is still an active area of research, with the approaches being based on the imaginary parts [23,26,30] and potentially in addition the real parts [31] of the eigenvalues, or alternatively symmetry considerations of the eigenvectors [27,28].
In this work, we propose a different approach for obtaining stability information from the Hill matrix, which circumvents both issues mentioned above. Using the Koopman framework we motivate a novel dynamical systems interpretation of the Hill matrix, which allows to compute an approximation of the monodromy matrix directly (i.e., without computing a large number of eigenvalues and subsequent sorting). The proposed method to find the monodromy matrix from the Hill matrix involves the action of the matrix exponential of the Hill matrix applied to a smaller sparse matrix, followed by a projection to the n × n monodromy matrix. Finally, the stability of the periodic solution can directly be assessed from the n eigenvalues of the monodromy matrix, known to be the Floquet multipliers.
Parts of the research in this work were presented in a preliminary form at the ENOC2020+2 Conference [32], in particular concerning the proposed choice of the Koopman basis functions in Sect. 3 and the connection to the harmonic balance equations and the Hill matrix. The main (and original) contributions of this work are the novel stability method of Sect. 4, the considerations with respect to the matrix projection and the formal proofs for the theorems in Sect. 3.2.
The paper is structured as follows. Section 2 provides an overview of the notation and gives a theoretical background for the concepts that are central to this work, in particular concerning a selection of topics from Koopman theory as well as frequency-based methods for periodic systems. Section 3 introduces the chosen basis for a Koopman lift on time-periodic systems and states the three central theorems which relate this Koopman lift to the classical frequencybased methods. The proofs of these theorems can be found in Appendix B. Section 4 presents the novel stability method based on the findings from Sect. 3, and the projection to the monodromy matrix as well as the computational effort are discussed. These results are illustrated in Sect. 5 using numerical investigations on two exemplary dynamical systems. Finally, concluding remarks are given in Sect. 6.

Theoretical background
In this section, the reader is provided with an overview over the theory in the Koopman framework and the Floquet theory that is necessary for the later parts of the paper.

Notation and terminology
The frequency-based methods considered in this work rely heavily on being represented as (generalized or classical) Fourier series. Hence, a short overview over Fourier series and the notation that is employed will be given. Scalar quantities will be represented by Greek or Latin slanted lower case letters. This includes scalarvalued functions as elements of a function space. Vectors in Euclidean space will be represented by Latin bold lowercase letters and matrices will be represented by Latin bold uppercase letters. Tuples of functions are represented by bold font and the distinction to Euclidean space can be drawn from context.
The choice of index is connected to its meaning. The index l ∈ N is used for indexing over the states of a system, whereas the index k ∈ Z is used for indexing over frequency harmonics. While j may appear as arbitrary auxiliary index, the letter i usually denotes the imaginary number i 2 = −1 and is only used for indexing purposes if the indexing context is obvious. The Kronecker delta is denoted by δ jl .
If an inner product ·, · with the usual inner product properties (see, e.g., [33]) is defined on a vector space F, elements f, g ∈ F of this space are orthogonal if their inner product is zero. An orthonormal system is a set ξ j D j=1 , D ∈ N ∪ {∞} with ξ i , ξ j = δ i j , and a maximal orthonormal system whose span is a dense subset of the considered vector space is called an orthonormal basis. By slight abuse of notation, the inner product is extended in this paper to tuples of functions g ∈ F l , h ∈ F m element-wise via It can be easily verified that, as a generalization of the conjugate symmetry of the inner product, the relation holds, where U * denotes the conjugate transpose of the matrix U. Moreover, the sesquilinearity of the inner product extends to matrices via for constant complex matrices U, V of appropriate size.
If a space F has an orthonormal basis ξ j D j=1 stacked into a tuple ξ , it is well-known [33] that elements g ∈ F admit a (generalized) Fourier series where the matrix inner product notation was used in the last term. The constant, scalar coefficient g, ξ j is called the j-th Fourier coefficient of g. In particular, in the space of trigonometric functions of period T , the functions u k : [0 T ) → C, t → e ikωt , k ∈ Z constitute an orthonormal basis with respect to the inner product This is the classical Fourier series. For a finitedimensional subspace spanned by finitely many basis functions {u k } N u k=−N u , a function g can thus be expressed by with the matrix-valued inner product notation introduced earlier.

Koopman theory overview
A short introduction to the Koopman framework to set the notation and help the reader understand the following sections is given below. It is, however, not intended to give a comprehensive understanding of the current overall state of the art in the field that would be suitable for a wider range of application. For such a more general and in-depth treatment of the considered methodology, the authors rather recommend [10,34]. Consider a non-autonomous time-periodic finitedimensional dynamical system governed bẏ where t ∈ R is the time, x(t) ∈ X ⊆ R n is the state trajectory starting at x(t 0 ) = x 0 and f : X × R → X is a smooth vector field which is T -periodic in t. The family of maps φ t (x 0 , t 0 ) = x(t) characterizes the flow of the system and assigns to each (initial) configuration (x 0 , t 0 ) the resulting state at time t ≥ t 0 . The Koopman framework [10] considers output functions g(x, t), also called observables. Any Banach spaces of functions over the complex or real numbers are permitted in the general Koopman framework. In this work, we consider in particular the space F of complex-valued functions g : X × R → C which are real analytic on X and T -periodic in the last argument t. Given any function g, it may be of interest how its function values evolve along the trajectories of the system. For instance, in the Lyapunov framework, it is desired that function values of a Lyapunov candidate decrease over time for any starting point. The operator K t : F → F; g → g • φ t performs this shift along the trajectory for arbitrary functions g from the considered function space. The family of all these operators for any t is called the Koopman semigroup of operators. Indeed, the semigroup properties with respect to the time parameter can be verified easily.
For suitable function spaces F, this Koopman semigroup contains all information about the system without explicitly knowing the vector field f or the flow φ t . In particular, if F is chosen such that the identity function id is contained in the vector space, then the flow can be recovered easily by simply evaluating K t (id). As a trivial counterexample, consider the one-dimensional vector space of constant functions. Any constant function will not change its function value while being evaluated along an arbitrary trajectory of arguments. Therefore, in this case, the Koopman operator semigroup is well defined, albeit trivial. No information about the underlying system is retained. This example shows that an appropriate choice of function space is a crucial part of the Koopman framework.
Under the aforementioned assumptions for the particular function space F and the vector field f, the Koopman semigroup is continuous with respect to time and there also exists the operator L : , mapping an observable g to its total time derivativeġ along the flow witḣ The operator L is called the infinitesimal Koopman generator. Again, for suitable F, this representation alone is a sufficient way to describe the behavior of the Fig. 1 Schematic drawing of the infinitesimal Koopman generator L, its finite-dimensional approximation LN and the Koopman liftÂ dynamical system. In particular, if id ∈ F, the vector field f is easily recovered. In addition, the infinitesimal Koopman generator and the Koopman semigroup of operators are linear in the argument g, even if the governing differential equation is nonlinear. This comes at the cost of dealing with a mapping on an (infinite-dimensional) function space F instead of the (finite-dimensional) state space X .
For practical reasons, we will be forced to project the dynamics on a finite-dimensional subspace FN ⊂ F spanned byN linearly independent basis functions ψ j N j=1 . Any projection N : F → FN defines a finite-dimensional approximationL : FN → FN of L on FN by LN := N L. The approximation process and the subsequent approximation error are visualized in Fig. 1. As the subspace FN generally is not closed w.r.t. L, the result of Lg must be projected back onto FN , introducing some approximation error. As for F, the choice of the finite space FN and the projection onto it is crucial.
To evaluate a system trajectory in the usual statespace form, the initial condition is fixed first, the state evolution is computed afterward and the output function is computed last. For the Koopman infinitesimal generator (and its approximation on the finitedimensional function space), this is reversed. First, an observable, or output, is fixed, the observable is propagated along the flow and the initial condition is determined in the last step by evaluating (Lg)(x) for an initial condition x. To arrive back at a linear system representation in the usual state-space form, this behavior must be reversed again. This is achieved by evaluating the action of the (approximate) Koopman infinitesimal generator on the basis functions. Settinĝ and keeping this dynamics for increasing t, the linear autonomous systeṁ results. This linear system is called the Koopman lift and describes the dynamics represented in LN . With the lifted stateŝ it is a finite-dimensional linear approximation of the original system dynamics. The Koopman lift matrixÂ can be derived from the original nonlinear dynamics manually using a column vector := ψ 1 , . . . , ψN T of the basis functions of FN , computing d dt and identifying terms linear in elements of after applying the projection N from F onto FN . If this projection is orthonormal and is an orthonormal system, then the matrix entryÂ i, j at i-th row and j-th column is given byÂ i, j = dψ i dt , ψ j , where ., . denotes the corresponding inner product.
Depending on the basis structure chosen, various popular embedding techniques emerge as a Koopman lift for specific system classes. For instance, if a monomial basis is chosen for an autonomous system, the Carleman linearization [35] results as Koopman lift. For smooth, polynomial systems, this is often the first choice of basis dictionary [36]. Alternatively, delay coordinates are often employed in Koopman-based applications [2]. Based on the Takens embedding theorem, this can capture weakly nonlinear dynamics [37]. For periodic systems, the Fourier embedding [38] has been known, although its properties have mainly been analyzed in the frequency domain.

Harmonic balance method
Consider the non-autonomous time-periodic finitedimensional dynamical system (6) as above. Often one is interested in finding a T -periodic solution, i.e., a solution to the dynamical system (6) which fulfills x(t + T ) = x(t) for all t ≥ 0. This constitutes a boundary value problem (BVP) and there are methods to solve this type of BVP in the time domain and in the frequency domain. Shooting, multiple shooting and collocation methods all rely on an interplay between time-integration (or finite differencing) of the ODE and solving nonlinear functions for periodicity and continuity constraints [21].
In contrast, the HBM is a frequency-based method. Under suitable smoothness assumptions, the periodic solution has a convergent Fourier series. Hence the periodic solution can be approximated by its Fourier expansion up to order N HBM with unknown parameters via with ω = 2π T , u(t) = (e −i N HBM ωt , . . . , e i N HBM ωt ) being a vector of Fourier base functions and p −N HBM , . . . , p N HBM gathering the corresponding (unknown) coefficients. These coefficients p k are then determined by substituting this ansatz into the system equation (6). The comparison of coefficients for the Fourier expansions of dx p dt from the definition (10) and f(x p , t) for every order up to N HBM transforms the BVP into a system of n(2N HBM + 1) algebraic equations. Existence and convergence of these HBM approximations has been shown [22]. While the left-hand side of the equation as well as linear terms in f are easy to handle, the frequency component of the nonlinear terms can usually not be expressed in closed form. The individual equations for each order are thus usually determined and simultaneously solved using the fast Fourier transform with an alternating frequency and time (AFT) method [24]. The equations for each order can also be isolated by projecting onto the corresponding basis function from the collection in u through the classical inner product (4). Hence, the HBM approximates a periodic solution by solving the n(2N HBM + 1) algebraic equations collected in With this notation, the numerically cumbersome task of calculating the Fourier coefficients of the nonlinear components of f is hidden in the definition of the inner product.

Floquet theory: stability of periodic solutions
When a periodic orbit x p is found (via HBM or by other means), the next interesting question is that of its stability properties; that is, whether trajectories that start sufficiently close to the periodic orbit will approach it, stay in a vicinity of it or tend away from it with increasing time. To evaluate the stability properties, the dynamics of a perturbation y = x − x p from the periodic solution is considered. Substitution of this definition into the original system dynamics yieldṡ where is the Jacobian of the system evaluated along the periodic solution. The approximate linear time-varying (LTV) systeṁ has an equilibrium at zero, which corresponds to the periodic orbit of the original system, and the stability analysis of the periodic orbit reduces (in the hyperbolic case) to the stability analysis of this equilibrium. This will be the convention for the remainder of this paper, unless stated otherwise. The fundamental solution matrix (t) is the solution to the variational equatioṅ and any state can be obtained via y(t) = (t)y 0 . In particular, the fundamental solution matrix (T ) =: T evaluated after one period is called the monodromy matrix of the system and its eigenvalues {λ l } n l=1 are called Floquet multipliers [39]. The Poincaré map y k+1 = T y k provides snapshots for the evolution of the perturbation y, spaced at a time distance of T , i.e., y k = y(kT ). For the long-term behavior, it is sufficient to consider the evolution of these snapshots. Therefore, stability analysis of the periodic solution reduces to stability analysis of the Poincaré map. Hence, if all Floquet multipliers are of magnitude strictly less than one, the equilibrium of the perturbed LTV system and thus the periodic solution of the original system are asymptotically stable; if at least one eigenvalue has a magnitude strictly larger than one, they are unstable. If there exist Floquet multipliers with magnitude equal to one, but none with a magnitude larger than one, the equilibrium is non-hyperbolic and further investigation is necessary to give conclusive statements about stability of the originally considered periodic solution.
Alternatively to the Floquet multipliers, the stability properties of a time-periodic linear system can be characterized by the Floquet exponents. In the linearized perturbed system (13a), if the matrix T is diagonalizable, there exist n solutions y l (t) = p l (t)e α l t which form a basis of the solution space, where each function p l is T -periodic [39]. Hence, stability is characterized by the real parts of the Floquet exponents {α l } n l=1 . If at least one Floquet exponent lies in the open right half plane, i.e., if at least one real part is larger than zero, the equilibrium is unstable. The Floquet multipliers can be determined by substituting t = T in the Floquet solution and it follows that In contrast to the Floquet multipliers, the Floquet exponents are not uniquely defined. It is easy to see that if the pair (p l (t), α l ) generates a solution y l (t), the same solution is generated by Hence, in total, there are infinitely many valid Floquet exponents, which can be categorized into n distinct groups. All elements of one group have the same real part and differ in the imaginary part by multiples of iω. As stability is determined by the real part only, it is sufficient for stability analysis to know any one element from each of the n groups. All elements of one group map to the same Floquet multiplier. When a periodic orbit is determined using the purely time-domain-based shooting method, the monodromy matrix usually is a direct byproduct of the continuation method [17]. In this case, the numerically obtained monodromy matrix can be evaluated directly to obtain the Floquet multipliers and their stability information.
When the HBM is computed in the standard way, however, stability information about the identified limit cycle is unclear without further investigation. The Hill method [22,40] offers a frequency-domain-based way to approximate the Floquet exponents of the linearized perturbation equation.
The Floquet exponents are eigenvalues of the infinite Hill matrix H ∞ [30], which is constructed from the Fourier coefficients of the periodic system matrix J(t) = ∞ k=−∞ J k e iωkt and reads as Introducing a vector v ∞ of appropriate (infinite) length, the eigenproblem can be formulated. This infinite-dimensional problem has infinitely many discrete eigenvaluesα, which solve (17). They correspond identically to the Floquet exponentsα for all k ∈ Z as introduced above and can be sorted into n groups, where the entries of each group differ by multiples of iω [30]. However, in practice, only the eigenvalues of a finite-dimensional matrix approximation of H ∞ can be computed numerically.
The matrix of size n(2N u + 1) × n(2N u + 1) consists of the n(2N u + 1) most centered rows and columns of H ∞ and approximates the original infinite-dimensional Hill matrix. In the absence of truncation error, the eigenvalues of H would be a subset of the eigenvalues of H ∞ , i.e., the Floquet exponents. Due to the inevitable error which generally comes with truncation, however, this does not quite hold. The N eigenvalues of H, which do not identically coincide with Floquet exponents, will be called Floquet exponent candidates below.
The matrix H has a block Toeplitz structure except for the middle diagonal, and, for sufficiently large N u , the bands near the diagonal dominate as the Fourier coefficients of J tend to zero. Loosely speaking, some eigenvalues affiliated most with the central rows of H are less impacted by the truncation and provide a better approximation to the Floquet exponents than others [28]. Up until the change of the last century, this property was neglected and stability of the periodic solution was asserted based on the real parts of all Floquet exponent candidates [40]. This naive Hill method without any additional steps would often assert instability for stable solutions due to spurious Floquet exponent candidates without physical meaning, giving it a reputation of being inaccurate [22,29].
However, the accuracy of the Hill method can be improved significantly if only a subset of Floquet exponent candidates is considered, instead of all of them. Hence, the search for a selection criterion which determines the best approximation to the Floquet exponents from the Floquet candidates has received much attention in the literature [23,26,28,30].
For sufficiently large N u , it is proven that the candidates with minimal imaginary part in modulus converge to the true Floquet exponents [26]. Since the convergence may only occur for very large truncation orders, an addition to this method was very recently proposed [31]. This modified method first sorts the Floquet exponent candidates based on their real parts, before applying the imaginary part criterion to the most highly populated groups. However, in this real-partbased method, it is unclear how to proceed if two or more true Floquet exponents have the same real parts, and thus not enough groups are available.
An alternative criterion selects those candidates whose eigenvectors are most symmetric [27,28], as they should correspond most to the middle rows of the matrix H. This symmetry is computed based on a weighted mean. Even though there currently is no formal convergence proof for this symmetry-based sorting method, some numerical results indicate faster convergence than with the aforementioned eigenvalue criterion [12], while other results do not support this claim [31].
For all these criteria, there are currently no methods to efficiently and accurately compute only those eigenvalues which fit the criterion. Rather, all eigenpairs of the large matrix have to be computed first and then most of them are discarded. As the cost of solving an eigenvalue problem of a N × N matrix is of the order O(N 3 ), the computational cost of the approach is usually dominated by determining the eigendecomposition of a large matrix [29]. In addition, it is well known that the accuracy of the computed eigenvalues is not too high if all eigenvalues of a sparse matrix are sought. Sparsity of H cannot be reasonably exploited to decrease the computational cost of solving the complete eigenproblem [41].

A Koopman dictionary for time-periodic systems
For smooth autonomous systems in the Koopman framework, it is customary to choose as basis functions the so-called Carleman basis [2,42], i.e., a finite set of monomials ψ β (x) = x β , where β ∈ N n is a multiindex and standard multi-index calculation rules (see Appendix A) apply. As time-periodic functions are considered here, we propose in this paper to include as basis functions combinations of monomial terms as well as Fourier terms of the base frequency, i.e., basis functions of the form ψ β,k := x β e ikωt , where ω = 2π T . The functions ψ β,k |k ∈ Z, β ∈ N n 0 are an orthonormal system within the initially considered vector space F w.r.t. the inner product This inner product contains the standard inner product for Fourier series, and the derivatives serve as an inner product for the monomials. The inner product properties can be readily verified. Let N z , N u ∈ N be integers which describe the assumed maximum polynomial and frequency order, respectively. There is a set B = {β j } N β j=1 collecting all multi-indices with 1 ≤ β ≤ N z . These are all multi-indices that create monomials x β of degree N z and less. By (A4), it holds that N β = N z +n n − 1. For the sake of brevity, define N = N β (2N u + 1) and collects all basis functions that are not dependent on the state, but only on time, such that the evolution of u is not a product of the system dynamics, but is known a priori. All other state-dependent basis functions are collected into the vector Here, the basis functions are ordered by monomial exponent first, and then all frequency base functions for the same monomial are grouped together in ascending order. It is notable that any other orderings are also applicable and the matrices can be transformed into each other by similarity transforms. The ordering (22) has been chosen purely for convenience in later argumentation.
If the full basis as introduced above is considered, the state x itself is included in the basis. Hence, there is a selector matrix C z ∈ R n×N containing select rows of the identity matrix with C z z (x, t) = x. There exist other options to recover x from z (x, t). For instance, monomials of higher (uneven) order can be used by considering the corresponding root. Also, the matrix C z can be allowed to be time-dependent. This second option and its implications will be investigated in more detail in Sect. 4.3. Both vectors z and u are collected into one vector As introduced in Sect. 2.2, the Koopman lift describes the evolution of the basis functions under the approximate infinitesimal Koopman generator, and it is determined by expressing the time derivative along the flow of all basis functions in the finite-dimensional basis, projecting them back again to the subspace. Using the projection defined by the inner product (19), this giveṡ where r ∈ FN is the remainder that is orthogonal to the finite-dimensional subspace and will be projected out. Substituting (23) into the definition (9) for the Koopman lift and separating the two vectors of basis functions yieldṡ where A ∈ C N ×N and B ∈ C N ×(2N u +1) are constant coefficient matrices. The lower rows of the large matrix in (24b) determine the dynamics of u. However, since u is purely time-dependent and its time evolution is known a priori, the lower rows of (24b) are superfluous and the original nonlinear system is approximated by the LTI systeṁ for the very specific input u as in (21). The lifted state vector z(t) is an approximation to the state-dependent basis functions z (x(t), t) evaluated along the flow of the original system. This is only an approximation, however, due to the projection onto the finite-dimensional function space. In particular, for t ≥ 0, the lifted state z will cease to adhere to the constraints posed by z , meaning that there may not exist any vectorx ∈ R n which fulfills Many results of this work tie in to classical results that are obtained from Jacobian linearization of the state space. Therefore, the maximum monomial order in the basis will often be restricted to N z = 1. In this case, the vector of basis functions will be denoted by z,lin of size n(2N u + 1) and for the sake of brevity, these functions will be called linear basis functions, even though the influence of time is still decidedly nonlinear. To make the ordering unique, these linear basis functions are ordered in an ascending order, first by state and then by frequency such that the resulting vector reads as 3.1 Koopman lift on the perturbed system While a nonlinear dynamical system may have an arbitrary number of attractors of any kind (equilibrium, periodic, quasi-periodic, chaotic) located anywhere in the state space, the finite-dimensional linear Koopman lift (25a) has a very limited range of possible attractors, i.e., at most one single globally attractive solution, which is an equilibrium point. This means that the system (25a) will most likely not be able to describe the system dynamics of a nonlinear time-periodic system globally. For time-periodic systems, periodic solutions are expected. In this section, the Koopman lift of the perturbed system around such a periodic solution will be regarded to assess its local behavior. As in the standard HBM, a periodic solution is approximated by its Fourier expansion (10) up to order N HBM with unknown Fourier coefficients. The perturbed system is then given by as in (12). This allows the Koopman lift to be performed on the perturbed system, i.e., on functions of the state y evaluated along the flowf. Now, the origin of the lifted state does correspond to the origin equilibrium of the perturbed dynamics, or a periodic solution of the original dynamics. This also means that the Koopman lift matrices A and B now depend on the Fourier coefficients p of the periodic solution around which we are linearizing. The following sections shed light on the properties of the matrices A(p) and B(p).

Koopman-based harmonic balance and Hill equations
One of our key findings is that the B(p) matrix of the Koopman lift contains information about the parameterization of the periodic solution that is also encoded in the HBM. This is made more precise in the theorems given in this section. Let x p = k p k e ikωt be a real-valued periodic solution candidate of the time-periodic system (6). The harmonic residual of this candidate is given by Since x p as well as f are periodic, the residual is periodic as well, and therefore it has a Fourier series The HBM of order M returns solution candidates for which r k = 0 for all |k| ≤ M. Since r is real-valued, it also holds that r −k =r k , where r k = [r k,1 , . . . , r k,n ] T is a vector with n complex-valued entries. This definition of the residual allows to concisely relate the HBM to the Koopman lift.
Theorem 1 Letż = A(p)z + B(p)u be the lifted dynamics of frequency order N HBM of system (6) around an unknown periodic ansatz of the form (10). The N HBM -th order HBM equations (11), i.e., r k = 0, |k| ≤ N HBM , are given by C z B(p) = 0, where C z is the constant selection matrix that fulfills y = C z z (y, t) for all t.
The formal proof of this theorem is given in Appendix B.1. From an intuitive point of view, this condition is not surprising. For a periodic solution, the lifted dynamics has an equilibrium at zero, meaning that if y = 0 it should also hold thatẏ = 0. Since the lifted dynamics approximates z(t) ≈ z (y, t) and this approximation holds identically at the initial point (25b), it should hold thatẏ ≈ C zż is zero if y = 0. This equation can be evaluated for an arbitrary initial point on the periodic ansatz via a time-shift. With y(0) = 0 and z(0) = z(y(0),0) = 0, it turns out that the only remaining summand in the approximate ydynamics is C z B(p).
If only basis functions that are linear in the state (N z = 1) are considered, the above result still holds.
Moreover, we can state Theorem 2 about all entries of B and not just specific rows. (6) with linear basis functions z,lin of frequency order N u that are sorted as in (22), evaluated for the perturbed system around an unknown periodic ansatz of the form (10) up to frequency order at least N HBM = 2N u . Then, the matrix B(p) ∈ C n(2N u +1)×(2N u +1) consists of n stacked Toeplitz matrices. The l-th Toeplitz matrix B l contains as entries (ignoring duplicates) precisely the 4N u + 1 residuals r k,l (p), |k| ≤ 2N u that follow from the HBM w.r.t the l-th state.

Theorem 2 Letż = A(p)z + B(p)u be the lifted dynamics of system
If B(p) = 0, then all these residuals of the HBM vanish. Conversely, if p solves the HBM equations r k (p) = 0, |k| ≤ 2N u , then it holds that B(p) = 0.
The formal proof of this theorem is given in Appendix B.1.
In addition to the B matrix, the A matrix of the Koopman lift also holds frequency information about stability of the periodic solution. This is summarized in the following theorem. (6) with linear basis functions z,lin of frequency order N u that are ordered as in (22). Then the Hill matrix H, truncated to frequency order N u , for the periodic solution parameterized by p results from the matrix A(p) by the similarity transform H = UA(p)U T , where U is an orthogonal permutation matrix that satisfies U z,lin = (y T e i N u ωt , . . . , y T e −i N u ωt ) T .

Theorem 3 Letż = A(p)z + B(p)u be the lifted dynamics around a periodic solution of system
The formal proof of this theorem, again based on explicit evaluation of the inner product in the Koopman lift, is given in Appendix B.2.
With the three above theorems, qualitative insight about the accuracy of the presented Koopman lift can be gained. Locally (in the vicinity of a periodic solution), the lifted system contains the same dynamical information that is encapsulated in the Hill matrix of the same frequency order. Convergence results about the HBM [22] and the Hill method [26] can thus be related to the accuracy of the Koopman lift. Within the applied Koopman community, such a link is quite unusual. For many applications, convergence results for the finite-dimensional Koopman lift do not exist at all. The convergence results that do exist (e.g., [43,44]) state that, under special conditions, a given error tolerance can be reached using a large number of specific basis functions, without giving an explicit upper bound. Hence, for time-periodic systems of the form (6), the proposed class of basis functions is a favorable choice due to its added connection with the Hill matrix, which indicates that the Koopman lift retains valuable stability information.

Sorting-free stability method
The Koopman lift with Theorems 2 and 3 gives a dynamical systems interpretation for the Hill matrix, allowing for a novel stability method based on the Hill matrix. This is demonstrated in the following section.

Approximating the monodromy matrix
Consider the Koopman lift as in Theorems 2 and 3, i.e., consider the linear basis z,lin evaluated for the perturbed system around a periodic solution which is determined up to a frequency order of 2N u . From Theorem 2, we know that B = 0 and from Theorem 3, that A = U T HU. The linear dynamical systemż = Az+Bu resulting from the Koopman lift therefore reduces tȯ where C(t) ∈ C n×N is a possibly time-dependent projection matrix that satisfies C(t) z,lin (y, t) = y for all t ∈ [0, T ). There is a n N u -parameter family of choices for C if it is allowed to be time-dependent, which will be investigated further in Sect. 4.3. However, the naive choice would be to pick the entries in (26) that correspond to frequency zero. In this case, the matrix C is constant and given by where I n×n is the n × n identity matrix, ⊗ denotes the Kronecker product and the second matrix in (31) is a row vector ∈ R 2N u +1 with zeros everywhere except for the middle column.
Since for t = 0 all exponential terms are 1 and vanish in (26), the vector z,lin (y, 0) can be expressed as a matrix product z,lin (y, 0) = Wy (32) for all y. The matrix W ∈ R N ×n consists of (repeated) rows of the identity matrix. More specifically, it can be expressed as where the second matrix in (33) is a column vector ∈ R 2N u +1 filled with ones. As the system (30a)-(30c) is a linear time-invariant (LTI) system (except for the possibly time-dependent matrix C), its closed form solution can be explicitly computed as = C(t)e At Wy(0).
As a key finding of the current paper, the matrix C(t)e At W ∈ R n×n is an approximation of the fundamental solution matrix (t), which is the matrix that satisfies y(t) = (t)y(0). In particular, for t = T , the monodromy matrix is approximated via If the Hill matrix H is already known due to other computations, e.g., as a by-product of a frequency-based continuation method such as MANLAB [27], then the similarity transform can be substituted into the matrix exponential to yield whereC,W can also be computed directly viã are also applicable. Equations (35)-(37b) show that the choice of basis function ordering, and thus the exact numerical contents of the matrix, differ only formally and can be easily transformed into each other. For the sake of clarity, only the approximation (35) will be used in the following sections, unless otherwise stated. All results can, however, be transferred analogously to the formulation (4.1).

Stability and computational effort
As (35) is an approximation of the monodromy matrix, the Floquet multipliers can be approximated by its eigenvalues. This constitutes a novel projection-based stability method based on the Hill matrix. It is labeled sorting-free because the approach does not use the same steps as standard stability methods based on the Hill matrix [28,30,31]. In all standard approaches, the complete set of eigenvalues of the Hill matrix are determined in the first step. This is a computationally intense operation which may lack accuracy [41]. The result of this first operation is N = n(2N u + 1) Floquet exponent candidates, of which only n are the best approximations to the true Floquet exponents. There are various sorting algorithms that vary in computational expense, which aim to find these n best approximations as introduced in Sect. 2.4.
In contrast to these approaches, we now propose a sorting-free method. The Hill matrix is constructed identically to the aforementioned approaches. However, once the Hill matrix is obtained, the sortingfree projection-based stability method goes a different route. The approximation to the monodromy matrix is determined first via (35). In the most general case, a matrix exponential is of the same complexity range O(N 3 ) as an eigenvalue problem [41,45], however this is not the case for the specific operations presented here. We give two reasons: 1. The matrix exponential in (35) is always multiplied by the matrix W, which is smaller and relatively sparse. The matrix exponential can thus be computed by multiple evaluations of the action of the matrix exponential on a sparse vector. For this operation, efficient scaling and squaring approaches which utilize these properties exist [46]. 2. In many applications the constructed Hill matrix will be relatively sparse. This is the case if the fre- As a second step in the sorting-free approach, it remains to solve an eigenvalue problem for the approximation of T . This matrix is only of the size n × n, i.e., in general much smaller than the size n(2N u +1)×n(2N u +1) of the Hill matrix H. The result are n Floquet multipliers and no a posteriori sorting of candidates is necessary. The reduction of candidates is performed implicitly by the projection matrices C and W. This reduction through projection is essentially different from sorting, as the resulting Floquet multipliers from the projectionbased method do generally not coincide identically with any Floquet multiplier candidates obtained as eigenvalues of the Hill matrix, transformed from Floquet exponents to Floquet multipliers via (15). Alternatively to the Hill matrix approaches, the monodromy matrix can be determined via time-integration of the variational equation (14) [29]. Afterward, it again remains to solve the eigenvalue problem on the n × n monodromy matrix. The three general approaches are visualized in Fig. 2. While the starting point for the novel method is the Hill matrix as in the standard Hill methods, the final steps of our method are instead identical to the time-integration method. However, the second indicated step, which needs the most computational effort in all three approaches, is different. The presented sorting-free projection method therefore has the advantage of being a Hill-based method, which is favorable in an HBM setting, and at the same time only requiring to compute the smaller eigenvalue problem of the monodromy matrix, similar to the time-integration-based method.

Choice of projection matrix
Until now, it has been established that C(T )e AT W approximates the monodromy matrix if C z,lin = y. However, the choice of the matrix C to achieve this is not unique. In addition to the naive choice (31), many more matrices are admissible. Recall that the vector z,lin (26) is sorted such that it contains n consecutive blocks of length 2N u + 1, each block collecting all functions of a single state. Therefore, the l-th row of C, which singles out the l-th state, consists of n − 1 blocks of zeros, each of length 2N u + 1, and one block c l (t) ∈ R 1×(2N u +1) that is possibly filled with nonzero entries, yielding a non-square block-diagonal structure with c l (t) ∈ C 1×(2N u +1) ; l = 1, . . . , n. With this notation, the requirement C(t) z,lin (y, t) = y for all t ∈ [0, T ] can be simplified to by noting that y l can be pulled out on both sides of the equation. As all entries in u are linearly independent, the time dependency of c l is uniquely given via to cancel out the time dependency in (39), such that V(t)u(t) is constant. In this expression,ĉ l ∈ C 1× (2N u +1) is a constant row vector and V(t) is given by V(t) = diag e i N u ωt , . . . , e −i N u ωt . Due to (39) and y being real, two additional conditions onĉ l are All choices forĉ l that satisfy these conditions are admissible for the projection matrix. If a set ĉ l,k N u k=1 of N u arbitrary independent complex positive-frequency coefficients is given, they can be easily extended to an admissible choiceĉ l using the conditions (41), which implies that the admissible projection matrices form a n N u -parameter family. For the sake of readability, the constraints will be left explicit in the considerations below. In Fig. 6-8 of Sect. 5.2, two choices forĉ are compared and it turns out that this choice indeed strongly influences the approximation quality. For numerical computations, handling the sparse matrix C with only few unknown coefficients is cumbersome. It is easier to collect all unknown variableŝ c l , l = 1, . . . , n into a long row vectorĉ all witĥ If the true monodromy matrix T would be known, then (35) can be viewed as a fitting problem for the unknownĉ all . The squared matrix residual for l = 1, . . . , n. In this equation, the matrix Q ∈ C n(2N u +1)×(2N u +1) is divided into n square blocks and due to the diagonal structure of C(T ), only the l-th block of Q influences the l-th row of the least-squares problem (43). A matrix C that was determined for comparison purposes using this least-squares problem under knowledge of the true monodromy matrix will be referred to as C true below. This least-squares problem has n equations for N u independent unknowns. Therefore, to enable an optimal solution, it is advisable that the frequency order N u is chosen to be at least n.
As the aim of this method is approximating the true monodromy matrix, which is unknown in applications, the least-squares problem (44) cannot be utilized in practice. As a practicable condition, the matrix C(t) can be chosen such that it optimally satisfies the variational equation (14). The residual of (14) for the approximated fundamental matrix is and a cost function can be defined via Similarly to the block-based approach for the leastsquares problem, this function with a matrix argument can be transformed into a quadratic cost function (47) with the vectorĉ all as argument. This transformation of R(t) is detailed below. By product rule, the total time derivative of the first summand in (45) yields

(C(t)Q(t)) · = C(t) (A + D(t)) Q(t) =: C(t)L(t) ,
where D = I n×n ⊗ diag(i N u ωt, . . . , −i N u ωt) (cf. (40)) is the matrix that satisfiesĊ = CD andQ = AQ follows from its definition. Below, the dependency on t in the matrices is omitted for the sake of brevity. Substitution of Q, L into (45) yields To exploit the diagonal structure (38) of C, the matrices Q, L ∈ C n(2N u +1)×(2N u +1) are segmented into stacks of column vectors of length 2N u + 1 via and Q analogously. For the first summand in (49), this gives where each of the indicated entries is a time-dependent scalar in C. Similarly, the second summand can be separated into its entries to yield The still time-dependent c l (t) is generated from the constant coefficientsĉ l via the diagonal matrix V(t) as in (40). Collecting all unknown coefficients into one large vectorĉ all , the (i, j)-th scalar entry of the matrices in (49) can be expressed by In (53a), only the i-th block is nonzero. If the Frobenius norm is used for the residual, all absolute squared values of the entries of the matrix (49) are summed, and this gives the expression Finally, noting thatĉ all is not time-dependent, it can be pulled outside the integral to yield the quadratic cost function (47) with Therefore, the search for the best choice for the projection matrix C var can be formulated as a quadratic program where the equality constraint encodes the normalization condition (41a). The condition (41b) is not explicitly stated in the quadratic program. This is because minimizers of (56) always fulfill this condition. From an arbitrary candidateĉ all , one can construct a symmetric candidateĉ sym which fulfills (41b) and for which R has the same real part and zero imaginary part, meaning that R 2 can not be larger for the symmetric candidate than for the non-symmetric one. The proof of this is sketched in Appendix B.3.
As the matrix (t) is of size n(2N u +1)×n(2N u +1), and therefore rather large, efficient numerical determination of the integral (56a) is crucial to retain a numerically efficient stability method. Because the timedependency of is introduced by terms of the form e ikωt , the integral can be reformulated into a (infinite) sum of inner products with Fourier base functions if the power series expression for the matrix exponential is used. This suggests that the integral could be computed efficiently using FFT methods. However, due to the matrix exponential there is no limit in frequency of these terms and aliasing effects must be considered.
Alternatively, the integral can be determined using numerical quadrature schemes. This is accurate, but computationally expensive. Finally, it is also an option to require the integrand to be zero only at specific time instants. This reduces the quadratic program to a linear equation system. While this approach does not minimize the quadratic program (56) in an integral sense, the resulting approximate monodromy matrix seems to be relatively accurate in application.

Numerical examples
The theoretical results of the previous sections will be illustrated in this section using some numerical examples, which will allow us to demonstrate the numerical efficiency and accuracy of the proposed projectionbased Hill method. The Mathieu equation is utilized as a simple linear time-periodic system of the form (14) to explicitly illustrate the Koopman lift and the projectionbased stability approach. To demonstrate the computational advantage for larger numerical orders, the linearization of the vertically excited n-pendulum is then considered as a generalization of the Mathieu equation with arbitrary degrees of freedom.

Mathieu equation
The Mathieu equation is an example of a Hill differential equation [47], which has become very well-known since it results from linearization of a number of applications, among them rolling of container ships [48] and a vertically excited pendulum [49]. After bringing the system into firstorder-forṁ (58) is a linear periodic time-varying homogeneous system of the form (12) with system periodT = π ω , meaning that it is suitable to explore the stability methods of Sect. 4. Necessarily, the system is also T -periodic with T = 2T = 2π ω . It is known that the only parameter combinations of (58) which admit non-trivial T -periodic solutions are located at the stability boundaries [39,50]. The stability boundaries of the Mathieu equation can therefore be identified using the HBM or the shooting method. Below, the base frequency ω (i.e., including the first subharmonics of the system) has been chosen for investigation using the frequency-based methods. This is due to two reasons: First, as the stability boundaries may admit T -periodic or T 2 -periodic solutions, this choice allows to identify all stability boundaries using one unified HBM approach without further distinction. Second, in the projection-based Hill method, this choice of base frequency allows to better showcase the influence of the projection matrix.
The equivalence of the Koopman lift matrices, the Hill matrix and the HBM are explicitly verified for the smallest possible frequency order N u = 1. First, the Hill matrix is constructed in the standard way. The Fourier decomposition of J(t) can be read directly from (58). The nonzero Fourier coefficients with base frequency ω are and all others vanish. By construction using (18), the Hill matrix of frequency order 1 reads The block-Toeplitz-like structure with additional diagonal elements is visible in (61). Furthermore, it is obvious through the parameter b that terms of harmonic 2 influence the Hill matrix, even if the frequency order was chosen to be N u = 1 < 2. Further, assume there exists a (potentially nontrivial) periodic solution x p = j p j e i jωt . Because the original dynamics (58) is already linear, this also holds for the perturbed dynamicṡ With (62) expressed explicitly in Fourier series forṁ for the first state and for the second state. The coefficients for the basis functions can be identified in (64) by inspection. For the state-dependent basis functions with N u = 1, this yields The blocks of A have been visually separated to indicate the dependence on the individual states. The selection matrix and B 2 omitted for the sake of brevity. The expected Toeplitz structure and the HBM equations up to order 2N u are clearly visible. It is also easy to see that B 1 holds the HBM equations up to order 2 for the first (time-independent) row of (58). The accuracy of a stability traverse is investigated in Fig. 3 for the Mathieu equation with ω = 1 and b = 1.21. The parameter a takes the values −0.367 and −0.3673, respectively. Between these values, a stability change from unstable to stable occurs, where the pair of Floquet multipliers meet at 1, and continue on the unit circle as a complex conjugated pair. The Floquet multipliers labeled as true in Fig. 3 were determined using the time-integration method with a high relative and absolute tolerance of 10 −14 each. In contrast, all Hill-matrix-based methods rely on a Hill matrix of relatively low frequency order N u = 4, such that the differences in performance are more visible.
From the figures, it is apparent that the projectionbased method with an optimized projection matrix C var (i.e., by solving (56)) can be able to track the true Floquet multipliers more closely than any Floquet multiplier candidates obtained directly from the eigenval- The stability regions of the Mathieu equation are often visualized in a so-called Ince-Strutt diagram [51]. Figure 4 showcases the properties of the different approaches for drawing this stability map using the Hill matrix of a relatively low frequency order N u = 4. The color indicates the absolute value of the largest Floquet multiplier (in magnitude). Dark blue colors indicate a maximum value of exactly 1, i.e., stable regions. White color indicates a maximum value larger than 1, i.e., unstable regions.
For Hill equations, it is known that the product of both Floquet multipliers must always equal to 1 [39], i.e., the largest Floquet multiplier may never admit an absolute value smaller than one. In the diagrams, red and purple regions correspond to parameter combinations where the largest identified Floquet multiplier has an absolute value smaller than one, meaning that the identified Floquet multipliers can certainly not be the true ones. In every figure, the true stability boundaries are given in green by the solution of an accurate shooting method.
If all eigenvalues of the Hill matrix are considered for stability (naive approach, Fig. 4a), there are large regions that are wrongly classified as unstable, while none of the unstable regions are wrongly classified as stable. This is expected since in the unstable case, the best approximations will always be among the considered eigenvalues, but in the stable case the additional candidates add additional possibilities of asserted instability. In contrast, the classical Hill method with sorting procedures, which is the current state-of-the-art, classifies most of the stable regions correctly. This is visualized in Fig. 4b for the imaginary-part-based criterion [23,30] and in Fig. 4c for the symmetry-based criterion [27,28]. However, in the instability tongue that is separated by non-trivial solutions of period T , red artifacts that are classified as stable are visible in both approaches. These correspond to cases where the sorting algorithm did not choose the correct Floquet multipliers (of which one would be stable, i.e., realvalued and < 1, and the other unstable, i.e., real-valued and > 1), but rather chose two instances of the stable class, both with real parts < 1. With the naive projection matrix, the novel projection-based approach preserves most of the stability regions, while not exhibiting any stable artifacts. However, at the edges of the stable regions, in particular around a = 0.5 and b = 1, slightly larger values of the magnitude are visible. The stable regions are slightly underestimated. A very similar behavior can also be observed with the optimized projection matrix C var . The example with the naive projection matrix shows that the presented method can be very accurate, even in its simplest form which is easy to implement since only the matrix equation (4.1) is needed. The Mathieu equation considered in Sect. 5.1 can result from linearization of a vertically excited mathematical pendulum, also called the Kapitza pendulum [52]. As a scalable generalization for arbitrary degrees of freedom, the linearized dynamics of a vertically excited multiple pendulum is considered. A sketch of the considered mechanical system is given in Fig. 5. The pendulum consists of n p joints, each of mass m, with viscous absolute dampingd, linked by n p rods of length l. The minimal coordinates θ = θ 1 , . . . , θ n p T are the absolute angles of the individual joints. The suspension point of the pendulum moves vertically with y 0 (t) =ŷ 0 cos(2ωt). Gravitation acts in the vertical direction.
The equations of motion for the vertically excited pendulum can be derived similar to [53]. We define the auxiliary vectors s(θ ) T := n p sin θ 1 , (n p −1) sin θ 2 , . . . , sin θ n p (68a) c(θ ) T := n p cos θ 1 , (n p −1) cos θ 2 , . . . , cos θ n p (68b) as well as matrices S(θ ), C(θ ) ∈ R n p ×n p with Dropping the arguments for the sake of brevity, the potential energy V (θ, t) and the kinetic energy T (θ ,θ , t) are given by Using the Lagrange equations of the second kind (see, e.g., [54, p. 76]), the equations of motion are given by after some algebra. Introducing the abbreviations and linearizing around the origin, the first-order linearized dynamics of the vertically excited pendulum is which, after inversion of the left matrix, is of the formẏ = J(t)y that can be analyzed using the presented methods. The total number of system states is n = 2n p . The relation between the classical Mathieu equation (58) and the linearized vertically excited multiple pendulum is visible from (72). The total number n of states as considered throughout this paper is n = 2n p .
To analyze the convergence behavior of our proposed method, the accuracy of the Floquet multipliers of the equilibrium for one specific parameter set (a, b, d) is evaluated for various truncation orders N u and using the various approaches discussed in this paper. As a basis for comparison, the "true" Floquet multipliers are determined by integrating the variational equation (14) using the MATLAB function ode45 with absolute and relative tolerances of For the two standard Hill approaches, the eigenvalues and eigenvectors of the Hill matrix H were determined using the MATLAB procedure eig. The sorting procedure based on the imaginary part as described in [23,26,30] then singles out the n eigenvalues with least imaginary part in modulus, while the symmetrybased sorting procedure singles out the n eigenvalues whose eigenvectors have lowest weighted mean according to [27]. The Floquet multipliers are then determined from the Floquet exponents (i.e., the identified eigenvalues) using (15). The corresponding errors are depicted in Fig. 6 in dashed blue and dotted red for the imaginary part criterion and the symmetry criterion, respectively. The presented novel Koopman-based method is evaluated for two choices of the projection matrix C (see Sect. 4.3). The Floquet multipliers are determined as the eigenvalues of the approximate monodromy matrix (35). The matrix exponential in (35) is evaluated directly with its action onto the matrix W using the MATLAB function exmpv [46,55]. The dash-dotted yellow line shows the accuracy of the Floquet multipliers obtained from the monodromy approximation using the naive projection choice C 0 (31), while the solid purple line indicates the accuracy obtained by a more informed projection matrix which was obtained by inspection after running the optimization (56) for a few test cases. From the figure, it is visible that all considered approaches eventually converge to the true Floquet multipliers. The final error of order 10 −12 can be attributed to inaccuracies of the numerical integration. The imaginary-part-based criterion, while being the only one with a rigorous convergence proof for n → ∞, performs highly inaccurate for smaller N u . Moreover, the choice of projection matrix significantly influences the performance of the novel projection-based approach. While the naive choice of projection matrix does converge toward the correct value, the optimization-based projection matrix exhibits a significantly better convergence rate which can compete with the symmetry-based criterion or outperforms it. The data set of Fig. 6 is again plotted in Fig. 7, but against the computation time instead of the frequency order, giving a work-precision diagram. The strength of the projection-based methods for higher frequency orders is apparent in this figure. While convergence of the imaginary-part-based method with respect to the frequency order N u is better than that of the new projection-based method for the naive projection choice C 0 , the accuracy of the projection-based method is higher than that of the imaginary-part-based method for any computation time. The reason for this is that the eigenproblem becomes more expensive than the matrix exponential for higher matrix dimensions, i.e., for higher N u . The same property can also be observed regarding the symmetry-based approach and the novel projection-based approach with optimized projection matrix C 1 . The projection-based approach is able to achieve maximum accuracy already at approx. 0.1s compared to approx. 0.25s for the symmetrybased approach. This is a ratio of 2 5 , which is considerably lower than the corresponding ratio 15 18 of frequency orders. If the system dimension increases, this efficiency difference increases as well as the size n(2N u + 1) × n(2N u + 1) of the Hill matrix is influenced by a product of n and N u . Therefore, the superior efficiency of the projection-based approaches against the standard approaches is even more prominent. Figure 8 again compares computation time against accuracy for the various approaches, but for a pendulum with 15 instead of 6 links. Since the number of states is increased, the computation time for all the considered approaches is larger than in Fig. 7. However, the increase in number of states impacts the standard Hill methods which rely on solving the large eigenvalue problem more than the projection-based method. This is because of the key feature of the new method: a moderate growth of computational cost as sparsity can be exploited in the action of the matrix exponential. In fact, in Fig. 8, the break-even point between the symmetrybased method and the projection-based method with the optimized projection matrix already occurs at N u = 6. The imaginary-part-based Hill method fails to converge within the depicted computation time interval.

Conclusions
In this paper, a relationship between the HBM, the Hill method and the Koopman lift with specific basis functions has been addressed. It has been shown that the Hill matrix is the system matrix of a linear time-invariant dynamical system of higher order which approximates the perturbed dynamics. This relationship has been used to derive a novel stability method based on projecting the Hill matrix down to original system size, instead of computing all its eigenvalues. The resulting stability criterion in its naive form is remarkably simple to implement. It only needs a matrix exponential and two projection matrices, which are the same for all systems and easy to construct (see (31), (33)). It has been shown in the examples of Sect. 5 that this simple method can compete with the state-of-the-art algorithms in terms of Floquet multiplier accuracy over computational effort.
To further improve the accuracy of the proposed method, an optimization-based approach for the choice of the projection matrix through a quadratic program has been presented in Sect. 4.3. While the determination of the matrix integral in (56) is more computationally expensive, it has been shown in Sect. 5 that this optimized projection matrix allows for even better accuracy.
There are many topics for further research expanding from the results of this work: The approach itself offers opportunities for additional development both in its theoretical as well as for its computational properties. Further, there are many possible interesting applications of the approach to technical systems, and extensions to wider classes of systems would be valuable.
With respect to theoretical properties, there is currently no rigorous convergence proof for the novel approach for any choice of projection matrix. Such a convergence proof would aid in an a priori determination which method to employ for the best results.
Considering computational properties, while the presented approach does reduce computational effort compared to the state-of-the-art approaches, some bottlenecks remain. In particular, the efficient sparsitypromoting determination of the matrix exponential with projection from both sides as well as the integral in the optimization problem (56) would benefit from more efficient computation techniques. As the latter is constructed by products of the Fourier coefficients of J in A and additional Fourier terms, there could exist an operation to obtain this integral directly from the (FFT) Fourier series of J(t).
From an application point of view, going beyond the simple and mainly academical examples presented in this work, it would also be expedient to apply the novel stability method to systems of more practical relevance. To achieve this goal, it would be beneficial to integrate the novel method within an established continuation framework. The authors believe that the MANLAB continuation framework [25,28] would be especially suitable because the Hill matrix can be constructed in this framework without much additional effort. Large systems that were already analyzed extensively using the MANLAB framework (e.g., [56]) could benefit from the stability insights of our proposed method, while simultaneously serving as practical examples of the performance of the method.
Regarding possible extensions, while some relations were established in Theorem 1 for polynomial basis functions of the Koopman lift, most sections of this work rely on basis functions that are linear in the state. Some recent results [57,58] make statements about return time and period for autonomous systems based on the Carleman linearization. As our (full) basis provides a natural extension of the Carleman linearization, these results could possibly be integrated into the presented framework.
Finally, it would be of practical value to extend the method to a wider class of systems, for example, possibly including delay differential equations (DDEs) with delayed coordinates in the basis functions.
Author contributions Both authors contributed to the conceptualization, design and methodology. The formal analysis was performed by Fabia Bayer. The original draft of the manuscript was prepared by Fabia Bayer and both authors commented on previous versions of the manuscript. Remco Leine was respon-sible for funding acquisition and supervision. Both authors read and approved the final manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL. The authors declare that no funds, grants, or other support was received during the preparation of this manuscript.
Data availability All numerical data in this work have been generated from the considered system Eqs. (58), (72), with the approaches described in this work and the cited works, using MATLAB standard methods and the MATLAB function expmv [55]. Therefore, it is possible to completely reproduce the data from the information given in this work, except for the computation times in the work-precision diagrams. These computation times have been evaluated on a Windows PC, and more detailed computation time data can be obtained from the authors upon request.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Appendix A: Multi-index calculation rules
For convenience, the standard multi-index calculation rules (see, for example, [59, p. 319]) are revisited. A multi-index is a tuple β ∈ N d 0 with nonnegative integer entries. In particular, a multi-index β = (β 1 , . . . , β d ) has the 1-norm and a factorial given by the product of the factorials of all its com- The inner product with u is just a projection onto the Fourier basis functions between ±N u , and therefore, the row of B l that corresponds to frequency j is given by ψ j,l , u = (r −N u − j,l , . . . , r N u − j,l ).
Comparing two successive rows, i.e., the j-th and the ( j + 1)-th row of B l , it is obvious from (B13) that they contain the same entries, but shifted to the right by one element. Adequately, going from row j to row j +1, the entry r −N u −( j+1),l is added on the left, while the entry r N u − j,l is pushed out on the right. The matrix B l thus has a Toeplitz structure, where all residuals between r −2N u ,l (in the bottom left entry of B l with j = N u ) and r 2N u ,l (in the top right entry with j = −N u ) occur, yielding 4N u + 1 distinct entries in total.
In particular, if the coefficients p describe the HBM solution for order 2N u , then all these residuals vanish and it holds that B(p) = 0.
The arguments of the previous proof can be reused for Theorem 1 with slightly modified assumptions. Recall the theorem: Theorem 1 Letż = A(p)z + B(p)u be the lifted dynamics of frequency order N HBM of system (6) around an unknown periodic ansatz of the form (10). The N HBM -th order HBM equations (11), i.e., r k = 0, |k| ≤ N HBM , are given by C z B(p) = 0, where C z is the constant selection matrix that fulfills y = C z z (y, t) for all t.
Proof In contrast to the previously proven Theorem 2, in this theorem the basis functions are not restricted to be linear in the state. Due to sesquilinearity of the inner product (2b), however, the selection matrix C z can be pulled into the expression for B(p) via These expressions can again be evaluated row-wise with (B13) and since j = 0 for all entries of (B14), the condition C z B(p) = 0 is equivalent to r k = 0 for all |k| ≤ N HBM , which are exactly the HBM equations.

B.2 Theorem 3
Similarly to the theorems proven previously, it is central in the proof for Theorem 3 to evaluate the inner product explicitly. The theorem is restated for convenience.
Theorem 3 Letż = A(p)z + B(p)u be the lifted dynamics around a periodic solution of system (6) with linear basis functions z,lin of frequency order N u that are ordered as in (22). Then, the Hill matrix H, truncated to frequency order N u , for the periodic solution parameterized by p results from the matrix A(p) by the similarity transform H = UA(p)U T , where U is an orthogonal permutation matrix that satisfies U z,lin = (y T e i N u ωt , . . . , y T e −i N u ωt ) T .
Proof The matrix U is derived from the identity matrix by reordering its rows, and it reorders the entries of z,lin such that the resulting vector has entries descending in frequency, where all states corresponding to the same frequency are collected together. Denote˜ z := U z,lin =: ((˜ z ) T −N u , . . . , (˜ z ) T N u ) T and separate it into 2N u +1 blocks of length n with (˜ z ) k := ye −ikωt . First, we show that the Koopman lift performed with basis˜ z identically yields H as its system matrixÃ, and afterward we demonstrate how performing the lift with z,lin instead of˜ z demands the additional similarity transform.
As in the proof for Theorem 2, we start by determining the time derivative of the basis functions. For the vector of time derivatives, by the chain rule it holds thaṫ wheref is the nonlinear perturbed dynamics (27a). Evaluating the individual summands of (B17) yields  Proof Below, vectors d ∈ C (2N u +1) which fulfill (41b) will be called symmetric. The symmetric vectors form a vector space over R (but not over C). First, the structure of the Koopman lift matrix A is revisited. Applying the similarity transform of Theorem 3 to the Hill matrix explicitly yields for the matrix A the structure Hence, the l j-th block of the matrix A contains the Fourier coefficients of the l j-th entry of the system matrix J. In particular, since J(t) admits only real values, J l j,k = J l j,−k and each block A l j fulfills a variant of the symmetry property: The −k-th row is the complex conjugate of the k-th row, with the order of elements reversed. A little algebra shows that for matrices with these properties and symmetric vectors d ∈ C (2N u +1) , the symmetry is retained though multiplication: A l j d is symmetric if d is symmetric. If a matrix P = [d l j ] is constructed block-wise from symmetric column vectors d l j n l, j=1 , then the matrix AP is again constructed from symmetric column vectors since all block entries of AP can be decoupled into sums of products A li d i j , which each retain the symmetry as described above. In particular, since the matrix W in (33) is constructed in the considered fashion, the product AW again has block-wise symmetric entries. Iterative application of this multiplicative invariance to shows that the column vectors Q l j and L l j in (50) are also symmetric. Next, the scalar product b T d ∈ C of a symmetric vector d and an arbitrary vector b of appropriate size is considered. The operator flip b reverses the order of entries of a vector b. Element-wise evaluation shows Definingĉ l,sym = 1 2 ĉ l + flipĉ l and using the above relation yieldŝ and analogously for L i j . As R is constructed from summands of this form via (49)-(52), the proposition follows.
A purely real matrix R sym will necessarily have a Frobenius norm smaller or equal than the norm of a matrix with same real part and arbitrary imaginary part. Hence, Proposition 4 allows to easily construct a minimizer that satisfies all constraints in case a different minimizer is returned in the quadratic program (56).