Existence of quasi-periodic responses in quasi-periodically forced nonlinear mechanical systems

Forced responses of mechanical systems are crucial design and performance criteria. Hence, their robust and reliable calculation is of vital importance. While numerical computation of periodic responses benefits from an extensive mathematical basis, the literature for quasi-periodically forced systems is limited. More specifically, the absence of applicable and general existence criterion for quasi-periodic orbits of nonlinear mechanical systems impedes definitive conclusions from numerical methods such as harmonic balance. In this work, we establish a general existence criterion for quasi-periodically forced vibratory systems with nonlinear stiffness terms. Our criterion does not rely on any small parameters and hence is applicable for large response and forcing amplitudes. On explicit numerical examples, we demonstrate how our existence criterion can be utilized to justify subsequent numerical computations of forced responses.


Introduction
Quasi-periodic oscillations are vibrations containing multiple incommensurate frequencies. They have been reported in physics, chemistry and engineering. For example, Ditto et al. [21] study quasi-periodic oscillations of a magneto-elastic ribbon and Hauck and Schneider [30] observe quasi-periodic oscillations in an oxidation reaction. Engineering applications include Gendelman et al. [26] who utilize the emergence of quasi-periodic oscillations for vibration absorption in an electric circuit. Moreover, Coudeyras et al. [16] observe quasi-periodic vibrations impacting break squeal and Kim and Noah [37] examine quasiperiodic responses of a nonlinear Jeffcott rotor with bearing clearance.
It seems naturally to expect quasi-periodic vibrations in mechanical systems when the applied external loads are quasi-periodic. This expectation, however, is generally not true even for the periodic case, as we have investigated in depth for mechanical systems [8]. Consequently, the question arises when we can rigorously expect a quasi-periodic steady-state response of a forced nonlinear mechanical system.
Numerical time integration seems to be a natural candidate to compute the response of nonlinear mechanical systems. Due to the low damping in modern engineering structures, however, long integration times are required to observe a steady-state response. Moreover, quasi-periodic orbits lack a fixed period after which the steady-state response repeats itself. Thus, additional signal processing tools, such as Fourier transformations are required to even recognize a quasi-periodic steady state. For example Kreider and Nayfeh [40] observed "unexplained sidebands" in the measured forced response of a clamped-clamped beam. Later, Emam et al. [22] explained these sidebands with the occurrence of quasi-periodic motions. Moreover, Held and Jefferies [32], Cumming and Linsay [17] and He et al. [31] use Poincare sections (cf. Guckenheimer and Holmes [28]) to recognize quasiperiodicity. A periodic motion occurs as single point inside a Poincare section, whereas a quasi-periodic orbit with two incommensurate frequencies leaves a densely filled orbit inside the trap. Orbits with more incommensurate frequencies generally generate more complicated patterns inside a Poincare section (cf. Cumming and Linsay [17] Fig. 2b for an example). Other techniques have been used by, e.g., Walden et al. [60] who assess the number of incommensurate basis frequencies to decompose their measurements of a Rayleigh-Bérnard convection. Moreover, Laskar [42] presents a mathematically basis to decompose a signal into its quasi-periodic components including a convergence proof. Motivated by the computational success of the harmonic balance procedure in computing periodic orbits (e.g., Mickens [46]), Chua and Ushida [14] extended the harmonic balance method to the quasiperiodic case. In this procedure, an assumed quasiperiodic solution is expanded in a multi-frequency Fourier series and the governing ordinary differential equations of the nonlinear system are projected onto a finite set of Fourier modes. The arising nonlinear algebraic system is then solved with, e.g., the Newton-Raphson iteration. Recently, the harmonic balance procedure has been generalized by Schilder et al. [57], who also proposed a finite difference approximation. Their version is also implemented in the continuation package coco (cf. Dankowicz and Schilder [19]). The error due to the truncation of the generally infinite number of Fourier modes of quasi-periodic solutions, however, is unknown and the convergence of this method is a priori unclear. Furthermore, if the existence of a quasiperiodic is not shown, the expansion of the solution in a multi-frequency Fourier series is unjustified. Thus, the harmonic balance method is generally unreliable in those cases. If the mechanical system is close to a solvable limit, perturbation methods can rigorously compute quasi-periodic orbits. Examples of general perturbation methods are the method of averaging (cf. Sanders et al. [55]), normal forms (cf. Murdock [47]) or the method of multiple time scales (cf. Nayfeh [48]). Moreover, Samoilenko [54] establishes a perturbation approach of quasi-periodic tori relying on a Green's function. Similarly, the KAM-Theorem (cf. Guckenheimer and Holmes [28]) can guarantee the existence of quasi-periodic orbits in Hamiltonian systems under suitable conditions. The validity of all these methods, however, is restricted to a small neighborhood of the solvable, i.e., unphysical, limit. Whether the actual physical system is contained within this neighborhood is generally uncertain.
Fixed point theorems such as Brouwer's or Schauder's fixed point theorem (cf. Bobylev et al. [5] or Precup [51]) are powerful tools to prove the existence of solutions in the absence of small parameters. While the existence of periodic solutions is extensively treated in the literature (cf. Farkas [23], Rouche and Mawhin [52] or Bobylev [5]) the quasi-periodic or almost-periodic 1 case is less prominent (cf. Fink [24] and Corduneanu [15]).
Fink [24] guarantees the existence of an almostperiodic solution for one-dimensional first-order systems, if a monotonicity condition is met. Still relying on a monotonicity condition, Dafermos [18] extends the scalar results to abstract evolution equations in Hilbert spaces. His results can guarantee the existence of quasi-periodic responses in a fairly general setting for multi-dimensional first-order systems. The monotonicity condition, however, is not satisfied for simple mechanical systems such as the linear damped oscillator as we examine in Appendix A. Hence, the results of Dafermos [18] are inapplicable in the structural dynamics context. Berger and Chen [3] give conditions for the existence of almost-and quasi-periodic solutions for second-order systems. Their results have been extended by, e.g., Carminati [11], Bolt et al. [4] and Mahwin [45], who also cover the case of bounded and almost-periodic solutions. Moreover, their theorems can guarantee the 1 Almost-periodic functions are a generalization of quasiperiodic functions. Intuitively speaking, they are obtained by allowing for an infinity number of incommensurate vibration frequencies. For a more elaborate and rigorous introduction, we recommend the books of Fink [24], Krasnosels'kii et al. [39] and Corduneanu [15]. We note that all quasi-periodic functions are almost-periodic. uniqueness of the bounded solution. These results, however, require again a monotonicity condition on the geometric nonlinearities (i.e., position dependent terms). In structural engineering this condition requires the stiffness to be negative definite (cf. Appendix A). This condition, however, is generally not met for oscillatory systems.
Alonso and Ortega [1], notably, derive an existence criterion for bounded solutions of second-order equations, which can be directly applied to oscillatory systems having a positive definite stiffness matrix. Their result guarantees the existence of a unique, asymptotically stable and bounded solution. The condition of their theorem, however, requires the Hessian of the nonlinear stiffness terms to the globally positive definite. This requirement excludes softening nonlinearities which have been observed experimentally by, e.g., Cho et al. [13] in a micro-mechanical system and by Yuan [61] in energy harvesting applications. Furthermore, the theorem by Alonso and Ortega [1] suffers from a restrictive linear global upper bound on the nonlinear terms.
In summary, the existence criteria for quasi-periodic solutions from the literature are inapplicable even for simple nonlinear oscillators. For the periodic case, however, various results can establish the existence of periodic solutions even for high amplitude oscillations. In our previous work [8], we identified a strengthened version of a theorem by Rouche and Mawhin [52] to be the most relevant for periodically driven mechanical systems. Here, we extend our result to the quasiperiodic case.

Set-up
We consider the general N -dimensional mechanical system of the form where the mass matrix M ∈ R N ×N is symmetric and positive definite and the damping matrix C ∈ R N ×N is symmetric. The vector S(q) collects all position dependent forces, such a linear and nonlinear stiffness forces or non-potential forces (cf. Bolotin [6] for examples). The external forcing f(t) is assumed to be quasiperiodic according to the following definition where the time evolution of the angles φ is generated by K incommensurate frequencieṡ A forced response of system (1) is most commonly computed by evaluating Eq. (1) along an assumed quasi-periodic solution q * ( t) to E. (1). Projecting q * ( t) on the torus T K Eq. (1) yields a partial differential equation on the bounded domain T K , which can be discretized by, e.g., Fourier modes to yield the harmonic balance method (cf. Chua and Ushida [14]). Alternatively, Schilder et al. [57] propose a finite difference approximation. Since these methods evaluate Eq. (1) along a quasi-periodic orbit, they are only justifiable if a quasi-periodic solution to Eq. (1) actually exists. Applying the harmonic balance to a system without clarifying the existence of a steady-state response, can lead to erroneous conclusions as we have demonstrated on an explicit, periodically driven mechanical system [8]. To aid the existing numerical tools, it is thus of fundamental importance to guarantee the existence of quasi-periodic solutions to Eq. (1). For periodically excited nonlinear mechanical systems numerous existence criteria are readily available (e.g., Farkas [23], Rouche and Mawhin [52] or Bobylev [5]). For the quasi-periodic case, however, the literature lacks applicable existence results for realistic nonlinear mechanical systems as previously detailed (see also Appendix A). Due to their resemblance, one could intuitively expect that the quasi-periodic case does not significantly differ from the periodic case. Unfortunately, this expectation is false and existence criteria from the periodic case do not straightforwardly extend to the quasi-periodic case. To demonstrate a delicate differences between periodic and quasi-periodic oscillations, we consider the linear oscillator For the periodic case (K = 1 in Eq. (4)), a periodic response exists if the forcing is non-resonant, i.e., n = ω 0 holds for all integers n (cf. Burd [10] Remark 2.1). Extending these results to the more general quasiperiodic case (K > 1) one would intuitively expect that as long as the frequencies of the quasi-periodic forcing , κ are not equal to the natural frequency ω 0 (i.e., | , κ | = ω 0 ), then a quasi-periodic solution to system (4) exists. This expectation, however, is false. Even worse, for any incommensurate frequency basis there exists a quasi-periodic forcing such that the linear oscillator has no quasi-periodic solution (4). More specifically, in Appendix B we show the following: Fact 2.1 For any incommensurate frequency basis ∈ R K with at least two frequencies (K ≥ 2), there exists a continuous and bounded forcing f (φ) such that system (4) has no quasi-periodic solution.
Proof We rely on the occurrence of small denominators, which arise in system (4) for any incommensurate frequency basis. We detail the derivations leading to the above fact in Appendix B.
In summary, applying a quasi-periodic excitation to the undamped linear oscillator (4) does not guarantee a quasi-periodic response (cf. Fact 2.1), while for almost all periodic forcings (i.e., except forcing at exact resonance) a periodic steady-state response arises. This observation demonstrates that results from the periodic system do not always straightforwardly extend to the quasi-periodic case.

Existence of quasi-periodic solutions
We prove the following general result: holds. (C2) The stiffness terms S(q) derive from a potential, i.e., there exists a continuously differentiable scalar function V (q) such that holds.
(C3) The stiffness terms S(q) are continuous and globally Lipschitz, i.e., there exists a constant L > 0 such that holds. (C4) The inner product of the stiffness terms S(q) with the coordinates q grows sufficient far away from the origin. Specifically, there exists a radius r and a minimal growths α > 0, such that holds.
Proof We base our proof of Theorem 3.1 on a fixed point theorem by Schäfer [56] (also see Smart [59]). We establish a homotopy between system (1) and the equation Mq + Cq ± Mq = 0. Conditions (C1)-(C4) guarantee that the solutions remains bounded for all homotopy parameters. We detail our derivations in Appendix C and give a concise proof referring to relevant preliminary results and definitions in Sect. C.5.
Remark 3.1 Since periodic orbits are a special case of quasi-periodic solutions (i.e., K = 1 in Definition 2.1), one can also establish the existence of periodic solutions via Theorem 3.1. Comparing with our previous result [8] conditions (C1) and (C2) are identical. In the quasi-periodic case, however, we require inner product of the coordinates q and the stiffness terms S(q) to grow sufficient far from the origin (cf. condition (C4)), while in the periodic case it is sufficient that the stiffness terms are element-wise above or below the mean forcing value. Furthermore, we require the existence of a global Lipschitz constant in condition (C3). These deviations arise from the different fixed point theorems employed in the proofs of both results.

Remark 3.2
We emphasize that our Theorem 3.1 guarantees the existence of a quasi-periodic solution to Eq. (1) with the frequency base arising from the external excitation. This enables to approximate the steady-state solutions by Fourier series with known frequencies basis . In practice this is the most common approach to compute quasi-periodic orbits (e.g., Schilder et al. [57] or Chua and Ushida [14]). Competing results in a similar setting of our Theorem 3.1, e.g., by Alonso and Ortega [1] Lemma 3.2, can only guarantee the existence of a bounded solution. Hence, an approximation by Fourier series is generally not justified by such a result.
Generally, results guaranteeing the existence of a bounded solution, such as Alonso and Ortega's Lemma, can be strengthened to establish that an existing bounded solution is quasi-periodic with the same frequency base (or module) as the external loads. Thereby, Fourier series approximations with the frequency basis are again justified. In practice, however, the necessary arguments rely on very restrictive assumptions. For example, if the bounded solution is unique and stable, then general results on module containment (cf. Fink [24]) guarantee that the bounded solution is quasi-periodic with the same frequencies as the forcing. For nonlinear mechanical systems of the form (1), however, already in the periodic case multiple co-existing steady-state solutions and unstable periodic orbits are frequently observed (e.g., Nayfeh and Mook [49]).
In the case of multiple coexisting solutions, it is non-trivial whether these are actually even almostperiodic. For such a result one would need to show that only finitely many, semi-separated solutions exist (cf. Fink [24]). However, the periodically forced damped Duffing oscillator, a prominent example of system (1), can have infinite number of periodic orbits if a homoclinic tangle occurs (cf. Guckenheimer and Holmes [28]). Even if one achieves to establish that only finitely many semi-separated solutions exist, then the frequency basis of these almost-periodic solutions is generally not the same as the module of the timevarying terms (see Fink [24] for an explicit example). Thereby, an expansion of the solution in a Fourier series with frequency basis is generally unjustified.

Remark 3.3
If the stiffness terms in system (1) are polynomials such as in the classical Duffing oscillator, then condition (C3) is not satisfied. However, the Lipschitz constant L in Eq. (7) can be chosen arbitrarily large, which allows to conclude that condition (C4) is generally uncritical for realistic structures. Since any physical system will not allow for arbitrarily large displacements and velocities, equations of the form (1) modeling a realistic system are generally only valid in a bounded domain U ⊂ R N . Outside this domain of validity, Eq. (1) does not represent the reality. Thus, it is justified to modify Eq. (1) outside U such that the stiffness terms are globally Lipschitz. Such a truncation can always be carried out if S(q) is continuous. Then, S(q) is uniformly continuous on the compact domain U . Hence, S(q) is locally Lipschitz and there exists a Lipschitz constant L such that condition (C3) holds for all q 1 and q 2 in U .
We explicitly construct an example of such an envisioned truncation. We select a ball of radius r containing U (U ⊂ B r ⊂ R N ) and denote the point on the boundary of B r with minimal distance to q by q * (q). With this notation, we define the nonlinearities S(q * (q)), |q| ≥ r, where K * is a positive or negative definite matrix. We note that S 2 satisfies condition (C3) and hence if S is continuous, then the system is Lipschitz continuous and satisfies condition (C3). Thus, if the remaining conditions (C1), (C2) and (C4) of our Theorem 3.1 are met, a quasi-periodic response exists.
Remark 3. 4 We have restricted our Theorem 3.1 to time independent nonlinearities, which are customary in structural dynamics. However, our results can be extended to include time dependent nonlinearities S(q, t), if the time dependence is quasi-periodic. Then, if conditions (C2)-(C4) are satisfied for all times, then the assertion of Theorem 3.1 still holds.

Remark 3.5
The existence of a quasi-periodic orbit, e.g., established by our Theorem 3.1 does not rule out the existence of other invariant objects such as periodic orbits, fixed points or chaotic attractors. Indeed, the coexistence of multiple stable or unstable attractors is a fundamental characteristic of nonlinear dynamical systems (e.g., Guckenheimer and Holmes [28]). Moreover, a general connection between quasi-periodic motion and chaos has been established by Ruelle and Takens [53] and Newhouse et al. [50]. A quasi-periodic route to chaos has been experimentally observed by, e.g., Held and Jeffries [32], Cumming and Linsay [17] and He et al. [31]. Since Theorem 3.1 can guarantee the existence of quasi-periodic motions, it can be applied in numerical investigation of the quasi-periodic route to chaos in the nonlinear mechanical system (1).

Remark 3.6
From the existence Theorem 3.1, we can conclude, that one can observe a quasi-periodic motion if the dynamical system (1) is initialized on one of the quasi-periodic orbits. Whether other initial conditions are attracted to this quasi-periodic orbit or whether it is stable with a basin of attraction cannot be inferred from our Theorem 3.1. Moreover, to the best of our knowledge, a robust and universal tool to conclude about the stability of a quasi-periodic orbit is not available. Although only stable orbits are observable, it is also of interest to compute unstable orbits. For example, Denis et al. [20] utilize unstable periodic orbits for system identification.
Condition (C4) specifies the sign of the inner product of the nonlinearity S(q) and the coordinate q. Geometrically this implies that far away from the origin the vector S(q) points either toward the origin (q S(q) < 0) or away from the origin (q S(q) > 0). These considerations play a major role in the construction of guiding functions and also index theory (cf. Krasonosel'skii et al. [39], Krasnosel'skii and Zabreȋko [38] or Bobylev et al. [5]). The application reported, however, focus on either periodic or bounded solutions (see also Remark 3.2).
If the nonlinearities are continuously differentiable, then condition (C4) of Theorem 3.1 can be replaced by a more intuitive condition. holds.
Proof We proof the above theorem in Appendix D.
Remark 3. 7 We note that condition (C4*) requires differentiable nonlinearities, which is not required by condition (C4). Condition (C4*) is, however, generally easier to verify, since positive or negative definiteness can be checked through direct eigenvalue computation, leading minor criteria (or Sylvester's criterion) or Cholesky decomposition (cf. Horn and Johnson [33]).
On our way to Theorem 3.1, we have established an analogous result for non-smooth but L 2 -integrable external forcing, i.e., Equation (12) is satisfied by, e.g., saw-tooth or square wave functions. For such non-smooth forcing functions, we have the following result: Proof Our proof of the above theorem is also based on a fixed point theorem by Schäfer [56] (also see Smart [59]). We detail our derivations in Appendix C and give an overview of the proof in Sect. C.4.
Moreover, basic results on the smoothness of solutions to ordinary differential equation (e.g., Hale [29] or Arnol'd [2]) for finite time guarantee the smoothness of the quasi-periodic solution Proof Assume the contrary, i.e., the quasi-periodic solution q * of system (1) is non-smooth at some time instance t * . This statement is a direct contradiction to the existence of s-times continuous differentiable solutions for finite times proven by, e.g., Hale [29] or Arnol'd [2]. Hence, the quasi-periodic solutions are smooth as claimed in corollary (3.1).

Numerical examples
We demonstrate the applicability of our Theorem 3.1 on two oscillatory systems. First, we investigate a chain of oscillators proposed by Breunung and Haller [7]. This system is a natural extension of the classical Duffing oscillator or the Shaw-Pierre example [58] to higher dimensions. Subsequently, we apply Theorem 3.1 to a discretized slender beam. After concluding about the existence of a quasi-periodic steady-state response, we calculate it with the automated package proposed by Jain et al. [35].

Chain oscillator
In the following, we utilize our Theorem 3.1 to compute a steady-state response of a system of N masses depicted in Fig. 1. Two adjacent masses are connected with spring elements and linear dampers. The first mass is suspended to the wall, while the last mass is suspended to the moving base u(t). Additional to the base point excitation, the external loads f j (t) are applied to the jth mass. For a first assessment, we assume purely external loads and select forcing frequencies in resonance with the first two eigenfrequencies. Subsequently, we show that our Theorem 3.1 also allows for complex loading scenarios, e.g., excitations arising from the natural roughness of a road surface.

External loading
For the case u(t) = 0, the equation of motion for the jth mass is given by where the coordinates q 0 and q N +1 are set to zero. The example (13) includes the classical Duffing oscillator in the case N = 1 and S 1 (q 1 ) = ω 2 0 q + κq 3 . Higherdimensional variants of this system have been investigated by, e.g., Shaw and Pierre [58], Breunung and Haller [7] or Jain et al. [35]. In our previous work [8], we showed that if all the stiffness terms were hardening (∂ S j ( p)/∂ p > 0) then periodic forcing leads to a periodic response of system (13). For the quasi-periodic case, we have a similar statement holds, then system (13) has a quasi-periodic steadystate response.
Proof We verify that system (13) satisfies the conditions of Theorem 3.1 in Appendix 5.
Fact 4.1 extends the result of our previous work [8] to the case of quasi-periodic forcing and also includes softening nonlinearities. As an example, we consider a chain of five masses with the parameters m j = 1, c j = 0.02 and softening springs of the form S j (d) = d − 0.5d 3 . For this configuration, Fact 4.1 guarantees the existence of a quasi-periodic steady-state response.
Selecting the quasi-periodic forcing f 1 = 0.006 (sin( 1 t)+sin( 2 t)) and f j = 0 for j = 2, 3, . . . , N , we use the automated package of Jain et al. [35] to compute a quasi-periodic steady state. In this algorithm, the steady-state response is approximated with a finite number of Fourier coefficients. The arising nonlinear algebraic equations are solved either with the Picard or the Newton-Raphson iteration, where the response of the corresponding linear system serves as initial guess. While the existence criterion of Jain et al. [35] fails for forcing in resonance with an eigenfrequency, our Theorem 3.1 guarantees the existence of a quasi-periodic response for any forcing amplitude and frequency. We chose the forcing frequency 1 to be close to the first eigenfrequency ω 1 ≈ 0.52 and 2 in the vicinity of the second eigenfrequency ω 2 ≈ 1. In Fig. 2, we depict the L 2 -norm of the quasi-periodic steady-state response of the chain oscillator (13).

Base excitation through road roughness
Many engineering systems are subject to complex or even random excitations. For example, various norms

Fig. 2
Quasi-periodic steady-state responses of system (13) with five masses and parameters m j = 1, c j = 0.02 and S j (d) = d −0.5d 3 to the external loads f 1 = 0.006(sin( 1 t)+sin( 2 t)) and f j = 0. We select the forcing frequencies 1 and 2 in resonance with the first two natural frequencies ω 1 ≈ 0.52 and ω 2 ≈ 1 (e.g., MIL-STD-810G, ISO 4866 or ISO 1032) require the assessment of vibration amplitudes to loads with a specified power spectral density. Often it is assumed that these complex loading scenarios can be accurately modeled with a finite set of Fourier modes (cf. Loprencipe and Zoccali [43]). In such a setting the arising external loads are naturally quasi-periodic. In the following, we use our Theorem 3.1 to reliably compute an arising steady-state response to a complex excitation generate by a road surface. We assume that the base point u(t) (cf. Fig. 1) moves with a constant speed v along a rough street. We envision, e.g., parts of a vehicle. Road profiles are often classified by their spatial spectral density according to the norm ISO 8608 [34] (cf. Loprencipe and Zoccali [43]). Furthermore, the smoothened power spectral densities given in [34] are commonly used for dynamical analysis of vehicles. We follow this approach to generate a road roughness profile according to ISO 8608. The smoothened spatial power spectral density from the standard [34] is given bŷ where is the frequency in rad/m and the amplitudeŝ u( ) are given in meters. The normalization constant 0 is 1 rad/m and the constant C is quantifies the road quality. For our analysis, we assume that the road belongs to the class with the second highest quality and set C to be 64 m 3 . Following Loprencipe and Zoccali [43], we create an artificial road profile by selecting a finite number of frequencies k whereby the amplitudes are given by Eq. (15). Then, the road profile is given by where s is the distance traveled and φ k denote phases. We select the frequencies k by drawing random samples from the specified frequency interval (15). To emphasize the presence of low frequency terms, we assume an underlying beta distribution (e.g., Jonson et al. [36]) with the parameters α = 1 and β = 5. Moreover, we sample the angles φ k from a uniform distribution. We depict such a road profile with ten frequencies (K = 10 in Eq. (16)) in Fig. 3a.
We consider the chain of oscillators depicted in Fig. 1 with five masses. Assuming springs with linear and cubic coefficients the equations of motion are where we selected a linear force displacement relationship for the last spring S N +1 with stiffness k r and v denotes the constant speed (cf. Fig. 1). The forcing f r (t) arising through the base point excitation is analytic and quasi-periodic for any finite number of frequencies K . The frequency spectrum of the forcing arising from the road profile shown in Fig. 3a is depicted in Fig. 3b. Selecting the positive mass  Randomly generated road profile generated according to the norm ISO 8608 [34] and frequency spectrum of the forcing arising through base excitation with then base frequencies (K = 10). We consider system (17) Figure 3b shows that in the frequency band between 10 and 20 rad/s three eigenfrequencies and four forcing frequencies are in an immediate vicinity. Additionally, at 6 rad/s one of the forcing frequencies is in resonance with the lowest eigenfrequency. Due to these multiple resonances, we expect considerable vibration amplitudes.
We compute the arising steady-state response via the harmonic balance procedure with ten base frequencies, whereby the arising nonlinear algebraic equations are solved the Newton-Raphson iteration. As initial guess, we select the response of the linearization of system (17) at the origin. In Fig. 4, we depict the arising quasi-periodic steady-state response for the first one hundred meters. Additionally, we include the linear solution and a time series obtained by numerical time integration with matlab's ode45 algorithm serving as benchmark solution. To this end, we consider the system to be at rest one thousand seconds before reaching s = 0, i.e., we select the initial condition q(t = −10 3 ) =q(t = −10 3 ) = 0. Figure 4 reveals that the harmonic balance solution matches very well with the numeric time integration, while the linear response differs. The maximal displacement of q 3 along the arising steady-state response is 97.2 mm, while a linear analysis leads to a maximal displacement of 83.6 mm. Hence, linear analysis underestimates the response amplitude by more than 15 percent. Thus, it is essential to perform a nonlinear analysis on our example. Therein, Theorem 3.1 justifies the expansion the coordinates of system (17) in a Fourier series.

A slender beam
After applying our existence criterion 3.1 to the springmass systems (13), we advance to a beam serving as more realistic engineering structure. Beam elements are frequently used to model complex engineering structures and are readily implemented in common computational packages such as Abaqus, Autodesk or ANSYS. These computational packages discretize the model's governing partial differential equations in space. Through the discretization high-dimensional ordinary differential equations arise. In the following, we apply our existence criterion 3.1 to such a high-dimensional oscillatory system. Subsequently, it justifies numerical computations of quasi-periodic responses.
We consider the clamped-clamped beam of length L depicted in Fig. 5. The material properties and geometry such as cross-sectional area A and moment of inertia I are assumed to be constant. Nayfeh and Mook [49] derive a partial differential equation governing the vertical vibrations of slender beams. They assume that the Euler-Bernoulli hypothesis is met (i.e., plane cross sections initially perpendicular to the neutral axis remain and use a third-order approximation of the Green-Lagrange strain tensor (van Kármán strains). With the additional assumption of a small radius of gyration r , i.e., r = √ I /A 1, Nayfeh and Mook state the nondimensional partial differential equation where w(x, t) denotes the vertical displacement as a function of the non-dimensional horizontal position x and the non-dimensional time t (cf. Fig. 5). The subscripts of w in Eq. (18) refer to partial differentiation with respect to time t or space x and the constant μ is a damping parameter. We discretize the partial differential equation (18) in the spatial coordinate following a finite difference scheme with N equally space nodes x j ( j = 1, . . . , N ). The parameter h denotes the distance between neighboring nodes. We use a forward scheme for to odd spa-tial derivative w x and a central scheme for the even spatial derivatives which yields To approximate the spatial derivative at the node x N = L, we use a backward-difference scheme, i.e., w x (L) = (w(x N )−w(x N −1 ))/ h+O(h)). Moreover, we approximate the integral in Eq. (18) via the trapezoidal rule yielding where we have used the discretization (19). Moreover, we concentrate the distributed load f (x, t) onto the closest mesh point by defining the nodal forcing functions The boundary conditions (18) imply that the vertical displacements of the first and last two elements are zero, i.e., w( With the notation (22), a discretized version of the partial differential equation (18) is given bÿ After deriving Eq. (23) as a discretization of the partial differential (19), we will now show that our Theorem 3.1 guarantees the existence of a quasi-periodic response of system (23). We emphasize that our result is especially relevant for Eq. (23). Competing perturbations methods commonly extract a steady-state response from linearizing Eq. (23) around the origin (w =ẇ = 0). Subsequently, those approaches argue that a steady-state response of the linearization implies a close by steady-state response of the full nonlinear system if the nonlinearity is small enough. The nonlinearity in Eq. (23), however, is dominant , since the slenderness parameter r is small. Hence, the validity of regular perturbation methods is restricted to a tiny neighborhood of the origin. Our result 3.1 does not suffer from such a shortcoming.
As noted in Remark 3.3, system (23) can be modified such that condition (C3) holds. In the current setting, a truncation of the nonlinearity for large displacements is justified, since the partial differential equation (18) was obtained by expanding the Green-Lagrange strain tensor in a Taylor series up to order three. Such an approximation is generally only valid for small enough displacements.
To verify condition (C4), we note that Eq. (8) for the discretized beam (23) is given by Since K 1 is a tridiagonal Toeplitz matrix with positive eigenvalues (e.g., Gover [27]), it is positive definite. Thus, the second summand in Eq. (25) is positive for all displacements, i.e., w K 1 ww K 1 w > 0. To analyze the first term in Eq. (25), we split the matrix K 2 as follows where all entries of the matrices B(1) and B(N − 4) are zero expect for the first, respectively, the last entry on the main diagonal which is equal to one. Both matrices (B (1) and B(N − 4)) are positive semi-definite. Furthermore, we note that the square of the symmetric and positive definite matrix K 1 is also positive definite.
Since We calculate a forced response of the discretized beam (23) with the numerical package proposed by Jain et al. [35] and choose N = 10 node points for the spatial discretization of the partial differential equation (18). The arising system (23) has six degrees of freedom. Moreover, we select the parameter values r = 0.01, μ = 0.05 and L = 1 and assume quasi-periodic forcing of the form f (x, t) = 0.06r 2 (sin(2π x/L) + 1)(sin( 1 t) + sin( 2 t)). We sweep the forcing frequencies 1 and 2 close to the first two eigenfrequencies ω 1 ≈ 7.13 rad/s and ω 2 ≈ 19 rad/s and depict the L 2 -norm of the arising quasi-periodic steady-state response in Fig. 6.
The results depicted in Fig. 6 arise from two incommensurate frequencies. Our Theorem 3.1, however, Fig. 6 Quasi-periodic steady-state responses of system (23) with six degrees-of-freedom and non-dimensional parameters μ = 0.05, r = 0.01 and L = 1 and forcing f (x, t) = 0.06r 2 (sin(2π x/L) + 1)(sin( 1 t) + sin( 2 t)). We select the forcing frequencies 1 and 2 in resonance with the first two natural frequencies ω 1 ≈ 7.13 and ω 2 ≈ 19 guarantees the existence of a quasi-periodic steadystate response for any finite number of incommensurate frequencies.
To illustrate the applicability of Theorem 3.1 in the case of more than two incommensurate forcing frequencies, we consider the excitation where ω k denote the eigenfrequencies of the discretized beam (23). For γ close to one, the forcing (27) is in resonance with the first three eigenfrequencies. As previously, Theorem 3.1 guarantees the existence of a quasi-periodic steady-state response also for the forcing (27). Once again, we use the numerical package proposed by Jain et al. [35] and compute a quasiperiodic steady-state response. While the package handles two-dimensional tori fast and efficiently, its memory requirements become challenging and computation times increase significantly for three incommensurate frequencies.
We depict the maximal displacement of the arising quasi-periodic steady-state response of the discretized beam (23) for the forcing (27) in resonance with the first three eigenfrequencies in Fig. 7. Moreover, we include a solution by direct time integration with Matlab's ode45 algorithm serving as a benchmark solution, which we initialize on the computed quasi-periodic orbit at t = 0. While the maximal vibration amplitude of the nonlinear response agrees well  Fig. 7 Quasi-periodic steady-state responses of system (23) with three incommensurate base frequencies (K = 3). We consider system (23) with six degrees-of-freedom, non-dimensional parameters μ = 0.05, r = 0.01 and L = 1 and select the forcing (27) to be in resonance with the first three natural frequencies with the benchmark solution, the linear response differs significantly. Thus, a nonlinear analysis of system (23) is essential to accurately predict the response amplitudes.

Conclusions
We have presented a criterion guaranteeing the existence of a quasi-periodic response for wide class of externally driven multi-degree-of-freedom vibratory systems. Our Theorem 3.1 guarantees the existence of a quasi-periodic orbit for arbitrary large forcing and response amplitudes, whereby it overcomes the major shortcoming of competing perturbation results. Thus, it can provide a priori justification for formal perturbation methods or numerical methods such as harmonic balance. Thereby, the results of such heuristic approaches become valid and computational resources can be efficiently employed. Our criterion requires the damping to act on all coordinates, the stiffness terms to be derived from a potential and satisfy a growths conditions far away from the origin. These conditions are similar to our results in the periodic case [8] except for the growth condition (C4) which is stricter compared to the sign condition in the periodic case. Moreover, the nonlinear terms are required to be globally Lipschitz, which is generally satisfied if Eq. (1) accurately model a physical system (cf. Remark 3.3). We also cover the case of discontinu-ous external forcing, such as square waves or saw-tooth functions (cf. Theorem 3.3).
We utilize our results to revisit a chain of oscillators investigated by Breunung and Haller [7] and a slender beam model proposed by Nayfeh and Mook [49]. Both systems satisfy the conditions of our Theorem 3.1 and hence an approximation of the solution in a Fourier series is justified. We compute the arising steady-state vibrations with the numerical package proposed by Jain et al. [35]. Moreover, we consider a complex loading scenario, where the vibrations are excited by the natural roughness of a road. In this scenario, quasi-periodic functions arise naturally and our Theorem 3.1 justifies subsequent computations with the harmonic balance procedure.
We have limited our discussion to position depending nonlinearities, since these are customary in the vibrations literature. For more complex damping forces or oscillatory systems derived from Lagrangian with position dependent mass matrix arising in, e.g., robotics, an extension to account for velocity dependent nonlinearities is desirable.
After computing a steady-state response, it is important to assess its stability, since only stable solutions are observable. Numerical calculation of the monodromy matrix and calculating its eigenvalues (often referred as Floquet theory) is a well-established and universal procedure to assess the stability of periodic orbits. To the best of our knowledge such a procedure is absent for quasi-periodic orbits.
Furthermore, we have focused on quasi-periodic solutions, since these type of vibrations have been reported for a variety of mechanical and even more general systems (cf. Introduction 1). From a theoretical perspective, an extension of our results to the case of almost-periodic or bounded solutions seems tempting.
Acknowledgements I thank George Haller for fruitful discussions throughout this work and Shobhit Jain for discussing the slender beam example in depth. Furthermore, I apologize to Walter Lacarbonara (Organizer) and Lawrence Virgin (Session Chair) for prematurely announcing these results in the abstract for the first NODYNCON conference before obtaining them.
Funding Open Access funding provided by ETH Zurich.

Availability of data and material
The data shown in the figures are available from the author upon request.

Code availability
The code to compute the data shown in the figures is available from the author upon request.

Declarations
Competing interests The author declares that he has no competing interests.
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/.

A Existence results applied to vibratory systems
In this section, we examine known results guaranteeing the existence of quasi-periodic solutions. Thereby, we demonstrate that the available existence results are inapplicable in a structural dynamics context. We differ between Dafermos' results [18] for first-order systems and the results of Bolt et al. [4] covering second-order systems. We show that both results are inapplicable to the simple single-degree-of-freedom oscillator in the form A.1 First-order systems Dafermos [18] applies his abstract result to the equatioṅ where we restrict our attention to N -dimensional Euclidean spaces. Originally, Dafermos [18] considers more abstract Hilbert spaces. He proves the existence of an weak solution to Eq. (29) which is unique and almost-periodic, if there exists a positive constant α > 0 such that holds. We evaluate condition (30) for the linear oscillator (28) reformulated in the form of Eq. (29) and obtain where we have set y 1 − z 1 =: w 1 and y 2 − z 2 =: w 2 .
Considering the case w 2 = 0, Eq. (31) cannot be satisfied for any α > 0 and w 1 = 0. Hence, Dafermos' results on first-order systems is not relevant in the structural dynamics context.

A.2 Second-order systems
Extending earlier work of Berger and Chen [3], Bolt et al. [4] study systems of the form where the time dependence of all coefficients is continuous and bounded or almost-periodic. Besides technical assumptions on the coefficients in Eq. (32), the main condition of Bolt et al. [4] is requirement that there exists a positive constant α > 0 such that holds. While their result appears to be very general and also guarantees a unique bounded solution, it is inapplicable in the structural dynamics context. This is due to the negative sign in front of the nonlinearity F(x, t) in Eq. (32). For the linear oscillator (28) condition (33) requires where we have set B(t) = 0 and b(t) = c. Equation (34) clearly holds only if the stiffness k is negative. The example considered by Bolt et al. [4] (gyroscopic stabilization) falls exactly in this class. In the structural dynamics context, however, stiffness matrices are commonly positive definite. Hence, the results of Bolt et al. [4] are insignificant for oscillatory systems. For the linear oscillator (28) one can also absorb the damping c into the matrix B(t). Then, Eq. (34) would change and condition (33) requires the linear oscillator (28) to be either overdamped or the stiffness k to be negative. Again, both cases are uncommon in structural engineering.

B Non-existence of a quasi-periodic solution to the undamped-linear oscillator (4)
In the following, we show that for any incommensurate frequency basis with at least two frequencies (i.e., K ≥ 2) there exists a forcing f (φ(t)) such that no quasi-periodic response to system (4) exists. First, we assume resonant forcing, i.e., holds. Then, for the periodic forcing f (t) = cos( κ * , t) system (4) is the undamped harmonic oscillator forced at resonance. Since all solutions of the undamped harmonic oscillator forced at resonance are unbounded, no quasi-periodic solution exists. Next, we focus on the non-resonant case, i.e., holds and construct a forcing such that no quasiperiodic solution to Eq. (4) exists. To this end, we consider the forcing which only depends on the first two frequencies of the vector . The forcing (37) is bounded and continuous.
Assuming the existence of a quasi-periodic solution q * to Eq. (4) we utilize the linearity of system (4) and obtain a steady-state response for every summand of the forcing individually, i.e., q * = ∞ n=1 q * n ,q * n + ω 2 0 q * n = 1 n 2 cos( κ n , t), (38) which yields cos( κ n , t) Then, by superposition the assumed steady-state response of system (4) is given by At the time t = 0 the displacement q * of the assumed quasi-periodic solution is given by In the following, we obtain a lower bound on Eq. (41).
To this end, we derive an upper bound onto the denominator of (41) and rewrite the last part of the denominator as where the integers κ n 1 and κ n 2 denote the first, respectively, second entry of the vector κ n = κ n 1 , κ n 2 , 0, . . . . We use the following Theorem by Minkowski to establish an upper bound onto Eq. (42):

Theorem B.1 For any irrational number θ and any α = mθ +n for all integers n and m, there are infinitely many integers p such that
holds.
Proof For a proof of the above theorem, we refer to Cassels [12] (Theorem IIA p. 48).
To show that Theorem (B.1) applies to Eq. (42), we note that the fraction θ = 1 /| 2 | is irrational, since is incommensurate. Moreover, the non-resonance (36) yields which implies that the second condition of Theorem (B.1) is satisfied. Equation (43) holds for infinitely many integers. Hence, for every n = 1, 2, . . . we can find some | p n | > n such that Eq. (43) holds. We utilize this observation by setting κ n 1 = p n in Eq. (42). Furthermore, we note that min l∈Z | pθ − α − l| in Eq. (43) returns the distance of the irrational number pθ − α to the integer l closest to pθ − α. We select κ n 2 such that sign( 2 )κ n 2 (cf. Eq. (42)) is exactly this integer, i.e., Thus, for Eq. (42) we obtain where we have used Eq. (45) and the assertion of Theorem B.1. We proceed by deriving an upper bound onto the first part of the denominator of Eq. (41) and obtain where we have used the upper bound (46). The bounds (46) and (47) imply for assumed steady-state response evaluated at t = 0 (cf. Eq. (41)) which is diverges. Therefore, the assumed quasiperiodic response for system (4) does not exist and the claim of Fact 2.1 is proven.

C Proof of the main theorem
We base result guaranteeing the existence of a quasiperiodic solutions to system (1), on the following fixed point theorem by Schäfer: Theorem C.1 Let X be a normed space and A a continuous mapping of X into X which is compact on each bounded subset U ∈ X . Then, either (i) the equation has a solution for κ = 1, or (ii) the set of all solutions x to Eq. (49) for 0 < κ < 1 is unbounded.
Proof Smart [59] states and proves the above theorem. We also refer to the original proof by Schäfer [56].  [51] or Smart [59]). Schäfer's fixed point theorem C.1, however, does not require the underlying space X to be complete. As Smart [59] details, the requirement on A to map every bounded set into a compact set implies the existence of a complete metric space.
The key observation is to consider the closure of the space A(X ) instead of X itself. Indeed, assuming that A(X ) is compact implies that the closure of A(X ) (i.e., A(X )) is compact. Moreover, assuming that the operator A maps X onto X , A maps A(X ) onto A(X ) by construction. Hence, we obtain the compact and therefore complete metric space A(X ) and the operator A mapping A(X ) to itself.
Having constructed an underlying complete metric space A(X ), standard versions of either Banach's fixed point theorem (A is a contraction), Brouwer's fixed point theorem (A(X ) ⊂ R N ) or Schauder's fixed point theorem (cf. Smart [59]) formulated in complete spaces can subsequently guarantee the existence of a fixed point. An extensive treatment of continuation theorems of fixed points in normed spaces (not necessarily complete) can be found in Gaines and Mawhin [44] and Mawhin [25]. Before proving Theorem 3.1, we establish some preliminary results. First, we define the relevant metric spaces (cf. section C.1) and we rewrite the mechanical system (1) in the operator form (49) (cf. section C.2). Subsequently, we derive the relevant bounds to show that the second case of Schäfer's theorem does not hold. Thus, by the virtue of Theorem C.1 Eq. (1) has a solution.

C.1 The relevant function spaces
We denote the set of all continuous quasi-periodic functions (cf. Definition 2.1) with C(T K ). The set of quasiperiodic functions with a continuous quasi-periodic time derivative is denoted by C 1 (T K ). Defining the norm ||x(φ(t))|| C 0 = sup the set C 1 (T K ) is a Banach space (cf. Corduneanu [15]), which we will denote by B(C 1 (T K ), || · || C 1 ). Analogously, we define the Banach space B(C(T K ), ||·|| C 0 ), which can also be introduced by considering trigonometric polynomials in the form Then, B(C(T K ), || · || C 0 ) is the closure with respect to the C 0 -norm of all trigonometric polynomials of the form (51) (cf. Corduneanu [15]). Hence, x is an element of B(C(T K ), || · || C 0 ) if and only if it can be uniformly approximated by a series of trigonometric polynomials T L , i.e., there exists trigonometric polynomials of the form (51) such that holds. Moreover, we introduce the map which is a norm on the set C(T K ) (cf. Corduneanu [15]). We will denote the norm (53) by || · || L 2 , which we extend by including the first time derivative as follows Then, the set C 1 (T K ) together with the L 1 2 -norm is a normed space which we denote by N (C 1 (T K ), ||·|| L 1 2 ).
We emphasize that N (C 1 (T K ), ||·|| L 1 2 ) is not complete and hence not Banach with respect to the L 1 2 -norm. Extending C 1 (T K ) by all quasi-periodic functions with finite L 2 -norm and a time derivative bounded in the L 2 -norm, we obtain the set L 1 2 (T K ). We note that L 1 2 (T K ) together with the L 1 2 -norm (54) is a Banach space. Similar to the Banach space B(C(T K ), || · || C 0 ), the space B(L 2 (T K ), || · || L 2 ) can also be introduced by considering the closure of the trigonometric polynomials (51) with respect to the L 2 -norm. Therefore, any function x ∈ L 2 (T K ) can be uniformly (in the L 2norm) approximated by an infinite series of trigonometric functions of the form (51). In summary, all functions in either C 1 (T K ) or L 1 2 (T K ) can be represented by an infinite Fourier series.
Schäfer's Theorem C.1 requires the operator A to be compact on bounded sets. Before establishing compactness of a specific operator arising for the mechanical system (1), we clarify compactness of bounded subsets of the metric spaces B(L 1 2 (T K ), || · || L 1 2 ) and ). We collect our first result in the following proposition: Proof We base the proof of the above Proposition on the following theorem

relatively compact if and only if the following conditions are satisfied
• Y is bounded, i.e., there exists a constant C such that holds for all u ∈ U . • for every ε > 0 there exists a δ > 0 such that holds for all u ∈ U , whenever |δ| < δ.
Proof For a proof of the above theorem, we refer to Precup [51] or Brezis [9].
We show that the above theorem applies to any bounded set U ⊂ B(L 1 2 (T K ), || · || L 1 2 ). Since U is bounded in the L 2 norm by assumption, for any U there exists some C U > 0 such that holds. Equation (57) implies that condition (55) holds. Moreover, we note that any function x ∈ L 2 (T K ) satisfies Parseval's identity Now, we claim that the L 2 -norm along the orbit φ = t on the tours T K is the same as L 2 -norm along any other orbit of the formφ =˜ t if˜ ∈ R K is incommensurate. Indeed, we define Sinceφ is defined on a K -dimensional torus, Parseval's identity holds and we obtain For any function x(φ) ∈ B(L 1 2 (T K ), || · || L 1 2 ), the identity (60) holds also for the velocityẋ(φ). This observation implies d dt x(φ(t)) where we have used the chain rule. The last equality in Eq. (61) implies that the directional derivative of x in the direction˜ is bounded in the L 2 -norm. We rely on this observation to verify condition (56) of Theorem C.2. We permute to generate K linearly independent vectors of the form and collect the vectors k as columns of the matrix := [ 1 , 2 , . . . , K ]. By construction, the vectors k are linearly independent and hence the matrix is invertible, i.e., || −1 || < C holds for some finite constant C > 0. We observe that the vector δ in condition (56) can be represented as a linear combination of the basis vectors k , i.e., Equation (63) also includes an upper bound onto the entries η k of the vector η. We substitute Eq. (63) into condition (56) and obtain where we have used the Cauchy-Schwarz inequality.
Following Proposition (C.1), we have the following Corollary .
By Proposition C.1Ũ is compact. Since a closed subset of a compact space is compact (cf. Lang [41]), is compact.

C.2 The operator A
To identify the operator A of the homotopy (49), we rewrite the mechanical system (1) into the form (49).
To this end, we note that condition (C4) implies that the continuous function q S(q) is either positive or negative for |q| > r . Depending on the sign of q S(q) for large displacements, we define the following matrix Then, we consider the following homotopy which is equivalent to system (1) if κ equals one. To transform the differential equation (68) into the operator form (49), we collect properties of the continuous stiffness terms into the following proposition: Proposition C.2 Assume the stiffness terms S(q) are continuous. Then, S maps C(T K ) to C(T K ) continuously with respect to the C 0 -norm.
Proof By assumption S(q(φ)) maps the torus T K to R N continuously, i.e., S : C(T K ) → C(T K ) holds.
We have the analogous statement for the L 2 case, if the global Lipschitz condition (C3) holds.

Proposition C.3 Assume that the stiffness terms S(q)
satisfy condition (C3). Then, S maps L 2 (T K ) to L 2 (T K ) and is continuous with respect to the L 2 -norm.
Proof First, we note that S(q(φ)) maps the torus T K to R N , hence S(q(φ)) is quasi-periodic with the same frequency base as the coordinates q. For the L 2norm, we obtain S(q(φ)) 2 where we have used Eq. (7). For continuity, we select Equation (70) implies that S(q(φ)) is continuous for any q(φ) ∈ L 2 (T K ). This completes the proof of Proposition (C.3).
We apply Schäfer fixed point Theorem C.1 in the normed spaces N (C 1 (T K ), || · || L 1 2 ) and B(L 1 2 (T K ), || · || L 1 2 ). In both spaces the coordinate q( t) can be represented as a trigonometric series with coefficients q κ l . Then, Propositions C.3 and C.2 establish that S(q) can be expanded in an infinite trigonometric series of the form (51) with the same frequency base as the coordinate q( t). We denote the coefficients of the trigonometric series approximating S evaluated along q(φ) by S κ l (q) and rewrite Eq. (68) as where we have used the definitions (67) and (72). We denote the minimal entry of the diagonalized damping matrix by D min := min 1≤ j≤N (|D j j |) > 0, which is greater than zero due to condition (C1). With this notation, we derive an upper bound onto the inverse of the diagonal matrixH(κ l ) by noting where we have differentiated between low frequencies 2 κ l , 2 < 1 and higher frequencies. In summary, the inverse of the matrixH(κ l ) is bounded from above by max(4, 2/(D 2 min )). For the later purpose, we derive an upper bound onto the norm of || κ l , (H(κ l )) −1 || by the following The bound (75) implies that the inverse of the transfer functions H(κ l ) is well defined for any excitation frequency κ l , , i.e., we obtain Similarly, for || κ l , (H(κ l )) −1 || we obtain Equation (77) implies that the inverse of H(κ l ) exists for all κ l ∈ Z K , therefore we can rewrite Eq. (71) as which is in the form (49) by defining the operator After identifying the operator A, we establish that A satisfies the requirements of Schäfer's theorem.
To this end, we proceed by showing that the opera- ). Moreover, we establish continuity in the L 1 2 -norm and compactness on bounded subsets. By inspection the right-hand side of Eq. (80) can be written as a Fourier series with the base frequencies . Hence, the operators maps quasi-periodic functions to quasi-periodic functions.
For the normed space N (C 1 (T K ), ||·|| L 1 2 ) we collect our results in the following lemma: Proof We derive the following bounds in the C 0 -norm where we have used the bound (77). Similarly, for the time derivative we obtain where we have used the bound (78). We note that the right-hand side of Eqs. (81) and (82) are bounded for any q ∈ C(T K ). Hence, A maps C 1 (T K ) to C 1 (T K ). Furthermore, selecting q 1 − q 2 L 1 which implies continuity of the operator A in the L 1 2norm.
To show compactness for any bounded set U , we apply the Fréchet-Kolmogorov theorem C.2. We denote the constant bounding the set U in the L 1 2 -norm by C U . Then, for the condition (55) we have where we have used Proposition C.3, i.e., the fact that S(q) is bounded in the L 2 -norm (cf. Eq. (69)). The upper bound (84) implies that condition (55) of Theorem C.2 is satisfied. For condition (56), we obtain where we have used the upper bound (66) in the proof of Proposition C.1. Selecting |δ| < δ = ε/ (h 1 (||K|| + L) K C C U ) Eq. (85) implies that condition (56) holds. Hence, A is compact on any bounded subset of N (C 1 (T K ), ||·|| L 1 2 ). This completes our proof of Lemma C.1. Proof We begin by showing that A maps B(L 1 ). To this end, we obtain where we have used the bounds (77) ).
The remainder of Lemma C.2, i.e., continuity in the L 1 2norm and compactness on bounded sets, is analogous to the proof of Lemma C.1.

C.3 Bound on the solution
In this section, we show that all solutions to Eq. (49) with the operator (80) are bounded. To his end we note that any quasi-periodic function in C(T K ) is also almost-periodic, i.e., the following holds Property C.1 (Almost-periodic in the sense of Bohr) The continuous and bounded function f(t) is almostperiodic, if for any t 0 , τ ∈ R and ε > 0 there exists T = T (ε) > 0, such that holds (cf. Corduneanu [15]).
For the periodic case T (ε) is the period. Therefore, T (ε) is also considered as quasi-period. Property (C.1) only holds for continuous quasiperiodic functions. To make use of property (C.1) in the non-smooth case, we approximate the nonsmooth functions q ∈ L 1 2 (T K ) by continuous functions. Indeed, any function q ∈ L 1 2 (T K ) can be approximated arbitrarily well by a continuous function, e.g., the trigonometric polynomial of the form (51). Therefore, for any solution q ∈ L 1 2 (T K ) to Eq. (49) we can find some p ∈ C 1 (T K ) such that holds. With the notation f * (p, q, κ) := p − q + κ(A(q) − A(p)), we observe that p approximately satisfies the operator equation (49) Equation (87) implies the error term f * is small for all q ∈ Ł 1 2 (T K ), i.e., we have the following upper bound ||f * (p, q, κ)|| L 2 Rewriting the homotopy (88) as a differential equation, we obtain (91) In the continuous case, i.e., for q ∈ N (C 1 (T K ), || · || L 1 2 ) the approximation with smooth functions is not necessary and hencef(t) = f(t) holds.

Lemma C.3 Assume that conditions (C1) and (C2)
hold and f(t) ∈ L 2 (T K ). Then, for all solutions p ∈ C 1 (T K ) to (90), the upper bound holds, where δ can be made arbitrarily small.
Proof We begin by defining the following scalar quasiperiodic auxiliary functions Then, the function is quasi-periodic for any fixed t 0 . Thus, property C.1 implies that inside any interval of length T (ε) the function ξ(t, t 0 ) is ε-close its value at t 0 , i.e., zero. We collect these time instances in the set Q(p, t 0 , ε) To obtain an upper bound onto the velocities of sys-Following Lemma C.3, we derive an upper bound onto the L 2 -norm of the positions p of all quasi-periodic solutions to Eq. (90). To this end, we denote the maximum value of the continuous functions |p S(p)| in the bounded interval of |p| < r by S m := sup |p|≤r |p S(p)|.
Furthermore, the positive definite mass matrix satisfies for some positive constants C * M ≥ C M > 0. Finally, we denote the minimum of the two positive constants C M and α (cf. condition (C4)) by β, i.e., β := min(α, C M ) > 0. (108) With the introduced notation (106), (107) and (108) With the definition of β (cf. Eq. (108)) we refine the lower bounds (113), respectively, (115) to fixed point theorem C.1 applies. The upper bounds stated in Lemmata C.3 and C.4 and Eq. (87) guarantee, that the set of solutions to Eq. (49) is bounded by where ε is an arbitrarily small constant. Thus, the second case of Schäfer's fixed point does not hold and hence Eq. (68) has a solution for κ = 1. Since Eq. (68) is equivalent to system (1) for κ = 1, system (1) has a solution bounded in L 1 2 (T K ) by Eq. (122). This proves the claim of Theorem C.3, which is equivalent to Theorem 3.3.

C.5 Existence of smooth solutions
If the forcing is continuous, then the result from Theorem C.3 can be strengthened to guarantee a continuously differentiable quasi-periodic solution.
Theorem C. 4 Assume that conditions (C1)-(C4) are met and the forcing is quasi-periodic and continuous, i.e., f ∈ C(T K ). Then, system (1) has a solution in C 1 (T K ).
Moreover, Lemma C.1 guarantees that A maps N (C 1 (T K ), || · || L 1 2 ) to N (C 1 (T K ), || · || L 1 2 ), is continuous and compact on bounded sets. Thus, all requirements of Schäfer's fixed point theorem C.1 are met. The upper bounds provided by Lemmata C.3 and C.4 imply that the set of solutions to Eq. (49) is bounded by where the constant δ in the bounds (92) and (109) can be set to zero. Once again, the second case of Schäfer's fixed point theorem C.1 does not hold. Hence, Eq. (68) has a solution for κ = 1. Since (68) is equivalent to system (1) for κ = 1, system (1) has a solution in C 1 (T K ). More specifically, there exists a solution q * to Eq. (1) for which q * ∈ C 1 (T K ), and ||q * || L 1 holds. This proves Theorem C.4, which is a restatement of Theorem 3.1.
Remark C.2 Our Theorem (C.4) guarantees that at least one solution to Eq. (1) is continuous and bounded.
With the knowledge of such a solution q * ∈ C 1 (T K ), we can apply Schäfer's fixed point theorem C.1 in the Banach space B(C 1 (T K ), · C 1 ). Since the first case holds (i.e., q * exists), the second cases does not hold, i.e., the set of solutions to the homotopy (49) is bounded in the C 1 -norm.

D Proof of Theorem 3.2
To proof Theorem 3.2, we show that condition (C4*), implies that condition (C4) holds. To this end, we assume that the nonlinearity S(q) evaluated at q = 0 is zero. We can make this assumption without the loss of generality. First, we assume that the Hessian is positive definite for q >r and note that within the bounded domain |q| <r there exists a constant 0 < C min < ∞ such that the real and symmetric Hessian matrix ∂ 2 V (q)/∂q 2 satisfies |q| <r , In the following, we show that Eq. (11) ensures that Eq. (8) holds with r = 2r 1 + C min C S , and α > C S /2, where the constant C S is a lower bound on the magnitude of the Hessian for |q| > r (cf. Eq. (11)). To this end, a Taylor series expansion of the nonlinearity yields q S(q) = q S(0) where we have denoteds =r /|q| as the value for which the argument of the Hessian matrix sq leaves the bounded domain |q| <r . Due to the choice of the radius r in Eq. (126),s is upper bounded bỹ s ≤ C S 2 (C S + C min ) ≤ 1.
For |q| >r , i.e., for the second integral in Eq. (127), the lower bound (11) from Theorem 3.2 holds, whereas for the first integral the lower bound (125) holds. Together these bounds imply for Eq. (127) q S(q) > (−sC min + (1 −s)C S )|q| 2 where we have used the lower bound on the parameter s from Eq. (128). Equation (129) implies, that condition (8) is satisfied for the radius defined in Eq. (126) and α = C S /2 > 0. The proof in the case of a negative definite Hessian matrix for |q| > r is analogous. The only difference is that instead of a lower bound on the quadratic form (125) one needs to work with an upper bound of the quadratic form (125) inside a bounded domain.

E Theorem 3.1 applies to the oscillator chain
In the following, we show that the chain of oscillators depicted in Fig. 1 satisfies the conditions of our existence Theorem 3.1. In our previous work [8], we have shown that the damping matrix is positive definite and have derived a potential for the stiffness terms of Eq. (13). Thus, conditions (C1) and (C2) of Theorem 3.1 hold.
For the Lipschitz condition (C3), we point out Remark 3.3 where we have constructed a Lipschitz continuous truncation of system (1). Selecting the radius of the truncation large, i.e., in the overflow our computation package (Matlab), this truncation does not necessarily need to be implemented.