# Non-smooth modal analysis of piecewise-linear impact oscillators

Anders Thorin, Pierre Delezoide, Mathias Legrand

► To cite this version:

Anders Thorin, Pierre Delezoide, Mathias Legrand. Non-smooth modal analysis of piecewise-linear impact oscillators. SIAM Journal on Applied Dynamical Systems, 2017, 16 (3), pp.1710-1747. 10.1137/16M1081506 . hal-01298983v2

HAL Id: hal-01298983

https://hal.science/hal-01298983v2

Submitted on 8 Aug 2017

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Non-smooth modal analysis of piecewise-linear impact oscillators

Anders THORIN¹, Pierre DELEZOIDE², Mathias LEGRAND¹

¹ Department of Mechanical Engineering, McGill University, Canada

² Independent researcher, France

**ABSTRACT** Periodic solutions of autonomous and conservative second-order dynamical systems of finite dimension *n* undergoing one unilateral contact condition are investigated in continuous time. The unilateral constraint is complemented with a purely elastic impact law which preserves total energy. The dynamics is linear when there is no contact. The number *k* of impacts per period arises as a natural parameter of the proposed formulation.

Interestingly, the existence of the targeted periodic solutions is essentially governed by a system of only *k* − 1 nonlinear equations with *k* unknowns, regardless of the number of degrees-of-freedom. This serves to prove that the phase-space is populated by one-dimensional continua of periodic solutions generating invariant manifolds which can be understood as *nonsmooth modes of vibration* in the context of vibration analysis. Additionally, these equations provide an efficient and systematic way of calculating nonsmooth modes of vibration. They also demonstrate the existence of interesting properties: symmetries of trajectories and emergence of unique features such as manifolds supporting constant-frequency orbits. All results are illustrated through a simple in-line spring–mass system whose last mass undergoes the unilateral impact law. Stability is briefly discussed and a few neutrally stable modes are depicted.

**KEYWORDS** impact dynamics, nonlinear normal modes, vibro-impact oscillators, nonlinear oscillators, nonlinear modal analysis, nonsmooth dynamics

<table><tr><td>1</td><td>Background</td><td>1</td></tr><tr><td>2</td><td>Assumptions and formulation</td><td>2</td></tr><tr><td>2.1</td><td>General framework</td><td>2</td></tr><tr><td>2.2</td><td>Formulation for periodic solutions</td><td>4</td></tr><tr><td>3</td><td>Main results</td><td>5</td></tr><tr><td>3.1</td><td>Preliminary notations</td><td>5</td></tr><tr><td>3.2</td><td>Important theorems</td><td>6</td></tr><tr><td>4</td><td>Proof of Theorems 3.1, 3.2, 3.4</td><td>8</td></tr><tr><td>4.1</td><td>Free flight dynamics</td><td>8</td></tr><tr><td>4.2</td><td>Formulation of the impact law</td><td>9</td></tr><tr><td>4.3</td><td>Sequence of impacts</td><td>9</td></tr><tr><td>4.4</td><td>Initial conditions invariant after k successive impacts (proof of Theorem 3.1)</td><td>10</td></tr><tr><td>4.5</td><td>Gap closure at impact times (proof of Theorem 3.2)</td><td>11</td></tr><tr><td>4.6</td><td>Explicit expression of the positions and velocities (proof of Theorem 3.4)</td><td>11</td></tr><tr><td>5</td><td>Dimension of the solution space and nonsmooth modes [proof of Theorem 3.3]</td><td>12</td></tr><tr><td>6</td><td>Symmetries in the modal orbits</td><td>12</td></tr><tr><td>6.1</td><td>Symmetry of a k-periodic sequence</td><td>13</td></tr><tr><td>6.2</td><td>Extension of the impact times sequence to Z</td><td>13</td></tr><tr><td>6.3</td><td>Symmetries of positions in the time domain</td><td>14</td></tr><tr><td>6.4</td><td>Dimension of the solution space in the symmetric case</td><td>15</td></tr><tr><td>7</td><td>Illustration on a simple serial oscillator</td><td>16</td></tr><tr><td>7.1</td><td>Simplified governing equations</td><td>16</td></tr><tr><td>7.2</td><td>One-impact-per-period solutions</td><td>18</td></tr><tr><td>7.3</td><td>Two-impact-per-period solutions</td><td>19</td></tr><tr><td>7.4</td><td>Multiple-impact-per-period solutions</td><td>25</td></tr><tr><td>8</td><td>Insight on stability of nonsmooth modes</td><td>26</td></tr><tr><td>9</td><td>Concluding remarks</td><td>27</td></tr><tr><td>10</td><td>References</td><td>27</td></tr><tr><td>A</td><td>Dimensions of modes in the symmetric case</td><td>29</td></tr></table>

## 1. Background
In linear dynamics, the solution space is a linear space spanned by linear modes of vibration. For nonlinear systems governed by Ordinary Differential Equations, the above linear space structure no longer holds and methods able to systematically characterize the global dynamics are not available. However, it is known that fixed points, also called equilibrium points, limit cycles, and chaos may exist. The centre manifold theorem [13, 16] together with Lyapunov's center theorem [6, p. 5] prove that for second-order dynamics away from resonance conditions, two-dimensional manifolds of the phase space exist and are tangent to the linear modes of the system linearized in the vicinity of its fixed points. Such invariant two-dimensional manifolds in the phase space were later defined as *nonlinear modes of vibration* [37, 40] in the vibration community.
Nonlinear modes have been primarily applied to smooth systems, i.e. systems described by ODEs. Though not as informative as their linear counterparts, nonlinear modes still reveal significant features such as potential vibratory resonances, stability properties, and step up the ground for modal reduction [40]. However, many mechanical systems, referred to as nonsmooth systems, involve unilateral constraints inducing velocity discontinuities and cannot be described with ODEs; existing methods can therefore not be used to compute their nonlinear modes. Yet, the vibratory characterization of nonsmooth systems has many applications. For instance, rotor-stator contact interactions in rotating machinery involves unilateral contact occurrences between the casing and the blades [42], which needs to be understood for obvious safety and maintenance reasons. Optimal progression rate of vibro-impact systems such as pile driving or percussive drilling [31] or, on a smaller scale capsule systems (capsubots) [23], rely on capturing the nonsmooth vibratory dynamics. The dynamics of impact dampers implemented to reduce vibrations should also be properly understood [21]. Supplementary applications can be found in [4].

The present work is a contribution towards the vibration analysis of such nonsmooth systems. It focuses on the following subclass of nonsmooth systems: *n*-degree-of-freedom linear oscillators subject to a single unilateral contact condition complemented with an elastic impact law. Since modal analysis is targeted, the autonomous and conservative dynamics is considered.

Similar systems, often referred to as vibro-impact oscillators, have been investigated for several decades from the standpoint of vibratory analysis. A majority of such works consists either in i) regularizing the impact condition(s) with the addition of springs [4, 7, 15, 19, 27, 28, 36, 41] or ii) investigating the dynamics of a forced vibro-impact oscillator [12, 34, 35, 43, 44]. Regularizing is known to introduce numerical stiffness [1, 2] and is not considered in this work (the relationship between regularized and infinitely rigid impact was studied using a two-degree-of-freedom system with an elastic stop with one impact per period in ref. [30]).

Existence results for such nonsmooth systems are provided in [29] for grazing bifurcations as well as in [32] for oscillators with bilateral rigid barriers where each contact is activated at most once per period. Dynamic behaviour in the neighbourhood of grazing trajectories is investigated in [8, 11]. Most occurrences in the literature of the terminology multiple impacts are related to impact oscillators with several unilateral constraints (see some of the references above or [24, 26]); here, multiple impacts per period refers to the same contact condition being activated several times during one period. Periodic motions with multiple impacts per period have sometimes been mentioned for forced systems [17, 22, 25, 33] but never systematically explored. Recently, a spring-mass system subject to a purely elastic impact law with one impact per period has been thoroughly investigated [20]. The present work significantly extends the results of the latter—thorough stability analysis set aside—to more complex geometries, non-diagonal mass matrices (such as in Finite-Element models) and multiple impacts per period, by providing a general analytic and minimal expression of periodic solutions for such systems for any number of $k \in \mathbb{N}^*$ impacts per period.

Because the dynamics is deterministic, the problem is formulated as finding the initial conditions $\mathbf{x}_0 \in \mathbb{R}^{2n}$ from which emanates *T*-periodic solutions which satisfies a succession of “free flights” governed by linear dynamics, punctuated by *k* impacts. Enforcing periodicity boils down to finding an initial condition invariant by a linear mapping which is explicitly exhibited. Although challenging, it nicely reduces to solving two independent systems of equations, nonlinear in the *k* times of impacts and linear in *k* scalar unknowns, irrespective of the number of degrees-of-freedom *n*. These *k* scalar quantities then determine potential initial conditions $\mathbf{x}_0$ via an explicitly known isomorphism; $\mathbf{x}_0$ is a solution of the problem when it satisfies a simple (but restrictive) inequality. The assumptions and the general problem of interest are first described (Sec. 2). Then, the main results are succinctly provided (Sec. 3) and proven (Sec. 4 to 6). Results are illustrated on a simple system (Sec. 7) and stability is briefly discussed (Sec. 8). Animations of motions with *k* impacts per period are available online using the hyperplink provided in [39].

**Remark 1.1.** Nonlinear modes of discrete systems are, by definition, continua of periodic orbits of a nonlinear system. By extension, nonsmooth modes are defined as nonlinear modes of nonsmooth systems. In the same spirit, linear modes can be conceived as continua of periodic elliptical orbits of a linear system.

## 2. Assumptions and formulation
## 2.1. General framework

The following assumptions are considered throughout this work: the n-dof mechanical system of interest lies in a one, two or three-dimensional space. No damping and no external forces are considered. The system is constrained by a single unilateral condition preventing penetration in a rigid still obstacle. The system is assumed to be linear when contact is not activated. In particular, displacements are small so that the unilateral constraint is expressed in its linearized form.

Let $\mathbf{u} \in \mathbb{R}^n$ be the vector of generalized positions and $\mathbf{v} = \dot{\mathbf{u}}$ be the generalized velocities. As a general principle, subscripts $\bullet_1$ indicate quantities related to $\mathbf{u}$ while $\bullet_2$ indicates $\mathbf{v}$-related quantities. The gap $g(\mathbf{u})$, i.e. the signed distance function separating the contacting dof and the rigid obstacle, is chosen positive when there is no contact. Its time-derivative, the contact velocity, is denoted $\gamma$. The unilateral constraint preventing interpenetration writes

$$g(\mathbf{u}) \geq 0 \qquad (2.1)$$

which does not explicitly depend on time since the obstacle does not move. Because it is assumed linear in $\mathbf{u}$ there exists $\mathbf{w}_1 \in \mathbb{R}^n$ such that

$$g(\mathbf{u}) = \mathbf{w}_1^\top \dot{\mathbf{u}} + g_0 \qquad (2.2)$$

where the real number $g_0$ stands for the signed distance between the unconstrained rest position of the system and the position of the obstacle, see Figure 1. Accordingly:

$$\gamma = \frac{dg(\mathbf{u})}{dt} = \mathbf{w}_1^\top \dot{\mathbf{u}} = \mathbf{w}_1^\top \mathbf{v}. \qquad (2.3)$$

<figure><img src="image_2.png" /><figcaption>Figure 1: Positions of the obstacle illustrated in a one-dimensional case.</figcaption></figure>

At a given time, exactly two distinct configurations are possible:

1. Contact is not activated, that is $g(\mathbf{u}) > 0$ and the dynamics obeys

$$M\ddot{\mathbf{u}} + K\mathbf{u} = 0 \qquad (2.4)$$

where $M \in \mathbb{R}^{n,n}$ is the mass matrix and $K \in \mathbb{R}^{n,n}$, the stiffness matrix; $M$ and $K$ are symmetric positive-definite. This situation is referred to as “free flight”. By convention, the unconstrained rest positions are chosen to be zero.

2. Contact is activated, that is $g(\mathbf{u}) = 0$ and the dynamics is governed by [38]:

$$M(\mathbf{v}^{+} - \mathbf{v}^{-}) = \mathbf{w}_{1} \Lambda \qquad (2.5)$$

where $\Lambda$ is the contact impulse along the outward-pointing normal direction of the contact surface, $\mathbf{v}^{-}$ and $\mathbf{v}^{+}$ are the pre- and post-impact velocities: if $t$ is a time of impact, $\mathbf{v}^{+}(t) = \lim_{\epsilon \to 0^{+}} \mathbf{v}(t+\epsilon)$ and $\mathbf{v}^{-}(t) = \lim_{\epsilon \to 0^{-}} \mathbf{v}(t+\epsilon)$. Since this work involves a single unilateral constraint, $g(\mathbf{u}) \in \mathbb{R}$ and $\mathbf{w}_1$ is a column vector. Sticking is not considered, so necessarily $\Lambda \geq 0$. The impulse $\Lambda$ is the integral of a generalized contact force $\lambda$.

In contrast to the space-continuous framework for which existence and uniqueness hold without any impact law [10, 18], discrete mechanical systems do not incorporate the constitutive description of the impact process and the formulation should be supplemented with an *impact law* [5, 9]. It maps $(g(\mathbf{u}), \gamma^{-})$ to $\gamma^{+}$ [5] where $\gamma^{-}$ and $\gamma^{+}$ are the pre- and post-impact normal contact velocities, respectively. Since unforced periodic solutions are targeted, the only possible choice is a conservative impact law of the form

$$g(\mathbf{u}) = 0 \quad \Rightarrow \quad \gamma^{+} = -\gamma^{-}. \qquad (2.6)$$
which can also be expressed in terms of the generalized velocities $\mathbf{v}^-$ and $\mathbf{v}^+$ through the introduction of the matrix $\mathbf{N}_2 \in \mathbb{R}^{n,n}$:

$$g(\mathbf{u}) = 0 \implies \mathbf{v}^+ = \mathbf{N}_2 \mathbf{v}^- \qquad (2.7)$$

Proposition 2.1, proven in Sec. 4.2, shows that $\mathbf{N}_2$ is fully determined by the mass matrix $\mathbf{M}$ and the gradient $\mathbf{w}_1$ of the gap.

**Proposition 2.1** *The impact law (2.7) is such that*

$$\mathbf{N}_2 = \mathbf{V}^{-\mathrm{T}}(\mathbf{I} - 2\mathbf{r}_2\mathbf{r}_2^{\mathrm{T}})\mathbf{V}^{\mathrm{T}} \qquad (2.8)$$

where $\mathbf{r}_2$ is the unit vector defined by $\mathbf{r}_2 = \mathbf{V}^{-1}\mathbf{w}_1(\mathbf{w}_1^{\mathrm{T}}\mathbf{M}^{-1}\mathbf{w}_1)^{-1/2}$ and $\mathbf{V}$ comes from the Cholesky decomposition $\mathbf{M} = \mathbf{V}\mathbf{V}^{\mathrm{T}}$.

**2.2. Formulation for periodic solutions** We focus on solutions corresponding to $t \mapsto g(\mathbf{u}(t))$ of finite support, excluding phenomena such as chattering or sticking. This naturally introduces two quantities of interest: the unknown period $T$ and the (finite) number of impacts per period (ipp) $k$. Solutions with $k=0$ are smooth solutions of the linear problem and are of course not of interest here, hence $k \in \mathbb{N}^*$. Without loss of generality, contact is assumed to be activated at $t=0$, i.e. $g(\mathbf{u}(0)) = 0$. With $T$ the unknown period of a solution $\mathbf{u}$, the times of impact in $[0, T]$ are denoted by $t_1, \dots, t_k$ with $t_k = T$ and, by convention, $t_0 = 0$. The dynamics is governed by Problem (2.9) expressed as: Given $T$ and $k$, find $\mathbf{u}$, $\mathbf{v} = \dot{\mathbf{u}}$ and $t_1, \dots, t_k \in [0, T]$ satisfying

$$\left\{
\begin{array}{ll}
\forall i \in [1, k], \quad \forall t \in (t_{i-1}, t_i), \quad \mathbf{M}\ddot{\mathbf{u}}(t) + \mathbf{K}\mathbf{u}(t) = \mathbf{0} & \text{free flight,} \qquad (2.9a) \\
\forall i \in [1, k], \quad g(\mathbf{u}(t_i)) = 0 & \text{gap closure,} \qquad (2.9b) \\
\forall t \in [0, T], \quad g(\mathbf{u}(t)) \ge 0 & \text{non-penetration,} \qquad (2.9c) \\
\forall i \in [1, k], \quad \mathbf{v}^+(t_i) = \mathbf{N}_2\mathbf{v}^-(t_i) & \text{impact law,} \qquad (2.9d) \\
\mathbf{u}(0) = \mathbf{u}(T), \quad \mathbf{v}(0) = \mathbf{v}(T) & \text{periodicity.} \qquad (2.9e)
\end{array}
\right.$$

The contact force $\lambda$ is determined at the times of impact $t_1, \dots, t_k$ through Eq. (2.5) and vanishes everywhere else.

**Definition 2.1** *Nonsmooth modes of vibration are continuous families of periodic solutions to Problem (2.9).*

**2.2.1. Formulation in terms of initial conditions** Being deterministic, Problem (2.9) can be written in terms of initial positions $\mathbf{u}(0)$ and initial velocities $\mathbf{v}(0)$. The motion consists in a succession of $k$ free flights of linear dynamics, separated by $k$ impacts dictated by the impact law. Introducing $\mathbf{x} = [\mathbf{u}^{\mathrm{T}}, \mathbf{v}^{\mathrm{T}}]^{\mathrm{T}}$, the free dynamics between arbitrary times $\tau_1$ and $\tau_2$ can be written in the form $\mathbf{x}(\tau_2) = \tilde{\mathbf{S}}(\tau_2 - \tau_1)\mathbf{x}(\tau_1)$. The matrix $\tilde{\mathbf{S}}(t) \in \mathbb{R}^{2n, 2n}$ is derived explicitly in Sec. 4. Similarly, the impact law can be written as a mapping of $\mathbf{x}$:

$$g(\mathbf{u}) = 0 \quad \implies \quad \mathbf{x}^+ = \tilde{\mathbf{N}}\mathbf{x}^- \quad \text{with} \quad \tilde{\mathbf{N}} = \begin{bmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{N}_2 \end{bmatrix}. \qquad (2.10)$$

The overall dynamics results in the following composite linear mapping acting on $\mathbf{x}_0 := \mathbf{x}(0)$:

$$\left\{
\begin{array}{ll}
\forall i \in [2, k], \quad \forall t \in [t_{i-1}, t_i), \quad \mathbf{x}(t) = \tilde{\mathbf{S}}(t - t_{i-1})\tilde{\mathbf{N}} \dots \tilde{\mathbf{N}}\tilde{\mathbf{S}}(t_1 - t_0)\mathbf{x}_0 & (2.11a) \\
\mathbf{x}(t_k) = \tilde{\mathbf{N}}\tilde{\mathbf{S}}(t_k - t_{k-1})\tilde{\mathbf{N}} \dots \tilde{\mathbf{N}}\tilde{\mathbf{S}}(t_1 - t_0)\mathbf{x}_0. & (2.11b)
\end{array}
\right.$$

with $t_0 = 0$.

The question of finding periodic solutions to Problem (2.9) with impact times $t_1, \dots, t_k$ is formulated in three necessary conditions (NC) to be satisfied by $\mathbf{x}_0$: Given $T \in \mathbb{R}_+^*$ and $k \in \mathbb{N}^*$, find $\mathbf{x}_0 \in \mathbb{R}^{2n}$ and $t_1, \dots, t_k$ such that:

$$\left\{
\begin{array}{ll}
\tilde{\mathbf{N}}\tilde{\mathbf{S}}(t_k - t_{k-1})\tilde{\mathbf{N}} \dots \tilde{\mathbf{N}}\tilde{\mathbf{S}}(t_1 - 0)\mathbf{x}_0 = \mathbf{x}_0 & \text{NC1 (dynamics + periodicity)} \qquad (2.12a) \\
\mathbf{x} \text{ determined from } \mathbf{x}_0 \text{ through (2.11) satisfies:} & \\
\forall i \in [1, k], \quad \tilde{g}(\mathbf{x}(t_i)) = 0 & \text{NC2 (gap closure)} \qquad (2.12b) \\
\forall t \in [0, T], \quad \tilde{g}(\mathbf{x}(t)) \ge 0 & \text{NC3 (non-penetration)} \qquad (2.12c)
\end{array}
\right.$$
with $\tilde{g}(\mathbf{x}) := g(\mathbf{u})$. Eq. (2.12a) enforces periodicity through a sequence of free flights over $(t_{i-1}, t_i)$ punctuated by perfectly elastic impacts at $t_i$, $i \in [1, k]$. Eq. (2.12b) guarantees that the impact law applies only when the separating gap is closed. Ineq. (2.12c) is necessary to ensure that the unilateral contact condition is not violated during the free flights.

The fact that NC3 is an inequality while NC1 and NC2 are equalities introduces the distinction between “potential solutions” and “admissible solutions” defined below.

**Definition 2.2 [Potential solution]** A potential solution is an initial condition $\mathbf{x}_0$ which satisfies NC1 and NC2 for some $t_1, \dots, t_k$. By extension, a time-function $\mathbf{x}$ determined by a potential solution $\mathbf{x}_0$ using relations (2.11) also qualifies as such.

**Definition 2.3 [Admissible solution]** An admissible solution is a potential solution which satisfies NC3.

# 3. Main results

**3.1. Preliminary notations** Efficiently exposing the main results necessitates a few quantities which require various definitions. The resulting heavy notations must not mask the fact that their computations are very inexpensive and they are all determined solely by the mass matrix **M**, the stiffness matrix **K**, the gap **g** and the times of impact $t_1, \dots, t_k$. Let us introduce the following quantities:

* **V** defined by the Cholesky decomposition of $\mathbf{M} = \mathbf{V}\mathbf{V}^\top$.

* An orthogonal matrix **Q** such that $\mathbf{V}^{-1}\mathbf{K}\mathbf{V}^{-\top} = \mathbf{Q}\mathbf{\Omega}^2\mathbf{Q}^\top$ where $\mathbf{\Omega} = \text{diag}(\omega_1, \dots, \omega_n)$.

* The period $T \in \mathbb{R}_+^*$, assumed not to be a multiple of a $2\pi/\omega_i$, that is $\forall i \in [1, n]$, $T \neq 0$ (mod $2\pi/\omega_i$).

* The number $k \in \mathbb{N}^*$ of impacts per period and $(t_1, \dots, t_k) \in \mathbb{R}_+^* k$ the times of impact with $t_{i-1} < t_i$. The duration between two successive impacts is denoted by $\sigma_i = t_i - t_{i-1}$ with $t_0 = 0$ and $T = t_k$.

* The matrix **P** and vectors **r**, **w** defined by

$$ \mathbf{P} = \begin{bmatrix} \mathbf{\Omega}\mathbf{Q}^\top\mathbf{V}^\top & \mathbf{0} \\ \mathbf{0} & \mathbf{Q}^\top\mathbf{V}^\top \end{bmatrix}, \quad \mathbf{r} = \begin{bmatrix} \mathbf{0} \\ \mathbf{Q}^\top\mathbf{r}_2 \end{bmatrix}, \quad \mathbf{w} = \mathbf{P}^{-\top}\begin{bmatrix} \mathbf{w}_1 \\ \mathbf{0} \end{bmatrix}. \qquad (3.1a) $$

* The matrix $\mathbf{S}(t)$ whose four blocks are diagonal:

$$ \mathbf{S}(t) = \begin{bmatrix} \cos(\mathbf{\Omega}t) & \sin(\mathbf{\Omega}t) \\ -\sin(\mathbf{\Omega}t) & \cos(\mathbf{\Omega}t) \end{bmatrix} \qquad (3.1b) $$

* The skew-symmetric matrix $\mathbf{\Pi}$ populated by the elements

$$ \mathbf{\Pi}_{i,j}(t_1, \dots, t_k) = \Upsilon_T(2(t_j - t_i) - T), \quad i < j \leq k. \qquad (3.1c) $$

where $r_{n+i}$ denotes the $(n+i)$-th entry of **r** and

$$ \Upsilon_T(\sigma) = \sum_{i=1}^{n} \frac{\sin(\omega_i \sigma / 2)}{\sin(\omega_i T / 2)} r_{n+i}^2. \qquad (3.1d) $$

* The symmetric matrix $\mathbf{\Sigma}$ populated by the elements

$$ \mathbf{\Sigma}_{i,j}(t_1, \dots, t_k) = 2\mathbf{w}^\top(\mathbf{S}(\frac{T}{2}) - \mathbf{S}(-\frac{T}{2}))^{-1}\mathbf{S}(t_i - t_j + \frac{T}{2})\mathbf{r}, \quad i \le j \le k. \qquad (3.1e) $$

* $\phi_T, \psi_T \in \mathbb{R}^n$ defined by

$$ \begin{bmatrix} \phi_T(\sigma) \\ \psi_T(\sigma) \end{bmatrix} = 2\mathbf{P}^{-1}(\mathbf{S}(T/2) - \mathbf{S}(-T/2))^{-1}\mathbf{S}(\sigma/2)\mathbf{r}. \qquad (3.1f) $$

* The linear application

$$ \varphi: \mathbb{R}^k \to \mathbb{R}^{2n} \\ \lambda \mapsto 2\mathbf{P}^{-1}(\mathbf{S}(T) - \mathbf{I})^{-1} \sum_{j=1}^{k} \lambda_j \mathbf{S}(T - t_j)\mathbf{r}. \qquad (3.1g) $$
**3.2. Important theorems** Solutions of Problem (2.12) satisfy three necessary conditions NC1, NC2 and NC3. NC1 (Eq. (2.12a)) consists in finding invariant elements of a mapping defined as the product of 2k matrices, k of which are matrix exponentials. This is a challenging problem even for k and n as small as k = 2 and n = 3. Theorem 3.1 drastically simplifies it to the simple computation of the kernel of a known k × k skew-symmetric matrix. Its proof is provided in section 4.

**Theorem 3.1** $\mathbf{x}_0$ satisfies NC1 iff there exists a vector $\boldsymbol{\lambda} \in \mathbb{R}^k$ such that:

$$ \begin{cases} \Pi(t_1, \dots, t_k)\boldsymbol{\lambda} = \mathbf{0} \\ \mathbf{x}_0 = \varphi(\boldsymbol{\lambda}). \end{cases} \qquad (3.2a) $$

$$ (3.2b) $$

Theorem 3.2 shows that the condition of gap closure NC2 (Eq. (2.12b)) can be expressed in terms of the same vector $\varphi^{-1}(\mathbf{x}_0) = \boldsymbol{\lambda} \in \mathbb{R}^k$.

**Theorem 3.2** $\mathbf{x}_0$ satisfies NC2 iff there exists $\boldsymbol{\lambda} \in \mathbb{R}^k$ such that:

$$ \begin{cases} \Sigma(t_1, \dots, t_k)\boldsymbol{\lambda} = -g_0 \mathbf{j} \\ \mathbf{x}_0 = \varphi(\boldsymbol{\lambda}). \end{cases} \qquad (3.3a) $$

$$ (3.3b) $$

with $\mathbf{j} = [1, \dots, 1]^T \in \mathbb{R}^k$.

Together, theorems 3.1 and 3.2 give a practical way of finding potential solutions: it is sufficient to find $\boldsymbol{\lambda} \in \mathbb{R}^k$ verifying Eqs. (3.2a) and (3.3a) for some $(t_1, \dots, t_k)$; $\mathbf{x}_0$ can then be explicitly computed using (3.2b). The proof of Theorem 3.2 is provided in section 4.

The following theorem, which is a direct consequence of Theorems 3.1 and 3.2, addresses the existence of potential solutions and thus of nonsmooth modes of vibration. It implies that, if $g_0 \neq 0$, then nonsmooth modes are two-dimensional surfaces of the phase space in the neighbourhood of a generic point, as in the smooth case.

Let $S$ denote the set of impact times for which a potential solution exists:

$$ S := \{(t_1, \dots, t_k) \mid \exists \boldsymbol{\lambda} \in \mathbb{R}^k \text{ solution of } \Pi(t_1, \dots, t_k)\boldsymbol{\lambda} = \mathbf{0} \text{ and } \Sigma(t_1, \dots, t_k)\boldsymbol{\lambda} = -g_0 \mathbf{j}\} \quad (3.4) $$

Then, there exists a one-parameter continuum of potential solutions, for some specific $t_1, \dots, t_k$.

**Theorem 3.3** [Dimension of the space of potential solutions] A generic point of $S$ is:

* a point of a curve of $S$ if $g_0 \neq 0$;

* an isolated point of $S$ if $g_0 = 0$.

The proof of this theorem is given in section 5.

**Corollary 3.1** [Dimension of nonsmooth modes of vibration seen as manifolds] If $\mathbf{x}$ is an admissible solution for some generic impact times $s = (t_1, \dots, t_k)$, then:

* $s$ is a point of a curve in the space of the impact times, and the neighbourhood of $s$ on the curve corresponds to admissible solutions, if $g_0 \neq 0$ (manifold of dimension 2 in the phase space);

* $s$ is an isolated point if $g_0 = 0$ (manifold of dimension 1 in the phase space).

Corollary 3.1 shows that when $g_0 \neq 0$, an admissible solution is necessarily part of a one-parameter continuum of periodic trajectories, which is by definition a nonsmooth mode. This one-parameter continuum of trajectories is a two-dimensional surface in the phase space, since every trajectory is a curve. In the end, it shows that periodic solutions of the initial nonsmooth problem form, when they exist, two-dimensional manifolds, as in the linear and smooth nonlinear counterparts. Moreover, Theorem 3.4 gives the explicit expressions of the $n$ positions and $n$ velocities as functions of time. Its proof is provided in section 4.

**Theorem 3.4** [Expression of modal positions and velocities] If $\mathbf{x}_0$ is a potential solution for impact times $t_1, \dots, t_k$, then the time-evolution of the components of $\mathbf{x} = [\mathbf{u} \ \mathbf{v}]^T$ are given $\forall i \in [0, k-1]$,
$\forall \tau \in [0, \sigma_{i+1}]$ by

$$ \mathbf{u}(\tau + t_i) = \sum_{j=1}^{i} \lambda_j \phi_T(2(\tau + t_i - t_j) - T) + \sum_{j=i+1}^{k} \lambda_j \phi_T(2(\tau + t_i - t_j) + T) \quad (3.5) $$

$$ \mathbf{v}(\tau + t_i) = \sum_{j=1}^{i} \lambda_j \psi_T(2(\tau + t_i - t_j) - T) + \sum_{j=i+1}^{k} \lambda_j \psi_T(2(\tau + t_i - t_j) + T). \quad (3.6) $$

where $\lambda_j$ is the j-th component of $\lambda := \varphi^{-1}(\mathbf{x}_0)$.

An insight of these theorems is shown in Figure 2 for $g_0 \neq 0$. For a generic point $s = (t_1, \dots, t_k)$ of the space of the impact times, $\Sigma(s)$ is invertible and Theorems 3.1 and 3.2 indicate that $s$ must satisfy $\Pi(s)\Sigma^{-1}(s)\mathbf{j} = \mathbf{0}$ for a potential solution to exist. Such an $s$ determines a vector $\lambda = -g_0\Sigma^{-1}(s)\mathbf{j}$ and the time-evolutions of positions and velocities are given by Theorem 3.4 as functions of $\lambda$. They correspond to trajectories in the phase space. If the unilateral constraints are satisfied by the potential solution, the latter qualifies as admissible solution. Corollary 3.1 shows that such a point $s$ is generically a point of a curve in the space of impact times; the neighbourhood of $s$ on the curve also corresponds to admissible solutions. Along the neighbourhood of $s$ in $S$, the corresponding admissible solutions describe a $k$-ipp nonsmooth mode in the phase space.

<figure><img src="image_3.png" /><figcaption>**Figure 2:** Summary of the main results for a generic point $s \in S$ with $g_0 \neq 0$. When $s = (t_1, \dots, t_k)$ travels along the red curve $S$, the periodic solution describes a nonsmooth mode with $k$ impacts per period.</figcaption></figure>

**Remark 3.1.** As shown in the proof of Theorem 3.3, finding $\lambda$ and thus $\mathbf{x}_0$ requires solving $k-1$ equations. However, the set of equations depends on the parity of $k$. Because $\Pi$ is a skew-symmetric matrix, solving the $k-1$ equations can be significantly simplified using the Pfaffian operator, especially when $k$ is odd. Finally:

* If $k$ is odd, for almost every $s$, $\mathrm{rank}(\Pi(s)) = k-1$ and it can be shown that its kernel is spanned by the vector $\mathbf{k}$ of the Pfaffians of the diagonal $k-1 \times k-1$ submatrices of $\Pi$ with alternating signs. Hence, $\mathbf{k}$ writes as a function of the $k$ impact times. Nevertheless $\Pi\lambda = \mathbf{0}$ so there exists a $\mu \neq 0$ such that $\lambda = \mu\mathbf{k}$. Eq. (3.3a) implies that all the rows of $\Sigma(t_1, \dots, t_k)\mathbf{k}$ have identical value; this provides the $k-1$ equations. Then $\lambda = \mu\mathbf{k}$ where $\mu$ is determined by $\Sigma(t_1, \dots, t_k)\mathbf{k} = -\frac{g_0}{\mu}\mathbf{j}$ and the sought initial solution is given by $\mathbf{x}_0 = \varphi(\lambda)$.

* If $k$ is even, for almost every $s$, $\mathrm{rank}(\Pi(s)) = k$ and there is no $\lambda \neq \mathbf{0}$ solution of Eq. (3.2a). Hence one must solve $\mathrm{pf}(\Pi) = 0$ for $(t_1, \dots, t_k)$. When this equation is satisfied, for almost every $s$, $\mathrm{rank}(\Pi(s)) = k-2$, and $\Pi(t_1, \dots, t_k)\lambda = \mathbf{0}$ provides $k-2$ independent equations. Again, a total of $k-1$ equations are to be solved.

**Remark 3.2.** Theorem 3.3 illustrates the effectiveness of Theorems 3.1 and 3.2, but the latter are also constructive. For example, let us treat the problem for $k=1$. Then, $\Pi = [0]$ because it is skew-symmetric, and the only condition on $\lambda \in \mathbb{R}$ and $t_1 = T$ is $\Sigma(T)\lambda = -g_0\mathbf{j}$. This leads to the existence of a continuum of potential solutions whenever $\Sigma(T) \neq 0$, that is almost everywhere in $\mathbb{R}_+^*$. Interestingly,
this result does not depend on $n$. Theorems 3.1 and 3.2 hence allow to solve the periodic problem with one impact per period in a few lines only. This is further detailed on a simple example in subsection 7.2.

**Remark 3.3.** No general result regarding the last necessary condition NC3 are presented here. In practice, the first step is to compute potential solutions and then check NC3. However, it follows from Theorem 3.1 that a general admissible solution is an element of a continuum of admissible solutions. As a consequence, the corresponding nonsmooth mode can be numerically retrieved from one admissible solution using continuation techniques. It is important to note that potential solutions can be straightforwardly filtered out by checking that the sign of the pre-impact velocities corresponds to non-penetration.

**4. Proof of Theorems 3.1, 3.2, 3.4** As seen in Eq. (2.11), the dynamics is a combination of free flights, captured by $\tilde{S}$, and impacts, accounted for with $\tilde{N}$. Matrices $\tilde{S}$ and $\tilde{N}$ respectively correspond to two endomorphisms of fundamentally different nature: the free flight can be regarded as $n$ rotations in $n$ orthogonal planes (the linear modes) of the phase space, while, as seen from Proposition 2.1, the impact dynamics is a reflection with respect to a hyperplane and only affects the coordinates appearing explicitly in the expression of $g(u)$.

Uncoupling the dynamics via a modal basis results in an impact law which affects modal positions. Through an appropriate change of basis, this can be simplified: the dynamics is uncoupled (as with any modal basis) while the impact law does not affect the position components (which is not the case a priori). Geometrically, the new basis to be introduced is constructed in such a way that the hyperplane defining the impact law intersects the $n$ invariant subspaces of the dynamics along basis vectors. This is done by rotating the pair of vectors spanning each stable plane such that the first vector is aligned with the intersection between the stable plane and the hyperplane, see Fig. 3.

<figure><img src="image_4.png" /><figcaption>**Figure 3:** In the new basis, the hyperplane invariant by the impact law $N$ intersects the $n$ orthogonal planes invariant by the free-flight matrix $S$ along a line spanned by the position basis vectors.</figcaption></figure>

**4.1. Free flight dynamics** Let $y = V^T u$ where $V$ corresponds to the Cholesky decomposition $M = VV^T$. The free flight dynamics, Eq. (2.4), then reads:

$$ \ddot{y} + V^{-1}KV^{-T}y = 0. \qquad (4.1) $$

By symmetry and definite positiveness, $V^{-1}KV^{-T} = Q\Omega^2Q^T$ for some orthogonal matrix $Q$ and a diagonal matrix $\Omega = \text{diag}(\omega_1, \dots, \omega_n)$ Let $z_1 = \Omega z$, $z_2 = \dot{z}$ and $z = [z_1, z_2]^T$, it follows that $\dot{z}_2 = -\Omega z_1$, or:

$$ \dot{z} = \begin{bmatrix} 0 & \Omega \\ -\Omega & 0 \end{bmatrix} z. \qquad (4.2) $$

Defining the transformation matrix $P$ from the new to the old basis, $Px = z$, matrices $\tilde{S}$ and $\tilde{N}$ become in the new basis $S := P\tilde{S}P^{-1}$ and $N := P\tilde{N}P^{-1}$, with

$$ P = \begin{bmatrix} \Omega Q^T V^T & 0 \\ 0 & Q^T V^T \end{bmatrix}. \qquad (4.3) $$
After this change of basis, the solutions of the dynamics during a free phase are simply given by

$$z(t) = S(t)z_0 = \exp\left(t \begin{bmatrix} \mathbf{0} & \mathbf{\Omega} \\ -\mathbf{\Omega} & \mathbf{0} \end{bmatrix}\right)z_0 = \begin{bmatrix} \cos(\mathbf{\Omega}t) & \sin(\mathbf{\Omega}t) \\ -\sin(\mathbf{\Omega}t) & \cos(\mathbf{\Omega}t) \end{bmatrix}z_0 \quad (4.4)$$

where $\cos(\mathbf{\Omega}t)$ and $\sin(\mathbf{\Omega}t)$ are diagonal matrices obtained by mapping cos and sin to each component of $\mathbf{\Omega}$. The matrix $S(t)$ is orthogonal ($S(t)S(t)^{\top} = \mathbf{I}$). Assuming $t$ is not a multiple of a $2\pi/\omega_i$ ($\omega_i t \notin 2\pi\mathbb{Z}$ for all $i$), it can then be shown that $\mathbf{I} - S(t)$ is non-singular and

$$(\mathbf{I} - S(t))^{-1} = \frac{1}{2} \begin{bmatrix} \mathbf{I}_n & \cot \frac{\mathbf{\Omega}t}{2} \\ -\cot \frac{\mathbf{\Omega}t}{2} & \mathbf{I}_n \end{bmatrix} \quad (4.5)$$

implying the two following identities:

$$(\mathbf{I} - S(t))^{-1} + (\mathbf{I} - S(t))^{-\top} = \mathbf{I} \quad (4.6)$$

$$S(t)(\mathbf{I} - S(t))^{-1} = (S(-t) - \mathbf{I})^{-1}. \quad (4.7)$$

### 4.2. Formulation of the impact law
We now proceed with the demonstration of Proposition 2.1. Through Eq. (2.5), the normal contact velocity $\gamma = J_g(\mathbf{u})\mathbf{v} = \mathbf{w}_1^\top\mathbf{v}$ is such that

$$\gamma^{+} - \gamma^{-} = \mathbf{w}_{1}^{\top}(\mathbf{v}^{+} - \mathbf{v}^{-}) = \mathbf{w}_{1}^{\top}\mathbf{M}^{-1}\mathbf{w}_{1}\Lambda = G\Lambda. \quad (4.8)$$

The scalar quantity $G = \mathbf{w}_1^\top\mathbf{M}^{-1}\mathbf{w}_1$ is non-zero because $\mathbf{M}^{-1}$ is positive definite. Inserting (2.6) into (4.8) yields $-2\gamma^{-} = G\Lambda$ hence $\Lambda = -2G^{-1}\gamma^{-} = -2G^{-1}\mathbf{w}_1^\top\mathbf{v}^{-}$ and thus, using Eq. (2.5), $\mathbf{v}^{+} = (\mathbf{I} - 2\mathbf{M}^{-1}\mathbf{w}_1G^{-1}\mathbf{w}_1^\top)\mathbf{v}^{-}$, providing the expression for $\mathbf{N}$ introduced in Eq. (2.7). With $\mathbf{r}_2 = \mathbf{V}^{-1}\mathbf{w}_1G^{-1/2}$, it comes that $\mathbf{N}_2 = \mathbf{V}^{-\top}(\mathbf{I} - 2\mathbf{r}_2\mathbf{r}_2^\top)\mathbf{V}^\top$ which ends the proof. Vector $\mathbf{r}_2$ is a unit vector since $||\mathbf{r}_2||^2 = \mathbf{r}_2^\top\mathbf{r}_2 = G^{-1/2}\mathbf{w}_1^\top\mathbf{M}^{-1}\mathbf{w}_1G^{-1/2} = G^{-1/2}GG^{-1/2} = 1$.

Introducing $\mathbf{r} = [\mathbf{0} \quad \mathbf{r}_2^\top\mathbf{Q}]^\top \in \mathbb{R}^{2n}$, it comes that

$$\mathbf{z}^{+} = \mathbf{N}\mathbf{z}^{-} \quad \text{with} \quad \mathbf{N} = \mathbf{I} - 2\mathbf{r}\mathbf{r}^{\top} \quad (4.9)$$

which is known as a Householder matrix [14], correspond to a reflection with respect to the hyperplane $\mathbf{r}^{\perp}$. It is easily seen that $\mathbf{N}$ does not affect the first $n$ elements of $\mathbf{z}^{-}$, corresponding to displacements.

**Remark 4.1.** The impact law is conservative since $\mathbf{N}_2^\top\mathbf{M}\mathbf{N}_2 = \mathbf{M}$.

**Remark 4.2.** Impact law (2.8) is not restricted to a mass impacting a rigid obstacle. It also handles two impacting masses, for instance.

### 4.3. Sequence of impacts
In this subsection, we introduce sequences which will be used to prove Theorems 3.1 and 3.2. Define the following matrix sequences for $l \in [1, k-1]$:

$$\mathbf{U}_1 = \mathbf{N}\mathbf{S}(\sigma_1), \qquad \mathbf{U}_{l+1} = \mathbf{N}\mathbf{S}(\sigma_{l+1})\mathbf{U}_l \in \mathbb{R}^{2n,2n} \quad (4.10)$$

$$\mathbf{C}_1 = \mathbf{r}, \qquad \mathbf{C}_{l+1} = [\mathbf{S}(\sigma_{l+1})\mathbf{C}_l \quad \mathbf{C}_1] \in \mathbb{R}^{2n,l} \quad (4.11)$$

$$\mathbf{L}_1 = \mathbf{r}^{\top}\mathbf{S}(\sigma_1), \qquad \mathbf{L}_{l+1} = \begin{bmatrix} \mathbf{L}_l \\ \mathbf{r}^{\top}\mathbf{S}(\sigma_{l+1})\mathbf{U}_l \end{bmatrix} \in \mathbb{R}^{l,2n} \quad (4.12)$$

$$\mathbf{T}_1 = 1, \qquad \mathbf{T}_{l+1} = \begin{bmatrix} [2\mathbf{r}^{\top}\mathbf{S}(t_{l+1} - t_1)\mathbf{r} & \dots & 2\mathbf{r}^{\top}\mathbf{S}(t_{l+1} - t_1)\mathbf{r}] & \begin{bmatrix} \mathbf{0}_{l,1} \\ 1 \end{bmatrix} \end{bmatrix} \in \mathbb{R}^{l,l}. \quad (4.13)$$

It can be shown that for all $l \in [1, k-1]$,

$$2\mathbf{C}_l^\top\mathbf{C}_l = \mathbf{T}_l + \mathbf{T}_l^\top. \quad (4.14)$$

Assuming $\mathbf{U}_l = \mathbf{S}(t_l) - 2\mathbf{C}_l\mathbf{L}_l$ leads to

$$\mathbf{S}(t_{l+1}) - 2\mathbf{C}_{l+1}\mathbf{L}_{l+1} = \mathbf{S}(\sigma_{l+1})\mathbf{S}(t_l) - 2\mathbf{S}(\sigma_{l+1})\mathbf{C}_l\mathbf{L}_l - 2\mathbf{r}\mathbf{r}^{\top}\mathbf{S}(\sigma_{l+1})\mathbf{U}_l \quad (4.15a)$$

$$= \mathbf{S}(\sigma_{l+1})(\mathbf{S}(t_l) - 2\mathbf{C}_l\mathbf{L}_l) - 2\mathbf{r}\mathbf{r}^{\top}\mathbf{S}(\sigma_{l+1})\mathbf{U}_l \quad (4.15b)$$

$$= (\mathbf{I} - 2\mathbf{r}\mathbf{r}^{\top})\mathbf{S}(\sigma_{l+1})\mathbf{U}_l = \mathbf{N}\mathbf{S}(\sigma_{l+1})\mathbf{U}_l = \mathbf{U}_{l+1} \quad (4.15c)$$
and since $U_1 = S(t_1) - 2C_1L_1$, it comes by induction that

$$U_l = S(t_l) - 2C_l L_l. \qquad (4.16)$$

It can also be shown by induction that

$$T_l L_l = C_l^\top S(t_l) \qquad (4.17)$$

hence, since $T_l$ is non-singular:

$$U_l - I = (I - 2C_l L_l S(-t_l) - S(-t_l)) S(t_l) \qquad (4.18)$$

### 4.4. Initial conditions invariant after k successive impacts (proof of Theorem 3.1)
For $t_l$ such that $\omega_i t_l \notin 2\pi\mathbb{Z}$ for all $i$,

$$U_l - I = \Delta_l (I - A_l B_l^\top) S(t_l) \qquad (4.19)$$

with $A_l = 2(I - S(-t_l))^{-1}C_l \in \mathbb{R}^{2n,l}$, $B_l^\top = L_l S(-t_l) \in \mathbb{R}^{2n,l}$ and $\Delta_l = (I - S(-t_l))$. The unknown initial conditions $z_0$ are such that $NS(\sigma_k)N...NS(\sigma_1)z_0 = z_0$, i.e. $z_0 \in \ker(U_k - I)$. From Eq. (4.19) and since $t_k = T$,

$$z_0 \in \ker(U_k - I_{2n}) \iff \left\{ \begin{array}{l} z_0 = -S(-T)A_k \lambda \\ \lambda = B_k^\top S(T) z_0 \end{array} \right. \iff \lambda \in \ker(I_k - B_k^\top A_k) \qquad (4.20)$$

and using property (4.7), $z_0$ and $\lambda$ also satisfy

$$z_0 = 2(S(T) - I)^{-1} C_k \lambda, \quad \lambda = L_k z_0. \qquad (4.21)$$

Then, following from Eq. (4.17),

$$I - B_k^\top A_k = I - 2L_k S(-T) (I - S(-T))^{-1} C_k = T^{-1} \Pi \qquad (4.22)$$

with $\Pi = T - 2C_k^\top \Delta_k C_k$, which means that finding initial conditions which are invariant by a succession of k impacts and k free flights (i.e. in $\ker(U_k - I)$) reduces to computing $\ker \Pi$. Matrix $\Pi$ is skew-symmetric because

$$\Pi + \Pi^\top = T + T^\top - 2C_k^\top \Delta_k^{-1} C_k - 2C_k^\top \Delta_k^{-\top} C_k \qquad (4.23)$$

$$= 2C_k^\top C_k - 2C_k^\top (\Delta_k^{-1} + \Delta_k^{-\top}) C_k \qquad \text{using Eq. (4.14)} \qquad (4.24)$$

$$= 0 \qquad \text{from Eq. (4.6).} \qquad (4.25)$$

Since $T$ is lower-triangular, because the jth column of $C_k$ is $S(t_k - t_j)r$ and since $S(t)$ and $\Delta_k$ commute, it comes that for $i < j$:

$$\Pi_{i,j} = -2r^\top (I - S(T))^{-1} S(t_i - t_j) r \qquad (4.26)$$

from which we can show that $\Pi_{i,j} = \Upsilon_T(2(t_j - t_i) - T)$, $\forall i < j$ with

$$\Upsilon_T(\sigma) = \sum_{i=1}^{n} \frac{\sin(\omega_i \sigma / 2)}{\sin(\omega_i T / 2)} r_{n+i}^2. \qquad (4.27)$$

To summarize, Eq. (2.12a) holds iff $U_k z_0 = z_0$ iff $\Pi \lambda = 0$ where $\Pi$ is a known skew-symmetric matrix function of the impact times $t_1, \dots, t_k$ with $t_k = T \notin \{0, 2\pi/\omega_1, \dots, 2\pi/\omega_n\}$. From Eq. (4.21), finding solutions $x_0 \in \mathbb{R}^{2n}$ reduces to first finding $\lambda \in \mathbb{R}^k$ and $x_0$ can then be recovered using

$$x_0 = 2P^{-1}S(-T)(I - S(-T))^{-1}C_k \lambda = 2P^{-1}(S(T) - I)^{-1} \sum_{j=1}^{k} \lambda_j S(T - t_j)r \qquad (4.28)$$

which ends the proof of Theorem 3.1.

**Remark 4.3.** The physical meaning of the condition $T \notin \{0, 2\pi/\omega_1, \dots, 2\pi/\omega_n\}$ is that $T$ cannot be the period of a linear mode. This condition is not really restrictive and does not prevent the investigation of the resonant behaviour.

**Remark 4.4.** The quantity $2(t_j - t_i) - T$ for $i < j$ is the signed sum of $\sigma_1, \dots, \sigma_k$ where the $j - i$ variables $\sigma_q$ for $i < q \le j$ have the sign - while the $k - j + i$ other variables have the sign +.
**4.5. Gap closure at impact times (proof of Theorem 3.2)** We now proceed with the proof of Theorem 3.2, which ensures gap closure at impact times in terms of $\lambda$ (NC2). The gap closure conditions at $t_1, \dots, t_k$ write $\tilde{g}(\mathbf{x}(t_l)) = 0 = \mathbf{w}_1^\top \mathbf{u}(t_l) + g_0$ (Eq. (2.12b)) or, with $\mathbf{w} = \mathbf{P}^{-\top} [\mathbf{w}_1^\top \quad \mathbf{0}]^\top$, $\forall l \in [1, k]$, $\mathbf{w}^\top \mathbf{U}_l \mathbf{z} = -g_0$. (4.29)

Using Eqs. (4.16) and (4.21), it comes that

$$ \mathbf{w}^{\top}(\mathbf{S}(t_l)\mathbf{z} - 2\mathbf{C}_l\mathbf{L}_l\mathbf{z}) = \mathbf{w}^{\top}(\mathbf{S}(t_l)(\mathbf{S}(t_k) - \mathbf{I})^{-1}\mathbf{C}_k - [\mathbf{C}_l \quad \mathbf{0}])\mathbf{\lambda} = -g_0 \quad (4.30) $$

yielding $k$ equations which can be gathered in a matrix equation

$$ \Sigma \lambda = -g_0 \mathbf{j} \quad (4.31) $$

with

$$ \Sigma_{i,j} = \begin{cases} 2\mathbf{w}^\top(\mathbf{S}(T/2) - \mathbf{S}(-T/2))^{-1}\mathbf{S}(t_i - t_j + T/2)\mathbf{r} & \text{if } i < j \\ 2\mathbf{w}^\top(\mathbf{S}(T/2) - \mathbf{S}(-T/2))^{-1}\mathbf{S}(t_i - t_j - T/2)\mathbf{r} & \text{if } i \le j \end{cases} \quad (4.32) $$

$$ \mathbf{w}^\top(\mathbf{S}(T/2) - \mathbf{S}(-T/2))^{-1}(\mathbf{S}(\sigma) - \mathbf{S}(-\sigma))\mathbf{r} = 0, \text{ implying that } \Sigma \text{ a symmetric matrix.} \quad (4.33) $$
and $\mathbf{j}^\top = [1 \quad \dots \quad 1]$. For any $\sigma$, $(\mathbf{S}(T/2) - \mathbf{S}(-T/2))^{-1}(\mathbf{S}(\sigma) - \mathbf{S}(-\sigma))$ is a diagonal matrix so that

To summarize, the initial conditions $\mathbf{z}$ satisfy the gap closure conditions (NC2) iff $\Sigma \lambda = -g_0 \mathbf{j}$ with $\lambda = \mathbf{L}_k \mathbf{z}$. Theorem 3.2 has been proven and $\Sigma$ was shown to be symmetric.

**4.6. Explicit expression of the positions and velocities (proof of Theorem 3.4)** The expressions of the positions $\mathbf{u}$ and velocities $\mathbf{v}$ as functions of the impact times $t_1, \dots, t_k$ are now derived. Given $k$ and $i \le k$, $\forall \tau \in [0, \sigma_i]$:

$$ \mathbf{z}(\tau + t_i) = \mathbf{S}(\tau)\mathbf{U}_i\mathbf{z}_0 \quad (4.34) $$

$$ = \mathbf{S}(\tau + t_i)\mathbf{z}_0 - 2\mathbf{S}(\tau)\mathbf{C}_l\mathbf{L}_l\mathbf{z}_0 \quad \text{from (4.16)} $$

$$ = 2 \sum_{j=1}^{k} \lambda_j \mathbf{S}(\tau + t_i - t_j)(\mathbf{I} - \mathbf{S}(-T))^{-1}\mathbf{r} - 2 \sum_{j=1}^{i} \lambda_j \mathbf{S}(\tau + t_i - t_j)\mathbf{r} \quad \text{from (4.7), (4.21)} $$

$$ = 2 \sum_{j=i+1}^{k} \lambda_j (\mathbf{S}(T/2) - \mathbf{S}(-T/2))^{-1} \mathbf{S}(\tau + t_i - t_j) \mathbf{r} \\ \qquad + 2 \sum_{j=1}^{i} \lambda_j (\mathbf{S}(T/2) - \mathbf{S}(-T/2))^{-1} \mathbf{S}(\tau + t_i - t_j + T/2) \mathbf{r} \quad \text{using (4.6)} $$

Setting

$$ \phi_T(\sigma) = 2[\mathbf{V}^{-\top}\mathbf{Q}\mathbf{\Omega}^{-1} \quad \mathbf{0}] (\mathbf{S}(T/2) - \mathbf{S}(-T/2))^{-1} \mathbf{S}(\sigma/2)\mathbf{r}, \quad (4.35) $$

$$ \psi_T(\sigma) = 2[\mathbf{0} \quad \mathbf{V}^{-\top}\mathbf{Q}] (\mathbf{S}(T/2) - \mathbf{S}(-T/2))^{-1} \mathbf{S}(\sigma/2)\mathbf{r} = 2 \frac{d\phi_{p,T}(t)}{dt}, \quad (4.36) $$

Eq. (4.34) can be recast into

$$ \mathbf{u}(\tau + t_i) = \sum_{j=1}^{i} \lambda_j \phi_T(2(\tau + t_i - t_j) - T) + \sum_{j=i+1}^{k} \lambda_j \phi_T(2(\tau + t_i - t_j) + T), \quad (4.37) $$

$$ \mathbf{v}(\tau + t_i) = \sum_{j=1}^{i} \lambda_j \psi_T(2(\tau + t_i - t_j) - T) + \sum_{j=i+1}^{k} \lambda_j \psi_T(2(\tau + t_i - t_j) + T), \quad (4.38) $$

which ends the proof of Theorem 3.4.

The main notations are summarized in the nomenclature at the end.
**5. Dimension of the solution space and nonsmooth modes [proof of Theorem 3.3]** We now prove Theorem 3.3 which justifies the notion of nonsmooth modes as two-dimensional surfaces in the phase space. It implies the existence of continua of potential solutions.

*Proof.* Assume $g_0 \neq 0$. Given a generic point $s \in S$, $\Sigma(s)$ is invertible and then $\lambda = -g_0\Sigma^{-1}(s)\mathbf{j}$. So in the neighbourhood of $s$, $S$ is the set of $(t_1, \dots, t_k)$ such that $\Pi(t_1, \dots, t_k)\Sigma^{-1}(t_1, \dots, t_k)\mathbf{j} = \mathbf{0}$. This gives a priori $k$ equations $\sum_{j=1}^k \Pi_{i,j}\lambda_j = 0$ for $i \in [1, k]$. But $\lambda \neq \mathbf{0}$ so there is an index $i_0$ such that $\lambda_{i_0} \neq 0$. And $\Pi$ is skew-symmetric so

$$0 = \lambda^\top \Pi \lambda = \sum_{i=1}^{k} \lambda_i \sum_{j=1}^{k} \Pi_{i,j} \lambda_j = \lambda_{i_0} \sum_{j=1}^{k} \Pi_{i_0,j} \lambda_j + \sum_{i \neq i_0} \lambda_i \sum_{j=1}^{k} \Pi_{i,j} \lambda_j. \quad (5.1)$$

Hence, the nullity of the $k-1$ rows $i \neq i_0$ of $\Pi\Sigma^{-1}\mathbf{j}$ implies the nullity of row $i_0$, which proves that, in the neighbourhood of $s$, $S$ is the set of solutions of $k-1$ equations in the $k$ unknowns $t_1, \dots, t_k$. These $k-1$ equations can be linearized around $s$ which yields $k-1$ linear equations $E_i = 0$ with $i \neq i_0$ in $t_1, \dots, t_k$. In the following, we assume w.l.o.g. that $i_0 = k$. The Jacobian matrix of $(t_1, \dots, t_k) \mapsto (E_1, \dots, E_{k-1})$ is generically of rank $k$, because otherwise the $k$ submatrices of dimension $k-1 \times k-1$ should be singular, corresponding to $k$ additional equations. Hence, generically, there exists a linear form $F$ such that $(E_1, \dots, E_{k-1}, F)$ is of rank $k$, hence the Jacobian matrix of $(t_1, \dots, t_k) \mapsto (E_1, \dots, E_{k-1}, F)$ is invertible. The inverse function theorem implies that $(E_1, \dots, E_{k-1}, F)$ yields locally a system of curvilinear coordinates and $S$ is the set of points whose first $k-1$ coordinates are zero, corresponding to curves. This ends the proof for the case $g_0 \neq 0$.

If $g_0 = 0$, $\Sigma\lambda = \mathbf{0}$ with $\lambda \neq \mathbf{0}$ so necessarily $\Sigma$ is singular. This gives one equation $\det(\Sigma(t_1, \dots, t_k) = 0$ and when this condition is satisfied, $\Sigma$ is generically of rank $k-1$. Its adjunct matrix is of rank 1 and the kernel of $\Sigma$ is spanned by one column $c$ of the transpose of its adjunct matrix. The other condition to satisfy is $\Pi c = \mathbf{0}$ which generically yields, as previously, $k-1$ independent equations. A total of $k$ equations to be solved in $k$ unknowns is available. There is generically a countable set of solutions. ■

The sub-generic case where the sequence $(\sigma_i)$ is invariant by a transformation $n \mapsto n-a$ for some $a \in \mathbb{N}$ yields $k-1$ linearly dependent equations $E_i = 0$, $i = 1, \dots, k-1$. However, from Proposition 6.6, the solution space is still of dimension 1, provided $g_0 \neq 0$.

We now prove Corollary 3.1 concerning admissible solutions, which is a direct consequence of Theorem 3.3. It proves the existence of nonsmooth modes as soon as an admissible solution exists, in the generic case.

*Proof.* Admissible solutions are potential solutions satisfying NC3 so if an admissible solutions has impact times $s \in \mathbb{R}^k$, then $s \in S$. If $g_0 \neq 0$, Theorem 3.3 implies that a generic $s$ lies on a curve of impact times corresponding to some potential solutions. Since NC3 is satisfied for $x$, it is satisfied in the neighbourhood of $s$ by continuity of $t \mapsto \tilde{g}(x(t))$ hence the curve of potential solutions is also locally a curve of admissible solutions. If $g_0 = 0$, from Theorem 3.3 $S$ forms, generically, a set of isolated points, hence $s$ is an isolated point of $\mathbb{R}^k$. ■

**6. Symmetries in the modal orbits** Theorems 3.1 and 3.2 state that for a given $k$, potential solutions $x_0 \in \mathbb{R}^{2n}$ can be sought as the images by $\varphi$ of vectors $\lambda \in \mathbb{R}^k$ solutions of

$$\{ \Pi(t_1, \dots, t_k)\lambda = \mathbf{0} \} \qquad (6.1a)$$

$$\{ \Sigma(t_1, \dots, t_k)\lambda = -g_0 \mathbf{j}. \} \qquad (6.1b)$$

When $g_0 \neq 0$, non-trivial solutions may exist only if $\Sigma(t_1, \dots, t_k)$ is invertible and then the two vectorial conditions reduce to $\Pi(t_1, \dots, t_k)\Sigma^{-1}(t_1, \dots, t_k)\mathbf{j} = \mathbf{0}$ which yields a priori $k-1$ equations in the $k$ unknowns $t_1, \dots, t_k$, and then a potential solution for $t_1, \dots, t_k$ is given by $x_0 = \varphi(\lambda) = (-g_0\Sigma^{-1}\mathbf{j})$.

However, some potential solutions $x_0$ correspond to free-flight durations $\sigma_i = t_i - t_{i-1}$ which are not all independent. For example, Fig. 4 illustrates the position of a contacting node of a 7-ipp solution for which the sequence of free-flight durations $\dots, \sigma_1, \sigma_2, \sigma_3, \sigma_4, \sigma_5, \sigma_6, \sigma_7, \dots$ is of the form $\dots, c, b, a, b, c, d, \dots$. For more detailed about Fig. 4, refer to subsection 7.4 (Fig. 4 is the same as Fig. 15③). The sequence $(\sigma_i)$ is invariant, modulo 7, by the transformation $n \mapsto 6-n$, as shown in Tab. 1. The position features two axes of symmetry corresponding to the symmetry in the sequence $(\sigma_i)$.

The current section shows in general how the symmetry in the sequence is reflected by symmetries in the motion.
$$
\begin{array}{c|ccccccc}
n & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\
\sigma_n & c & b & a & b & c & d & d \\
6-n & 5 & 4 & 3 & 2 & 1 & 7 & 6 \\
\sigma_{6-n} & c & b & a & b & c & d & d
\end{array}
$$

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Example of motions with symmetries of the free-flight durations sequence.</figcaption></figure>

**6.1. Symmetry of a k-periodic sequence** Let $(\sigma_n)_{n \in \mathbb{Z}}$ denote a sequence of period $k \in \mathbb{N}^*$.

**Proposition 6.1** *If $(\sigma_n)$ is invariant by the transformation $n \mapsto a-n$ for some $a \in \mathbb{Z}$, then $a$ is unique modulo $k$.*

*Proof.* The sequence $(\sigma_n)$, which is invariant by the transformation $n \mapsto a-n$, is also invariant by $n \mapsto b-n$ iff it is invariant by $n \mapsto b-a+n$, i.e. iff $b-a$ is a multiple of $k$. Hence if such an $a$ exists, it is unique modulo $k$. ■

The transformation $n \mapsto a - n$ is involutory so its cycle decomposition only involves transpositions. The following proposition lists the number of invariant elements of $n \mapsto a - n$, which explains several properties proven later.

**Proposition 6.2 [Number of invariant elements]** *The equation $a - i = i$ of unknown $i$ has:*

* exactly one solution modulo $k$ if $k$ is odd;

* no solution if $k$ is even and $a$ is odd;

* exactly two solutions modulo $k$ if both $k$ and $a$ are even.

*Proof.* If $k = 2r-1$, in $\mathbb{Z}/k\mathbb{Z}$, $2r=1$ implies that, the equivalence class [2] is invertible, of inverse $[r]$. Hence $a = 2i$ iff $ra = 2ri = i$ which gives exactly one solution modulo $k$. If $k = 2r$, $a = 2i$ has no solution if $a$ is odd. If $a$ is even, $i = a/2$ modulo $r$ so there are exactly two solutions: $a/2$ and $r+a/2$. ■

**6.2. Extension of the impact times sequence to $\mathbb{Z}$** In order to properly explore symmetries, the sequence $(t_i)$ is extended from $i \in \mathbb{N}$ to $i \in \mathbb{Z}$. The sequence of impact times $(t_i)_{i \in \mathbb{N}}$ introduced in subsection 2.2 is such that $t_0 = 0$ and $\forall i \in \mathbb{N}$, $t_{i+1} = t_i + \sigma_{i+1}$. The sequence of free-flight durations $(\sigma_i)_n$ can be defined over $\mathbb{Z}$ by periodicity, so the definition of $(t_i)_n$ can also be extended to $\mathbb{Z}$ with $t_0 = 0$ and $\forall i \in \mathbb{Z}$, $t_{i+1} = t_i + \sigma_{i+1}$.

It follows that the definition of $\mathbf{\Pi}$ and $\mathbf{\Sigma}$ can be extended too. Recall that the skew-symmetric matrix $\mathbf{\Pi}$ is defined for $0 < i < j \le k$ by $\Pi_{i,j} = \Upsilon_T(2(t_j - t_i) - \sigma)$. Introducing the sign function $\delta$ defined on $[-(k+1), k+1]$ with the convention $\delta(0) = 0$, this can be rewritten as

$$ \forall (i, j) \in [1, k]^2, \quad \Pi_{i,j} = \Upsilon_T(2(t_j - t_i) - \delta(j - i)T). \qquad (6.2) $$

Similarly, the symmetric matrix $\mathbf{\Sigma}$ is defined for $i < j \le k$ by

$$ \Sigma_{i,j} = 2\mathbf{w}^\top(\mathbf{S}(\frac{T}{2}) - \mathbf{S}(-\frac{T}{2}))^{-1}\mathbf{S}(t_i - t_j + \frac{T}{2})\mathbf{r} := \Gamma_T(2(t_j - t_i) - T). \qquad (6.3) $$
Introducing the sign function $\epsilon$ on $\mathbb{R}^+(k+1), k+1\mathbb{R}$ with the convention $\epsilon(0) = 1$, it can be rewritten as

$$ \forall (i, j) \in \mathbb{R}^+(k+1), k+1\mathbb{R}^2, \quad \Sigma_{ij} = \Gamma_T(2(t_j - t_i) - \epsilon(j - i)T) \qquad (6.4) $$

because $\Gamma_T(\cdot)$ is an even function. Finally, Definitions (6.2) and (6.4) can be extended from $\mathbb{R}^+(k+1)\mathbb{R}^2$ to $\mathbb{Z}^2$ by defining $\delta$ and $\epsilon$ on $\mathbb{Z}$ as:

$$ \forall p \in \mathbb{Z}, \quad \delta(p+k) = 2 + \delta(p) \quad \text{and} \quad \epsilon(p+k) = 2 + \epsilon(p) \qquad (6.5) $$

Then $(\Pi_{ij})$ and $(\Sigma_{ij})$ define two sequences which are $k$-periodic with respect to both indices.

**6.3. Symmetries of positions in the time domain** Let $(\sigma_i)_{i \in \mathbb{Z}}$ denote a $k$-periodic sequence invariant by $n \mapsto a-n$ for some $a \in \mathbb{Z}$. From the previous subsection, this sequence uniquely defines the sequence $(t_i)_{i \in \mathbb{Z}}$.

**Proposition 6.3** The sequence $(t_i + t_{a-i-1})_i$ is of constant value $t_{a-1}$.

*Proof.* $\forall i < j$,

$$ t_j - t_i = \sum_{p=i+1}^{j} \sigma_p = \sum_{p=i+1}^{j} \sigma_{a-p} = \sum_{p=a-j}^{a-i-1} \sigma_p = t_{a-i-1} - t_{a-j-1} \qquad (6.6) $$

so that $(t_i + t_{a-i-1})_i$ is constant of value $t_0 + t_{a-1} = t_{a-1}$.

<figure><img src="image_6.png" /></figure>

We can verify Proposition 6.3 with the example in Fig. 4, corresponding to $a = 6$:

$$ t_2 + t_{6-2-1} = t_2 + t_3 = (c+b) + (c+b+a) = c+b+a+b+c = t_5 \qquad (6.7a) $$

$$ t_4 + t_{6-4-1} = t_4 + t_1 = (c+b+a+b)+(c) = c+b+a+b+c = t_5. \qquad (6.7b) $$

**Theorem 6.4** Let $\mathbf{u}$ denote the position defined over $\mathbb{R}$ corresponding to free dynamics of durations $(\sigma_i)_{i \in \mathbb{Z}}$. If $g_0 \neq 0$ and $\Sigma$ is invertible, then

$$ \forall t \in \mathbb{R}, \quad \mathbf{u}(t) = \mathbf{u}(t_{a-1} - t). \qquad (6.8) $$

*Proof.* If $g_0 \neq 0$, it follows from Theorems 3.1 and 3.2 that for a given $(t_1, \dots, t_k)$, $u_0$ is uniquely determined by $x_0\varphi(\lambda) = \varphi(-g_0\Sigma^{-1}j)$. Then $t \mapsto \mathbf{u}(t)$ and $t \mapsto \mathbf{u}(t_{a-1} - t)$, which are both solutions of $\mathbf{M}\ddot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{0}$, correspond to the same $g_0$ and the same free dynamics durations: $\sigma_1 = \sigma_{a-1}$, $\sigma_2 = \sigma_{a-2}$, etc., therefore $\mathbf{u}(t) = \mathbf{u}(t_{a-1} - t)$ for all $t \in \mathbb{R}$.

<figure><img src="image_7.png" /></figure>

**Corollary 6.1** $\forall t \in \mathbb{R}$, let $\mathbf{u}(t)$ denote the positions corresponding to the free dynamics of durations $(\sigma_i)_{i \in \mathbb{Z}}$ with $g_0 \neq 0$. Then $\mathbf{u}$ has the following axes of symmetry:

$$ 2t = t_{a-1+pk}, \quad p \in \mathbb{Z}. \qquad (6.9) $$

These correspond to two axes of symmetry on every period.

*Proof.* From Theorem 6.4

$$ \mathbf{u}\left(t + \frac{t_{a-1+pk}}{2}\right) = \mathbf{u}\left(t + \frac{t_{a-1} + pT}{2}\right) = \mathbf{u}\left(t_{a-1} - \left(t + \frac{t_{a-1} + pT}{2}\right)\right) \qquad (6.10a) $$

$$ = \mathbf{u}\left(\frac{t_{a-1} - pT}{2} - t\right) = \mathbf{u}\left(\frac{t_{a-1} + pT}{2} - t\right) \qquad (6.10b) $$

where the $T$-periodicity of $\mathbf{u}$ is used in the last equality. This proves that $t_{a-1+pk}/2$ is an axis of symmetry of $\mathbf{u}$, for any $p \in \mathbb{Z}$. This corresponds to two axes of symmetry in $[0, T)$ because if $t_{a-1+pk}/2 \in [0, T)$ for some $p$, then exactly one of $t_{a-1+pk}/2 + T/2$ and $t_{a-1+pk}/2 - T/2$ is in $[0, T)$ too.

<figure><img src="image_8.png" /></figure>

**Corollary 6.2** Under the assumptions of Theorem 6.4, the following lines are axes of symmetry for the position $\mathbf{u}$, $\forall p \in \mathbb{Z}$:

* if $k = 2r - 1$ and $a = 2b$,

$$ t = \frac{t_{b-1} + t_b}{2} + pT \quad \text{and} \quad t_{b+r-1} + pT \qquad (6.11) $$
• if $k = 2r - 1$ and $a = 2b + 1$,

$$t = \frac{t_{b+r-1} + t_{b+r}}{2} + pT \quad \text{and} \quad t_b + pT \qquad (6.12)$$

• if $k = 2r$ and $a = 2b$,

$$t = \frac{t_{b-1} + t_b}{2} + pT \quad \text{and} \quad \frac{t_{b+r-1} + t_{b+r}}{2} + pT \qquad (6.13)$$

• if $k = 2r$ and $a = 2b$,

$$t = t_b + pT \quad \text{and} \quad t_{b+r} + pT \qquad (6.14)$$

*Proof.* If $k = 2r - 1$ and $a = 2b$, Proposition 6.3 with $i = b$ in Corollary 6.1 implies that

$$t = \frac{t_{a-1}}{2} = \frac{t_{2b-b-1} + t_b}{2} = \frac{t_{b-1} + t_b}{2} \qquad (6.15)$$

is an axis of symmetry. Another one is given by

$$t = \frac{t_{a-1}}{2} + \frac{T}{2} = \frac{t_{a-1}}{2} + \frac{\sigma_{b-r+1} + \cdots + \sigma_{b-1} + \sigma_b + \sigma_{b+1} + \cdots + \sigma_{b+r-1}}{2} \qquad (6.16a)$$

$$= \frac{t_{a-1}}{2} + \frac{\sigma_b + 2(\sigma_{b+1} + \cdots + \sigma_{b+r-1})}{2} = t_b + \sigma_{b+1} + \cdots + \sigma_{b+r-1} = t_{b+r-1} \qquad (6.16b)$$

• If $k = 2r - 1$ and $a = 2b + 1$, then $a + k$ is even and $(k + a)/2 = r + b$ which yields the axis of symmetry.

• If $k = 2r$ and $a = 2b$, as in Eq. (6.15), $t = (t_{b-1} + t_b)/2$ is an axis of symmetry. The other axis over the period is given by

$$\begin{aligned} t &= \frac{t_{a-1}}{2} + \frac{T}{2} = \frac{t_{a-1}}{2} + \frac{(\sigma_{b-r+1} + \cdots + \sigma_{b-1}) + \sigma_b + (\sigma_{b+1} + \cdots + \sigma_{b+r-1}) + \sigma_{b+r}}{2} \\ &= \frac{t_b + t_{b-1} + \sigma_b + 2(\sigma_{b+1} + \cdots + \sigma_{b+r-1}) + \sigma_{b+r}}{2} \\ &= t_{b+r-1} + \frac{\sigma_{b+r}}{2} = \frac{t_{b+r-1} - t_{b+r}}{2} \end{aligned} \qquad (6.17)$$

• If $k = 2r$ and $a = 2b - 1$,

$$t = \frac{t_{a-1}}{2} = \frac{t_{2b-b} + t_b}{2} = t_b \qquad (6.18)$$

which provides an axis of symmetry. And $a' = k - a = 2r - 2b - 1 = 2(r - b) - 1$ is odd too, so $t = t_{r-b} = t_{r-b+k} = t_{b+r}$ is also an axis of symmetry.

<figure><img src="image_9.png" /></figure>

The conclusions of Corollary 6.2 apply to the introductory example in Fig. 4 where $k = 7 = 2r - 1$ with $r = 4$ and $a = 6 = 2b$ with $b = 3$: there are indeed two axes of symmetry, of equations $t = (t_{b-1} + t_b)/2 = (t_2 + t_3)/2 = c + b + a/2$ and $t = t_{b+r-1} = t_6 = c + b + a + b + c + d$.

## 6.4. Dimension of the solution space in the symmetric case

The following proposition shows the consequences of symmetry of the sequence $(\sigma_i)$ on the matrices $\mathbf{\Pi}$ and $\mathbf{\Sigma}$.

**Proposition 6.5** If $(\sigma_i)$ is a symmetric sequence (invariant by $n \mapsto a - n$), then

$$\Pi_{i,j} = \Pi_{\tau(j),\tau(i)} \quad \text{and} \quad \Sigma_{i,j} = \Sigma_{\tau(j),\tau(i)} \qquad (6.19)$$

where $\tau: p \mapsto a - p - 1$ (modulo $k$).

*Proof.* Using Proposition 6.3 in Eq. (6.2):

$$\Pi_{i,j} = \Upsilon_T(2(t_j - t_i) - \delta(j - i)T) \qquad (6.20a)$$

$$= \Upsilon_T(2(t_{a-i-1} - t_{a-j-1}) - \delta((a-i-1) - (a-j-1))T) \qquad (6.20b)$$

$$= \Pi_{a-j-1,a-i-1} \qquad (6.20c)$$

In the exact same way, $\Sigma_{i,j} = \Sigma_{a-j-1,a-i-1}$.

<figure><img src="image_10.png" /></figure>

When $(\sigma_i)$ is a $k$-periodic symmetric sequence, the unknowns $t_1, \dots, t_k$ are no longer independent. However, this reduction of independent unknowns comes together with the reduction of the number of independent equations in Theorem 3.3 in such a way that the dimension of the solution space is unchanged. This result is shown in Appendix A and provided as the following Proposition.

**Proposition 6.6** If the $k$-periodic sequence $(\sigma_i)_{i \in \mathbb{Z}}$ is invariant by the transformation $i \mapsto a - i - 1$ and if $g_0 \neq 0$, the solution space of $\mathbf{\Pi}(t_1, \dots, t_k)\mathbf{\Sigma}^{-1}(t_1, \dots, t_k)\mathbf{j} = \mathbf{0}$ is of dimension 1.
## 7. Illustration on a simple serial oscillator
The mathematical results exposed in Section 3 are now illustrated on a simple oscillator: a one-dimensional spring–mass system with $n$ masses $m_1, \dots, m_n$ and $n$ springs of positive stiffnesses $k_1, \dots, k_n$ as depicted in Figure 5. For this system, the derivations slightly simplify. However, the methodology to calculate nonsmooth modes and the computational cost would be similar for more sophisticated systems.

### Supplementary material

Animations corresponding to the solutions presented in this section can be viewed online using the hyperplink provided in [39]; animations for other systems are also available there.

The positions of the masses are denoted $x_1, \dots, x_n$ and by choice, $x_1 = \dots = x_n = 0$ at the resting position. The motion of the last mass (in red) is constrained by an obstacle. The gap between the

<figure><img src="image_11.png" /><figcaption>Figure 5: Base model for the study of nonsmooth modes. Top: at rest. Bottom: in motion.</figcaption></figure>

contacting mass and the obstacle is given by $\tilde{g}(\mathbf{x}) = d - x_n$, corresponding to $g_0 = d$ where $d$ can be positive or negative. The perfectly elastic Newton impact law for the description of contact reads

$$d - x_n(t_i) = 0 \quad \Rightarrow \quad \dot{x}_n^+(t_i) = -\dot{x}_n^-(t_i). \qquad (7.1)$$

All illustrations are provided with $n = 5$, $m_i = 1/n$, $k_i = n$ and $g_0 = d = 1$.

### 7.1. Simplified governing equations

#### 7.1.1. Mass and stiffness matrices
The equilibrium equation for each mass writes, with the convention $x_0 = x_{n+1} = 0$ and $k_{n+1} = 0$:

$$\forall i \in [1, n], \quad m_i \ddot{x}_i = k_{i+1}(x_{i+1} - x_i) - k_i(x_i - x_{i-1}) \qquad (7.2)$$

or in a matrix form

$$\mathbf{M}\ddot{\mathbf{u}} = -\mathbf{K}\mathbf{u} \qquad (7.3)$$

with

$$\mathbf{u} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_{n-1} \\ x_n \end{bmatrix}, \quad \mathbf{M} = \begin{bmatrix} m_1 & & & & \\ & m_2 & & & \\ & & \ddots & & \\ & & & m_{n-1} & \\ & & & & m_n \end{bmatrix}, \quad \mathbf{K} = \begin{bmatrix} k_1 + k_2 & -k_2 & & & \\ -k_2 & k_2 + k_3 & -k_3 & & \\ & & \ddots & & \\ & & & -k_{n-1} & k_{n-1} + k_n & -k_n \\ & & & -k_n & k_n & -k_n \end{bmatrix}$$

which provides the **M** and **K** matrices of Problem (2.12). It is noteworthy that **M** is positive-definite as a diagonal matrix of coefficients $m_1, \dots, m_n > 0$ and so is **K**: $k_1, \dots, k_n > 0$ and the associated quadratic form $\sum_{i=1}^n k_i (x_i - x_{i-1})^2$ is positive, and zero iff $x_0 = \dots = x_n = 0$.

#### 7.1.2. Gap and impact law
The gradient $\mathbf{w}_1$ of the gap $g(\mathbf{u}) = d - x_n$ is given by $\mathbf{w}_1 = [0 \dots 0 -1]^\top$, so $G = m_n^{-1}$, $\mathbf{r}_2 = \mathbf{V}^{-1}\mathbf{w}_1\mathbf{G}^{-1/2} = -\mathbf{w}_1$ and

$$\mathbf{r} = \begin{bmatrix} \mathbf{0} \\ \mathbf{Q}^\top \mathbf{r}_2 \end{bmatrix} = \begin{bmatrix} 0 \\ \vdots \\ 0 \\ Q_{n,1} \\ \vdots \\ Q_{n,n} \end{bmatrix}, \quad \mathbf{w} = \mathbf{P}^{-\top} \begin{bmatrix} \mathbf{w}_1 \\ \mathbf{0} \end{bmatrix} = \frac{1}{m_n} \begin{bmatrix} Q_{n,1}/\omega_1 \\ \vdots \\ Q_{n,n}/\omega_n \\ 0 \\ \vdots \\ 0 \end{bmatrix}. \qquad (7.4)$$
**7.1.3. Positions and velocities** Following from Theorem 3.4, the positions and velocities **u** and **v** are given
$\forall \tau \in [0, \sigma_i+1]$ by

$$
\mathbf{u}(\tau + t_i) = \sum_{j=1}^{i} \lambda_j \phi_{p,T}(2(\tau + t_i - t_j) - T) + \sum_{j=i+1}^{k} \lambda_j \phi_{p,T}(2(\tau + t_i - t_j) + T), \quad (7.5a)
$$

$$
\mathbf{v}(\tau + t_i) = \sum_{j=1}^{i} \lambda_j \psi_{p,T}(2(\tau + t_i - t_j) - T) + \sum_{j=i+1}^{k} \lambda_j \psi_{p,T}(2(\tau + t_i - t_j) + T) \quad (7.5b)
$$

where $\phi_T(\sigma)$ and $\psi_T(\sigma)$ can be computed with

$$
\begin{bmatrix} \phi_T(\sigma) \\ \psi_T(\sigma) \end{bmatrix} = 2\mathbf{P}^{-1} \left(\mathbf{S}(T/2) - \mathbf{S}(-T/2)\right)^{-1} \mathbf{S}(\sigma/2) \mathbf{r} \quad (7.6)
$$

$$
= \begin{bmatrix} \mathbf{V}^{-\mathrm{T}} \mathbf{Q} \mathbf{\Omega}^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{V}^{-\mathrm{T}} \mathbf{Q} \end{bmatrix} \begin{bmatrix} \frac{\sin(\mathbf{\Omega} \sigma/2)}{\sin(\mathbf{\Omega} T/2)} & -\frac{\cos(\mathbf{\Omega} \sigma/2)}{\sin(\mathbf{\Omega} T/2)} \\ \frac{\cos(\mathbf{\Omega} \sigma/2)}{\sin(\mathbf{\Omega} T/2)} & \frac{\sin(\mathbf{\Omega} \sigma/2)}{\sin(\mathbf{\Omega} T/2)} \end{bmatrix} \begin{bmatrix} \mathbf{0} \\ Q_{n,1} \\ \vdots \\ Q_{n,n} \end{bmatrix}. \quad (7.7)
$$

It follows that

$$
\phi_{p,T}(\sigma) = - \frac{1}{\sqrt{m_p}} \sum_{j=1}^{n} \frac{\cos(\omega_j \sigma/2)}{\omega_j \sin(\omega_j T/2)} Q_{p,j} Q_{n,j}, \quad (7.8)
$$

$$
\psi_{p,T}(\sigma) = \frac{1}{\sqrt{m_p}} \sum_{j=1}^{n} \frac{\sin(\omega_j \sigma/2)}{\sin(\omega_j T/2)} Q_{p,j} Q_{n,j} = 2 \frac{d\phi_{p,T}(\sigma)}{dt} \quad (7.9)
$$

where $\phi_{p,T}(\sigma)$ and $\psi_{p,T}(\sigma)$ denote the $p$ element of $\phi_T(\sigma)$ and $\psi_T(\sigma)$, respectively.

**7.1.4. Components of the governing matrices** Using the definition of **Σ** with the values of **r** and **w** found in Eq. 7.4, it comes that **Σ** is populated by the following elements with $i \le j \le k$:

$$
\Sigma_{i,j} = \frac{1}{\sqrt{m_n}} \sum_{i=1}^{n} \frac{\cos(\omega_i (t_i - t_j + T/2))}{\omega_i \sin(\omega_i T/2)} Q_{n,i}^2 = -\phi_{n,T}(2(t_i - t_j) + T). \quad (7.10a)
$$

The entries of matrix **Π** reduce for *i* < *j* ≤ *k* to

$$
\Pi_{i,j} = \Upsilon_T(2(t_j - t_i) - T) = \sum_{i=1}^{n} \frac{\sin(\omega_i (t_i - t_j + T/2))}{\sin(\omega_i T/2)} Q_{n,i}^2 \quad (7.11)
$$

because $r_{i+n} = Q_{n,i}$ (see Eq. 7.4).

**7.1.5. Physical interpretation of the kernel isomorphic to the invariant subspace of u** The *n*-th velocity just before the *j*-th impact is given by

$$
v_n(t_j) = \mathbf{e}_{2n}^{\top} \mathbf{P}^{-1} \mathbf{z}_j = \frac{1}{\sqrt{m_n}} \mathbf{r}^{\top} \mathbf{z}_j \quad (7.12)
$$

where $\mathbf{z}_j$ denotes the state of the system just before the $j$-th impact written in the new basis: $\mathbf{z}_j = \mathbf{S}(\sigma_j)\mathbf{N}\dots\mathbf{N}\mathbf{S}(\sigma_1)\mathbf{z}_0$. However, $\boldsymbol{\lambda} = \mathbf{L}_k\mathbf{z}_0$ implies $\boldsymbol{\lambda}_j = \mathbf{r}^\top\mathbf{z}_j$ and then

$$
\lambda_j = \sqrt{m_n} \bar{v}_n(t_j). \tag{7.13}
$$

This provides a physical interpretation of $\boldsymbol{\lambda}$: its component $j$ is proportional to the velocity just before the
$j$-th impact. In particular, $\boldsymbol{\lambda}_j$ and $\bar{v}_n(t_j)$ have the same sign. This observation suggests an efficient way
of detecting potential solutions which are not admissible solutions, since admissible solutions necessarily
have non-negative pre-impact velocities.
## 7.2. One-impact-per-period solutions
In order to facilitate the coming discussions, the following definition of harmonic is used.

**Definition 7.1 [Harmonic]** A harmonic is a multiple of $T_i := 2\pi/\omega_i$ for some $i \in [1, n]$.

### 7.2.1. Periodicity: Theorem 3.1 (NC1)
If $k = 1$, $t_1 = T$ is the period. The skew-symmetric matrix $\Pi$ is of dimension $k = 1$, hence $\Pi = [0]$ and $j = [1]$, $\lambda = [\lambda]$. Any $\lambda \in \mathbb{R}$ satisfies $\Pi\lambda = 0$ and Theorem 3.1 stipulates that potential $x_0$ lie at most in a one-dimensional manifold.

### 7.2.2. Potential solutions: Theorem 3.2 (NC2)
According to Theorem 3.2, NC2 can be expressed as $\Sigma(T)\lambda = -g_0j = -dj$, which gives a condition on $\lambda = [\lambda]$. Here, $\Sigma = [-\phi_{n,T}(T)]$ and the condition becomes $\phi_{n,T}(T)\lambda = d$. The existence and uniqueness of potential solutions is determined by the values of $\phi_{n,T}(T)$ and the nullity of $d$.

* If $d = 0$, there is generically no non-trivial potential solution. If $T$ is solution of $\phi_{n,T}(T) = 0$, which occurs for a countable set of values, there is an infinity of potential solutions with the same period $T$ given by $\varphi(\lambda)$, $\lambda \in \mathbb{R}$.

* If $d \neq 0$, there is no potential solution if $\phi_{n,T}(T) = 0$ and a unique potential solution corresponding to $\lambda = d/\phi_{n,T}(T)$ whenever $\phi_{n,T}(T) \neq 0$. Generically, $\phi_{n,T}(T) \neq 0$ hence, as implied by Theorem 3.3, there is locally a one-parameter family of periodic trajectories, parametrized by $T$.

### 7.2.3. Admissible solutions for a non-vanishing unconstrained gap at rest
We focus on the case $d \neq 0$ in the neighbourhood of a $T$ such that $\phi_{n,T}(T) \neq 0$. The unique potential solution for $T$ is given by $\varphi(d/\phi_{n,T}(T))$. It is an admissible solution if the gap remains non-negative (NC3). Following the results from paragraph 7.1.5, a significant proportion of time intervals can be easily excluded using the necessary (but not sufficient) criterion $\lambda \ge 0$, or equivalently

$$ \operatorname{sign}(d) = \operatorname{sign}(\phi_{n,T}(T)). \qquad (7.14) $$

The sign of $\phi_{n,T}(T)$ is thus a meaningful quantity to track, in order to exclude straightforwardly non-admissible solutions. From Eq. (7.8), the following holds:

$$ \phi_{n,T}(T) = -\frac{1}{\sqrt{m_n}} \sum_{j=1}^{n} \cot(\omega_j T/2) \frac{Q_{n,j}^2}{\omega_j}. \qquad (7.15) $$

Since cot is a decreasing function on any of its intervals of definition, $T \mapsto \phi_{n,T}(T)$ increases on any open interval bounded by two consecutive harmonics. On such an interval, its left limit is $-\infty$ and its right limit is $+\infty$. It is continuous, hence necessarily has a unique root $T_0$ on every such interval. The condition $\lambda \ge 0$ comes down to seeking periods $T$ of admissible solutions in intervals delimited:

* on the left, by a harmonic and on the right, by the next root $T_0$ of $\phi_{n,T}(T)$, if $d > 0$;

* on the left, by a root $T_0$ of $\phi_{n,T}(T)$ and the right, by the next harmonic, if $d < 0$.

The sufficient condition NC3 ($x_n \le d$ on $(0, T)$ with $x_n$ given by Eq. (7.5)) has to be tested numerically on every such intervals. If it is satisfied for one $T$, then $\varphi(d/\phi_{n,T}(T))$ is the unique admissible solution and there is a one-parameter family of admissible solutions in the neighbourhood of $T$ through Corollary 3.1.

It is worthy to note that there is a complementarity between the open gap ($d > 0$) and pre-stressed system ($d < 0$) configurations, in the sense that for any given period, the impact velocity is always compatible with the unilateral constraints for either $d > 0$ or $d < 0$, but not simultaneously. This complementarity is illustrated in Fig. 6.

### 7.2.4. Closed-form expression of the nonsmooth modes and illustrations
From Eq. (7.5) consequence of Theorem 3.4, the positions and the velocities are known explicitly: $\forall \tau \in [0, T]$,

$$ x_p(\tau) = -\frac{\lambda}{\sqrt{m_p}} \sum_{j=1}^{n} \frac{\cos(\omega_j(\tau - T/2))}{\omega_j \sin(\omega_j T/2)} Q_{p,j} Q_{n,j} \qquad (7.16a) $$

$$ \dot{x}_p(\tau) = -\frac{\lambda}{\sqrt{m_p}} \sum_{j=1}^{n} \frac{\sin(\omega_j(\tau - T/2))}{\sin(\omega_j T/2)} Q_{p,j} Q_{n,j}. \qquad (7.16b) $$

On an interval of periods corresponding to admissible solutions, Eq. (7.16) gives a closed-form expression of 1-ipp nonsmooth modes parametrized by the period $T$ and time $\tau$. The motion of the $n$-th mass is
<figure><img src="image_12.png" /><figcaption>**Figure 6:** Graph of $T \mapsto \phi_{n,T}(T)$ [---]. Harmonics [----]. Grey area [ ] corresponds to $\lambda < 0$, i.e. zones with no admissible solutions if $d > 0$. Orange area [ ] corresponds to $\lambda > 0$, i.e. zones with no admissible solutions if $d < 0$. ① and ⑥ delimit an interval of periods corresponding to admissible solutions, see Fig. 7.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 7:** Continuum of solutions with 1 impact per period. All solutions correspond to positive pre-impact velocities ($\phi_{n,T}(T) > 0$), located in the hatched orange zone of Fig. 6.</figcaption></figure>

plotted in Fig. 7 for several periods in the interval delimited by ① and ⑥ in Fig. 6. The nonsmooth mode corresponding to this continuum of admissible solutions is projected in the $(x_{n-1}, x_n, \dot{x}_n)$ space in Fig. 8.

**Remark 7.1.**

* The sequence of free flights is 1-periodic, hence Corollary 6.1 with $a = 0$ implies that $t = t_0/2 = 0$ and $t = t_1/2 = T/2$ are axes of symmetry for the positions, hence the velocities are symmetric with respect to the points $(0, 0)$ and $(T/2, 0)$. This is clear in Fig. 7 and corresponds to the symmetry with respect to the plane $\dot{x}_n = 0$ in Fig. 8.

* All the velocities $\dot{x}_p$ vanish at $\tau = T/2$, which is a direct consequence of the symmetry.

* For $p \neq n$, $\dot{x}_p^+(0) = \dot{x}_p^-(T)$ because the lines of $Q$ are orthogonal: $\sum_j Q_{p,j} Q_{n,j} = 0$ which implies from Eq. (7.16) the continuity of positions of the non-constrained masses.

* For $p = n$, $\dot{x}_n(0^+) = \lambda_1/m_n = -\dot{x}_n(T^-)$ which corresponds to the impact law.

* For larger $T$, the probability of meeting the obstacle between impacts increases together with the number of oscillations during one period. This indicates that the density of nonsmooth modes decreases with the period.

## 7.3. Two-impact-per-period solutions

**7.3.1. Invariance by $u$: Theorem 3.1 (NC1)** The impact times are $t_1, t_2$, $T = t_2$ being the period. The free flight durations are $\sigma_1 = t_1$ and $\sigma_2 = T - t_1$. The skew-symmetric matrix $2 \times 2$-Π reads:

$$ \Pi = \begin{bmatrix} 0 & a \\ -a & 0 \end{bmatrix} \qquad (7.17) $$
<figure><img src="image_14.png" /><figcaption>**Figure 8:** Projection in $(x_{n-1}, x_n, \dot{x}_n)$ of the 1-ipp nonsmooth mode in the neighbourhood of $4T_3$. The plots correspond to various intervals of $T$. The green ellipse shows the motion ⑥ in Fig. 7, i.e. the grazing orbit of the third linear mode. The red curves correspond to motions ⑤, ③, ②, and ①. The yellow curve highlights the intersection between the nonsmooth mode and the hyperplane $x_n = d$.</figcaption></figure>

with

$$a = \Upsilon_{\sigma_1+\sigma_2}(\sigma_2 - \sigma_1) = \sum_{j=1}^{n} \frac{\sin(\frac{1}{2}\omega_j(\sigma_2 - \sigma_1))}{\sin(\frac{1}{2}\omega_j(\sigma_1 + \sigma_2))} Q_{n,j}^2 = \sqrt{m_n}\psi_{n,\sigma_1+\sigma_2}(\sigma_2 - \sigma_1). \quad (7.18)$$

Using Theorem 3.1, NC1 comes down to $\Pi\lambda = 0$ hence $\lambda = 0$ or $a = 0$. Non-trivial solutions can be found only if $a = 0$, which yields the following nonlinear equation in $\sigma_1, \sigma_2$:

$$\psi_{n,\sigma_1+\sigma_2}(\sigma_2 - \sigma_1) = 0. \quad (7.19)$$

The set of roots of $a$ is denoted by $\mathcal{S}$.

**Solution curves of Eq. (7.19)** There is generally no closed-form expressions for the solutions to Eq. (7.19). It is yet possible to describe the behaviour of the roots of $\psi_{n,\sigma_1+\sigma_2}(\sigma_2 - \sigma_1)$ when $\sigma_1 + \sigma_2$ is in the neighbourhood of a harmonic. Indeed, the function $(\sigma_1, \sigma_2) \mapsto \psi_{n,\sigma_1+\sigma_2}(\sigma_2 - \sigma_1)$ is defined when $\sigma_1 + \sigma_2$ is not a harmonic. It shares the roots of

$$\tilde{a}(\sigma_1, \sigma_2) := a(\sigma_1, \sigma_2) \prod_{j=1}^{n} \sin(\frac{1}{2}\omega_j(\sigma_1 + \sigma_2)) \quad (7.20)$$

which is a linear combination of $\sin(\omega_j(\sigma_2 - \sigma_1)/2)$ for fixed $T$, hence the equation $\tilde{a}(\sigma_1, \sigma_2) = 0$ is easier to solve numerically. Roots of $\tilde{a}$ are depicted in Fig. 9.

Obvious roots are given by $\sigma_1 = \sigma_2$, which corresponds to NC1 for 1-ipp motions as 2-ipp motions. The set $\sigma_2 = \sigma_1$ emerges clearly in Fig. 9. Some other particular roots can be very easily computed. When $\sigma_1 + \sigma_2$ is in the vicinity of a harmonic $pT_i$, $a(\sigma_1, \sigma_2)$ is dominated by the $i$-th element of the sum and its roots are in the neighbourhood of the roots of

$$\sin(\frac{1}{2}\omega_i(\sigma_2 - \sigma_1)) = \sin(\frac{1}{2}\omega_i(T - 2\sigma_1)) = \sin(\frac{1}{2}\omega_i T)\cos(\omega_i\sigma_1) - \cos(\frac{1}{2}\omega_i T)\sin(\omega_i\sigma_1) \quad (7.21)$$
<figure><img src="image_15.png" /><figcaption>**Figure 9:** Level set 0 of $\Upsilon_{\sigma_1+\sigma_2}(\sigma_2 - \sigma_1)$, which corresponds to potential solutions [—]. Line $\sigma_1 + \sigma_2 = 3T_3$ [—]. *Lower-right triangle:* Curve corresponding to an continuum of admissible solutions [—]. ① and ⑥ correspond to the curves reported in Fig. 10. *Upper-left triangle:* Dashed lines are iso-periods when $\sigma_1 + \sigma_2$ is harmonic. The dots correspond to the intersection between the roots and iso-period lines and they are regularly spaced on each line.</figcaption></figure>

which converges to $-\sin(\omega_i \sigma_1)$ when $T$ approaches $pT_i$. After all, the solutions of Eq. (7.19) are such that $\sigma_1 = pT_i/2$ with $p \in \mathbb{N}^*$ when $\sigma_1 + \sigma_2 \to pT_i$. The intersections between the harmonics lines $\sigma_2 = qT_i - \sigma_1$ and the roots of $a$ is the set

$$ \left\{ \left( \frac{p}{2} T_i, \left(q - \frac{p}{2}\right) T_i \right), i, p, q \in \mathbb{N}^* \right\}. \qquad (7.22) $$

This is illustrated by coloured dots in Fig. 9 and produces obvious roots of $\tilde{a}$. This result will become very handy in view of subsection 7.3.3.

**Remark 7.2.**

* Because $a(\sigma_1, \sigma_2) = -a(\sigma_2, \sigma_1)$, $\sigma_1 = \sigma_2$ is an axis of symmetry for the roots of $a$ (see Fig. 9).

* Since $\tilde{a}$ is analytic, the implicit function theorem applies almost everywhere and it is therefore possible to calculate locally a parametrization $T(\sigma_1)$.

* Generically, the roots of $a$ are arranged along curves, in agreement with Theorem 3.3.

### 7.3.2. Potential solutions: Theorem 3.2 (NC2)
The set of times of impact $(t_1, t_2)$ for which there exists a potential solution is a subset of the set of the roots of $a$. Using Theorem 3.2, NC2 writes

$$ \Sigma \lambda = -d\mathbf{j} \quad \text{where} \quad \Sigma = \begin{bmatrix} \alpha & \beta \\ \beta & \alpha \end{bmatrix}, \quad \mathbf{j} = \begin{bmatrix} 1 \\ 1 \end{bmatrix} \qquad (7.23) $$

with $\alpha = -\phi_{n,\sigma_1+\sigma_2}(\sigma_1 + \sigma_2)$ and $\beta = -\phi_{n,\sigma_1+\sigma_2}(\sigma_2 - \sigma_1)$. For $\sigma_1, \sigma_2$ solutions of Eq. (7.19), the existence and uniqueness of potential solutions is determined by the value of $\alpha, \beta$ and the nullity of $d$. The following illustrates Theorem 3.3 for $k = 2$.

**Non-vanishing algebraic gap $d \neq 0$**

* If $\alpha = \beta \neq 0$, any $\lambda$ such that $\lambda_2 = -\lambda_1 - d/\alpha$ satisfies (7.23), hence a solution space of dimension 1.

* If $\alpha = -\beta$, $\Sigma \lambda = -d\mathbf{j}$ has no solution.

* If $\alpha \neq \beta$, $\Sigma$ is invertible and there is a unique solution $\lambda = -d\Sigma^{-1}\mathbf{j}$. This is the generic case.

**Case $d = 0$**

* If $\alpha = \beta = 0$, then $\lambda$ satisfies $\Sigma \lambda = -d\mathbf{j}$ and the solution space is two-dimensional.

* If $\alpha = \beta \neq 0$, any $\lambda$ such that $\lambda_1 = -\lambda_2$ satisfies NC2 so the set of solutions is of dimension 1.

* If $\alpha = -\beta \neq 0$, $\lambda_1 = \lambda_2$ is sufficient and the solution space if of dimension 1.

* If $\alpha \neq \beta$, $\lambda = -d\Sigma^{-1}\mathbf{j} = 0$ and there is no non-trivial solution (generic case).

For $d \neq 0$, admissible solutions of the generic case $\alpha \neq \beta$ and the non-generic case $\alpha = \beta \neq 0$ are investigated in depth in subsection 7.3.3.

### 7.3.3. Admissible solution when $d \neq 0$
**Generic case** $\alpha \neq \beta$ The general setting $\alpha \neq \beta$ is considered and $\lambda$ satisfies both $\Pi\lambda = 0$ and $\Sigma\lambda = -d\mathbf{j}$. Hence, $\lambda = -d\Sigma^{-1}\mathbf{j}$ but $\mathbf{j}$ is obviously an eigenvector of $\Sigma$ corresponding to the eigenvalue $\alpha + \beta$ so $\lambda_1 = \lambda_2 := \lambda$ and $\lambda = 1/\alpha + \beta$. The physical interpretation is that the mass $n$ has the same impact velocity at $t_1$ and $T$. The expressions of the positions and velocities result from Eq. (7.5): $\forall \tau \in [0, \sigma_1]$,

$$ \begin{cases} x_p(\tau) = \lambda(\phi_{p,\sigma_1+\sigma_2}(2\tau + \sigma_2 - \sigma_1) + \phi_{p,\sigma_1+\sigma_2}(2\tau - \sigma_1 - \sigma_2)) \\ \dot{x}_p(\tau) = \lambda(\psi_{p,\sigma_1+\sigma_2}(2\tau + \sigma_2 - \sigma_1) + \psi_{p,\sigma_1+\sigma_2}(2\tau - \sigma_1 - \sigma_2)) \end{cases} \quad (7.24a) $$

$$ \dot{x}_p(\tau) = \lambda(\psi_{p,\sigma_1+\sigma_2}(2\tau + \sigma_2 - \sigma_1) + \psi_{p,\sigma_1+\sigma_2}(2\tau - \sigma_1 - \sigma_2)) \quad (7.24b) $$

and $\forall \tau \in [0, \sigma_2]$,

$$ \begin{cases} x_p(\tau + \sigma_1) = \lambda(\phi_{p,\sigma_1+\sigma_2}(2\tau - \sigma_1 - \sigma_2) + \phi_{p,\sigma_1+\sigma_2}(2\tau + \sigma_1 - \sigma_2)) \\ \dot{x}_p(\tau + \sigma_1) = \lambda(\psi_{p,\sigma_1+\sigma_2}(2\tau - \sigma_1 - \sigma_2) + \psi_{p,\sigma_1+\sigma_2}(2\tau + \sigma_1 - \sigma_2)). \end{cases} \quad (7.25a) \quad (7.25b) $$

$$ \dot{x}_p(\tau + \sigma_1) = \lambda(\psi_{p,\sigma_1+\sigma_2}(2\tau - \sigma_1 - \sigma_2) + \psi_{p,\sigma_1+\sigma_2}(2\tau + \sigma_1 - \sigma_2)). \quad (7.25b) $$

As previously, the admissible solutions necessarily correspond to non-negative pre-impact velocities: $\lambda_1 \ge 0$, $\lambda_2 \ge 0$ which reduces to $\lambda \ge 0$. Besides, the sign of $\lambda$ switches when $T$ crosses a harmonic. This sign is the product between the sign of $d$ and the sign of

$$ \alpha + \beta = -\phi_{n,T}(T) - \phi_{n,T}(T - 2t_1) = \frac{1}{\sqrt{m_n}} \sum_{j=1}^{n} \frac{Q_{n,j}^2 \cos\left(\frac{\omega_j}{2}T\right) + \cos\left(\frac{\omega_j}{2}(T - 2t_1)\right)}{\omega_j \sin\left(\frac{\omega_j}{2}T\right)}. \quad (7.26) $$

When $T$ is in a punctured neighbourhood of a harmonic $pT_i$ and $t_1$ is not a multiple of $T_i$, the term $i$ of the above sum dominates since $\sin(\omega_i T/2) \to 0$ when $T \to T_i$. It tends to $+\infty$ when $T \to T_i^+$ and $-\infty$ when $T \to T_i^-$. Therefore, when $d > 0$, $\lambda \ge 0$ for $T$ in the left-neighbourhood of a harmonic and $\lambda \le 0$ in the right-neighbourhood. When $d < 0$, $\lambda \ge 0$ in the right-neighbourhood of a harmonic and $\lambda \le 0$ in the left-neighbourhood.

**Remark 7.3.** Recall that the intersection points between the set $S$ of impact times leading to potential solutions and harmonic lines ($\sigma_1 + \sigma_2 = pT_i$) are given by Eq. (7.22). One strategy to search for admissible solutions is hence to check the admissibility of potential solutions in the left-neighbourhood (resp. right-neighbourhood) of these points if $d > 0$ (resp. $d < 0$).

<figure><img src="image_16.png" /><figcaption>**Figure 10:** Positions of the $n$-th mass of a continuum of admissible solutions with 2 impacts per period. Circled numbers refer to Figs. 9 and 11. Impact times [---]. Axes of symmetry [---].</figcaption></figure>

The position $x_n$ of a continuum of admissible solutions obtained for $T$ in the left-neighbourhood of $3T_3$ is pictured in Fig. 10. In Fig. 9, the admissible solution ⑥ is indeed at the intersection of the line $\sigma_2 = 3T_3 - \sigma_1$ and $S$. This continuum constitutes a nonsmooth mode illustrated in Fig. 11.

**Remark 7.4.** Theorem 6.4 holds here because the sequence of free-flight durations $\dots, \sigma_1, \sigma_2, \dots$ is invariant by the transformation $n \mapsto 2 - n$. The sets $t = t_1/2$ and $t = (t_1 + T)/2$ are indeed axes of symmetry for every $x_p$, in particular for $x_n$ as shown in Fig. 10.
<figure><img src="image_17.png" /><figcaption>Figure 11: 2-ipp nonsmooth mode. Circled numbers refer to Figs. 9 and 10.</figcaption></figure>

**Non-generic case $\alpha = \beta \neq 0$** As opposed to the generic case for which continua of solutions emanate from a curve in the space of impact times $(t_1, T)$ (or $(\sigma_1, \sigma_1 + \sigma_2)$), when $\alpha = \beta \neq 0$ continua emanate from a single point in this space. The condition $\Sigma\lambda = -d\mathbf{j}$ reads $\lambda_1 + \lambda_2 = -d/\alpha$ and the expression of the positions stems from Eq. (7.5):

$$x_p(\tau) = \begin{cases} \lambda_1\phi_{p,T}(2\tau + \sigma_2 - \sigma_1) - \left(\frac{d}{\alpha} + \lambda_1\right)\phi_{p,T}(2\tau - T) & \text{if } \tau \in [0, \sigma_1] \\ \lambda_1\phi_{p,T}(2\tau - T) - \left(\frac{d}{\alpha} + \lambda_1\right)\phi_{p,T}(2\tau + \sigma_1 - \sigma_2) & \text{if } \tau \in (\sigma_1, T] \end{cases} \quad (7.27)$$

Examples of such solutions are proposed in Fig. 12, which illustrates a continuum of admissible solutions sharing the same impact times $t_1$ and $T$. In contrast to the generic case where $\Sigma$ is invertible, it shows that positions do not display any axis of symmetry, except the blue one ($\lambda_1 = 0$, axes at 0 and $T/2$), red one ($\lambda_2 = 0$, axes at $t_1$ and $t_1 + T/2$) and the purple one ($\lambda_1 = \lambda_2$, axes at $t_1/2$ and $(t_1 + T)/2$). This highlights the necessity of an invertible $\Sigma$ in Theorem 6.4.

The geometric properties of the position curves in Fig. 12 are now interpreted.

* The pre-impact velocities $\lambda_1$ and $-\lambda_1 - d/\alpha$ are bounded by 0 and $-d/\alpha$. Since $\dot{x}_n$ is a continuous function of $\lambda_1$, the pre-impact velocities continuously vary from $(0, -d/\alpha)$ (blue curve) to $(-d/\alpha, 0)$ (red curve), which are trajectories with two grazing contacts and one (non-grazing) impact. The sum of the two pre-impact velocities is constant.

* Even though every position $x_n$ has generically no symmetry, the set of all admissible positions for $(t_1, T)$ has two axes of symmetry: $t = t_1/2$ and $t = (t_1 + T)/2$. Indeed, if $\mathbf{x}$ is an admissible solution for $\lambda_1$, then, for $\tau \in [0, t_1]$,

$$x_p(\sigma_1 - \tau) = -\left(\frac{d}{\alpha} + \lambda_1\right)\phi_{n,T}(2\tau + \sigma_2 - \sigma_1) + \lambda_1\phi_{n,T}(2\tau - T), \quad (7.28)$$

hence $\tau \mapsto \mathbf{x}(\sigma - \tau)$ is an admissible solution for $-\lambda_1 - d/\alpha \ge 0$. This proves that $t = \sigma_1/2 = t_1/2$ is an axis of symmetry for all positions; similarly, $t = (t_1 + T)/2$ is another axis of symmetry.
<figure><img src="image_18.png" /><figcaption>**Figure 12:** Continuum of 2-ipp admissible solutions with constant impact times. The continuum is delimited by the blue and the red curves, which both correspond to one impact and two grazing contacts. The thick purple curve has two axes of symmetry. All the other curves have no axis of symmetry, because Σ is not invertible (see Theorem 6.4).</figcaption></figure>

If $λ_1 = -d/(2α)$, then $λ_2 = λ_1$ and these two axes are axes of symmetry for the corresponding position.

* It can also be shown that the two limiting configurations $λ_1 = 0$ and $λ_2 = 0$ are associated to the same solution, but with a time-shift of $σ_1$ (and by interchangeability of $σ_1$ and $σ_2$, of $σ_2$ too).

* The continuum of positions has fixed points, i.e. points common to all trajectories. For example,

$$x_p(t_1/2) = λ_1φ_{n,T}(σ_2) - (\frac{d}{α} + λ_1)φ_{p,T}(σ_1 - T) = -\frac{d}{α}φ_{n,T}(σ_2) \quad (7.29)$$

which does not depend on $λ_1$.

* Lastly, positions $x_p$ and velocities $y_p$ are linear functions of $λ_1$, for all $p \in [1, n]$. This means that the manifold of solutions is a ruled surface of the phase space.

<figure><img src="image_19.png" /><figcaption>**Figure 13:** Continuum of 2-ipp admissible solutions with constant impact times and $t_1 = T_0/2$. The continuum includes a 1-ipp admissible solution [—]: it is a *bridge* between 1-ipp and 2-ipp solutions.</figcaption></figure>

Interestingly, some manifolds with constant impact times exhibit the additional property of *joining* a 1-ipp manifold to a 2-ipp manifold. Such a situation is referred to as *bridge* between 1-ipp and 2-ipp and arises when $α = β$. The positions then reduce to

$$x_p(τ) = \begin{cases} λ_1φ_{p,T}(2τ) - (\frac{d}{α} + λ_1)φ_{p,T}(2τ - T) & \text{if } τ \in [0, T/2], \\ λ_1φ_{p,T}(2τ - T) - (\frac{d}{α} + λ_1)φ_{p,T}(2τ) & \text{if } τ \in (T/2, T]. \end{cases} \quad (7.30)$$

An illustration is provided in Fig. 13. The corresponding manifold is shown in Fig. 14. Such solutions have the following attributes:

* when $λ_1 = -d/(2α)$, the solution is certainly a 1-ipp of period $T/2$,

* if this 1-ipp solution is admissible, it is generically part of a 1-ipp nonsmooth mode (Theorem 3.3), hence its trajectory both belongs to a 2-ipp NSM and a 1-ipp NSM,

* all positions have an axis of symmetry at $t = T_0/2$.
<figure><img src="image_20.png" /><figcaption>**Figure 14:** NSM bridge connecting 1-ipp and 2-ipp manifolds. All trajectories of this manifold share the same period. Opaque: part corresponding to $[0, T_0/2]$. Transparent: part corresponding to $[T_0/2, T_0]$. 1-ipp solution [—]. The purple curve surface intersects 1-ipp nonsmooth mode which is not represented for clarity purposes.</figcaption></figure>

**7.4. Multiple-impact-per-period solutions** In the more general case $k > 2$, the principle is identical to $k = 1$ or $k = 2$. Numerical difficulties emerge from the higher number of equations to be solved. From Theorems 3.1 and 3.2, it is known that potential solutions can be sought as impact times $s = (t_1, \dots, t_k)$ such that there exists a $\lambda \in \mathbb{R}^k$ solution of $\Pi(s)\lambda = 0$ and $\Sigma(s)\lambda = -d\mathbf{j}$. Following up on Remark 3.1 when $k$ is odd, then generically rank($\Pi(s)$) = $k-1$ and the kernel of $\Pi$ is spanned by the vector $\mathbf{k}(s)$ of the Pfaffians of the diagonal $k-1 \times k-1$ submatrices of $\Pi$ with alternating signs. It therefore suffices to find $s$ such that the $k$ lines of $\Sigma(s)\mathbf{k}(s)$ are all equal to the same arbitrary value $\mu \neq 0$ ($k-1$ equations) and for such an $s$, $\lambda$ is given by $\lambda = -d\mathbf{k}(s)/\mu$. In practice, one possibility to find admissible solutions is to seek $s$ with a root-finding numerical solver starting from a random initial guess. The vector $\lambda = -d\Sigma^{-1}(s)\mathbf{j}$ corresponding to the numerical solution $s$ must have only non-negative components. If satisfied, the only remaining condition is that the unilateral constraint should not be violated during the free flights, which has to be checked numerically as well. When such a point is found, Theorem 3.3 ensures that $s$ is a point of a curve $S$ of impact times corresponding to admissible solutions. This curve can be computed using continuation techniques [3].

<figure><img src="image_21.png" /><figcaption>**Figure 15:** Continuum of 7-ipp admissible solutions. Bottom-left: symmetric motion of free-flight durations $c, b, a, b, c, d, d$. The bottom horizontal arrows in ① and ④ illustrate the relationship between the two trajectories: they correspond to the same orbit but with different time-orientation.</figcaption></figure>

Fig. 15 depicts four positions of a 7-ipp nonsmooth mode, identified through the procedure described above. Diagrams ① and ④ represent limiting cases of the admissible curve: indeed, they both have a
<figure><img src="image_22.png" /><figcaption>**Figure 16:** Examples of neutrally stable solutions after 1000 periods for 1 ipp (top), 2 ipp (middle), 3 ipp (bottom). [—] Modal trajectory. [......] Perturbed solutions.</figcaption></figure>

zero pre-impact velocity, which would become negative if $s$ followed the curve of potential solutions. In other words, ① and ④ each have a grazing contact (the fifth and the seventh, respectively); on the outer-neighbourhood of the admissible domain, the motion violates the contact condition.

Interestingly, ③ exhibits an axis of symmetry. Indeed, the periodic sequence of free-flight durations $\sigma_1, \sigma_2, \dots, \sigma_7$ is of the form $c, b, a, b, c, d, d$, which is invariant by the transformation $n \mapsto 6-n$. Invoking Theorem 6.4 (or its Corollary 6.2), it follows that the positions have two axes of symmetry per period: $t = (t_2 + t_3)/6$ and $t = t_6$, as observed in diagram ③.

Lastly, ① and ④ correspond to the same curve in the phase space, one parametrized in increasing time and the other one, in decreasing time (see arrows in Fig. 15, ① and ④).

**8. Insight on stability of nonsmooth modes**
Stability of the periodic trajectories is not thoroughly investigated in this work, however we give a few indications on how to perform the stability analysis and illustrate a few neutrally stable orbits. Given an initial condition $\mathbf{x}_0$ leading to a $k$-ipp periodic solution $\mathbf{x}$, we choose the Poincaré section $g(\mathbf{x}) = 0$. The question is to quantify the effect of a perturbation $\delta\mathbf{x}_0$ of $\mathbf{x}_0$ on the state when $\mathbf{x}$ returns to the Poincaré map after $k$ impacts—assuming such a state exists in the neighbourhood of $\mathbf{x}(T)$. This corresponds to the mapping

$$ \mathbf{x}_0 + \delta\mathbf{x}_0 \mapsto \tilde{\mathbf{N}}\tilde{\mathbf{S}}(\sigma_k + \delta\sigma_k)\tilde{\mathbf{N}}\dots\tilde{\mathbf{N}}\tilde{\mathbf{S}}(\sigma_1 + \delta\sigma_1)(\mathbf{x}_0 + \delta\mathbf{x}_0) \quad (8.1) $$

where $\delta\sigma_l$ is an unknown change of duration of the $l$-th free flight. The first-order Taylor expansion of this assumed smooth mapping yields an equation of the form

$$ \delta \mathbf{x} = \tilde{\mathbf{N}} \tilde{\mathbf{S}}(\sigma_k) \tilde{\mathbf{N}} \dots \tilde{\mathbf{N}} \tilde{\mathbf{S}}(\sigma_1) \delta \mathbf{x}_0 + \left( \sum_{l=1}^{k} \tilde{\mathbf{N}} \tilde{\mathbf{S}}(\sigma_k) \dots \tilde{\mathbf{N}} \mathbf{S}'(\sigma_l) \tilde{\mathbf{N}} \dots \tilde{\mathbf{N}} \tilde{\mathbf{S}}(\sigma_1) \delta \sigma_l \right) \mathbf{x}_0. \quad (8.2) $$

The unknowns $\delta\sigma_1, \dots, \delta\sigma_k$ can be found by solving the linear system provided by the $k$ equations $\tilde{g}(\mathbf{x}((\sigma_1 + \delta\sigma_1) + \dots + (\sigma_l + \delta\sigma_l))) = 0$ so that in the end, the perturbation $\delta\mathbf{x}$ is related to $\delta\mathbf{x}_0$ by a matrix $\mathbf{A}(\mathbf{x}_0)$ such that

$$ \delta \mathbf{x} = \mathbf{A}(\mathbf{x}_0) \delta \mathbf{x}_0. \quad (8.3) $$

The eigenvalues of $\mathbf{A}(\mathbf{x}_0)$ determine the spectral stability of the periodic solutions emanating from $\mathbf{x}_0$. As in [20], the spectral radius $\rho$ of $\mathbf{A}(\mathbf{x}_0)$ is observed to be at least 1. If $\rho > 1$ for some $\mathbf{x}_0$, then the
## Nomenclature

Quantities related to the original basis

<table><tr><td>N<sub>2</sub></td><td>mapping of the impact law v<sup>+</sup> = N<sub>2</sub>v<sup>-</sup></td><td>n × n</td></tr><tr><td>r<sub>2</sub></td><td>normal vector to the hyperplane of the impact law (Prop. 2.1)</td><td>n × 1</td></tr><tr><td>u, v</td><td>generalized positions, generalized velocities</td><td>n × 1</td></tr><tr><td>Ñ</td><td>mapping of the impact law x<sup>+</sup> = Ñx<sup>-</sup></td><td>2n × 2n</td></tr><tr><td>S̃</td><td>mapping of a free flight of duration τ: x(τ + t) = S̃(τ)x(t)</td><td>2n × 2n</td></tr><tr><td>w<sub>1</sub></td><td>gradient of the gap g(u) = w<sub>1</sub><sup>T</sup>u + g<sub>0</sub></td><td>n × 1</td></tr><tr><td>x = [u v]<sup>T</sup></td><td>state vector</td><td>2n × 1</td></tr><tr><td>g(u), g̃(x)</td><td>gap in terms of positions u, of state x (g(u) = g̃(x))</td><td>1 × 1</td></tr></table>

### Quantities involved in the change of basis

<table><tr><td>P</td><td>change of basis matrix defined in (4.3) such that P<sub>x</sub> = z</td><td>2n × 2n</td></tr><tr><td>Q, Ω</td><td>orthogonal and diagonal matrices from V<sup>-1</sup>KV<sup>-T</sup> = QΩ<sup>2</sup>Q<sup>T</sup></td><td>n × n</td></tr><tr><td>V</td><td>triangular matrix from the Cholesky decomposition M = VV<sup>T</sup></td><td>n × n</td></tr></table>

### Quantities related to the new basis

<table><tr><td>N</td><td>mapping of the impact law z<sup>+</sup> = Nz<sup>-</sup> (N = I - 2rr<sup>T</sup>)</td><td>2n × 2n</td></tr><tr><td>r</td><td>normal vector to the hyperplane r = [0 Q<sup>T</sup>r<sub>2</sub>]<sup>T</sup></td><td>2n × 1</td></tr><tr><td>S</td><td>mapping of a free flight of duration τ: z(τ + t) = S(τ)z(t)</td><td>2n × 2n</td></tr><tr><td>w</td><td>gradient of the gap w = P<sup>-1</sup>[w<sub>1</sub> 0]<sup>T</sup></td><td>2n × 1</td></tr><tr><td>z</td><td>state vector</td><td>2n × 2n</td></tr></table>

periodic solutions induced by $x_0$ is unstable. If $\rho = 1$, it is neutrally stable. First numerical investigations show that both cases can be observed. Neutrally stable solutions are depicted in Fig. 16 for the system of section 7 for different numbers of impact per periods ($k \in \{1, 2, 3\}$). The motion stemming from a perturbed initial condition is shown to consist of the modal trajectory plus a small variation and a time shift.

## 9. Concluding remarks
Periodic solutions of n-degree-of-freedom nonsmooth systems i) involving piecewise-linear dynamics and ii) subjected to a single affine unilateral constraint have been investigated. Apart from the two above assumptions, the formulation is general and embraces non-diagonal mass matrices as in Finite Element Models, three dimensional systems, two impacting elastic bodies or self-contact, provided condition ii) holds, etc. Such periodic trajectories have been shown to be governed by $k$ equations, where $k$ is the considered number of impacts per period. They are organized as one-parameter continuous families, also called manifolds which commonly define nonsmooth modes of vibration. Symmetries in the sequence of free-flight durations have been shown to induce symmetries of the corresponding nonsmooth modes, seen as manifolds of the phase space.

Overall, the present work contributes to characterizing the vibratory properties of nonsmooth systems and directly makes possible the systematic computation of the response spectrum of a piecewise-linear vibro-impact oscillator, including branches of multiple-impact solutions. Results have been illustrated on a 5-dof spring-mass system. It was shown that some periodic solutions are neutrally stable, meaning that they can be potentially encountered and observed in a wide variety of real-world applications.

Future work will include a more detailed stability analysis and, ideally, an experimental comparison between the measured response spectrum of a vibro-impact oscillator and the computed one.

## Acknowledgements
The first author would like to thank Remco LEINE for a helpful discussion on impact laws, Stéphane JUNCA for his valuable suggestions and David URMAN for pointing out some relevant literature. Indications provided by reviewer #3 were greatly appreciated. The third author would like to thank Steven SHAW for insightful discussions on stability analysis.

## 10. References

[1] ACARY, Vincent, BROGLIATO, Bernard. *Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics*. Volume 35. Springer, 2008. [ISBN: 9783540753926]. [10.1007/978-3-540-75392-6].

[2] ACARY, Vincent, JAMES, Guillaume, PÉRIGNON, Franck. "Periodic motions of coupled impact oscillators." *ENOC 2014 - 8th European Nonlinear Dynamics Conference*. Vienna, Austria, 2014. [hal-01059824].
[3] ALLGOWER, Eugene, GEORG, Kurt. *Numerical continuation methods: an introduction*. Volume 13. Springer Science & Business Media, 2012. [ISBN: 9783642647642]. [10.1007/978-3-642-61257-2].

[4] ATTAR, Mostafa, KARRECH, Ali, REGENAUER-LIEB, Klaus. *Non-linear modal analysis of structural components subjected to unilateral constraints*. Journal of Sound and Vibration 389:380-410, 2017. [hal-01513216].

[5] BALLARD, Patrick. *The dynamics of discrete mechanical systems with perfect unilateral constraints*. Archive for Rational Mechanics and Analysis 154(3):199-274, 2000. [hal-00111308].

[6] BERTI, Massimiliano. *Nonlinear Oscillations of Hamiltonian PDEs*. Progress in Nonlinear Differential Equations and Their Applications. Birkhäuser Boston, 2007. [ISBN: 9780817646806].

[7] CHATI, Mandar, RAND, Richard, MUKERHJEECHATI, Subrata. *Modal analysis of a cracked beam*. Journal of Sound and Vibration 207(2):249-270, 1997. [hal-01510818].

[8] DI BERNARDO, Mario, BUDD, Chris, CHAMPNEYS, Alan Richard, KOWALCZYK, Piotr. *Piecewise-smooth dynamical systems: theory and applications*. Volume 163. Springer Science & Business Media, 2008. [ISBN: 9781846287084]. [10.1007/978-1-84628-708-4].

[9] DOYEN, David, ERN, Alexandre, PIPERNO, Serge. *Time-integration schemes for the finite element dynamic signorini problem*. SIAM Journal on Scientific Computing 33(1):223-249, 2011. [hal-00440128].

[10] ECK, Christof, JARUSEK, Jiri, KRBEČ, Miroslav. *Unilateral contact problems: variational methods and existence theorems*. Volume 270. CRC Press, 2005. [ISBN: 9781574446296].

[11] FREDRIKSSON, Mats, NORDMARK, Arne. "On normal form calculations in impact oscillators." *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*. Volume 456. 1994. The Royal Society. 2000: 315-329. [hal-01510821].

[12] GENEDELMAN, Oleg. *Exact solutions for discrete breathers in a forced-damped chain*. Physical Review E 87(6):062911, 2013. [hal-01510820].

[13] GUCKENHEIMER, John, HOLMES, Philip. *Nonlinear oscillations, dynamical systems, and bifurcations of vector fields*. Volume 42. Springer Science & Business Media, 1983. [ISBN: 9781461211402]. [10.1007/978-1-4612-1140-2].

[14] HOUSEHOLDER, Alston. *Unitary triangularization of a nonsymmetric matrix*. Journal of the Association for Computing Machinery 5(4):339-342, 1958. [hal-01316095].

[15] JAYAPRAKASH, KR, STAROSVETSKY, Yuli, VAKAKIS, Alexander, PEETERS, Maxime, KERSCHEN, Gaetan. *Nonlinear normal modes and band zones in granular chains with no pre-compression*. Nonlinear Dynamics 63(3):359-385, 2011. [hal-01510822].

[16] KELLEY, Al. *The stable, center-stable, center, center-unstable, unstable manifolds*. Journal of Differential Equations 3(4):546-570, 1967. [hal-01306971].

[17] LAMARQUE, Claude-Henri, JANIN, Olivier. *Modal analysis of mechanical systems with impact non-linearities: limitations to a modal superposition*. Journal of Sound and Vibration 235(4):567-609, 2000. [hal-01510826].

[18] LEBEAU, Gilles, SCHATZMAN, Michelle. *A wave problem in a half-space with a unilateral constraint at the boundary*. Journal of differential equations 53(3):309-361, 1984. [hal-01294216].

[19] LEE, Young, NUCERA, Francesco, VAKAKIS, Alexander, MCFARLAND, Donald, BERGMAN, Lawrence. *Periodic orbits, damped transitions and targeted energy transfers in oscillators with vibro-impact attachments*. Physica D: Nonlinear Phenomena 238(18):1868-1896, 2009. [hal-01510829].

[20] LEGRAND, Mathias, JUNCA, Stéphane, HENG, Sokly. *Nonsmooth modal analysis of a N-degree-of-freedom system with a purely elastic impact law*. preprint, 2015. [hal-01185980].

[21] LI, Kuinian, DARBY, Antony. *Modelling a buffered impact damper system using a spring-damper model of impact*. Structural Control and Health Monitoring 16(3):287-302, 2009. [hal-01510823].

[22] LIU, Yang, PAVLOVSKAIA, Ekaterina, WIERCIGROCH, Marian, PENG, Zhike. *Forward and backward motion control of a vibro-impact capsule system*. International Journal of Non-Linear Mechanics 70:30-46, 2015. [hal-01510830].

[23] LIU, Yang, WIERCIGROCH, Marian, PAVLOVSKAIA, Ekaterina, YU, Hongnian. *Modelling of a vibro-impact capsule system*. International Journal of Mechanical Sciences 66:2-11, 2013. [hal-01510831].

[24] MIKHLIN, Yuri, VAKAKIS, Alexander, SALENGER, Gary. *Direct and inverse problems encountered in vibro-impact oscillations of a discrete system*. Journal of Sound and Vibration 216(2):227-250, 1998. [hal-01352131].

[25] MOORE, Douglas, SHAW, Steven. *The experimental response of an impacting pendulum system*. International journal of non-linear mechanics 25(1):1-16, 1990. [hal-01510832].

[26] MOUSSI, El Hadi, BELLIZZI, Sergio, COCHELIN, Bruno, NISTOR, Ionel. *Nonlinear normal modes of a two degree of freedom oscillator with a bilateral elastic stop*, 2012. [hal-00784528].

[27] MOUSSI, El Hadi, BELLIZZI, Sergio, COCHELIN, Bruno, NISTOR, Ionel. *Nonlinear normal modes of a two degrees-of-freedom piecewise linear system*. Mechanical Systems and Signal Processing 64:266-281, 2015. [hal-01160325].

[28] NATSIAVAS, Sotirios. *Dynamics of multiple-degree-of-freedom oscillators with colliding components*. Journal of Sound and Vibration 165(3):439-453, 1993. [hal-01510834].

[29] NORDMARK, Arne. *Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators*. Nonlinearity 14(6):1517, 2001. [hal-01297283].

[30] PASCAL, Madeleine. *Dynamics and stability of a two-degree-of-freedom oscillator with an elastic stop*. Journal of Computational and Nonlinear Dynamics 1(1):94-102, 2005. [hal-01294059].

[31] PAVLOVSKAIA, Ekaterina, WIERCIGROCH, Marian. *Periodic solution finder for an impact oscillator with a drift*. Journal of Sound and Vibration 267(4):893-911, 2003. [hal-01510819].

[32] PILIPCHUK, Valery. *Impact modes in discrete vibrating systems with rigid barriers*. International Journal of Non-Linear Mechanics 36(6):999-1012, 2001. [hal-01295245].

[33] POPPLEWELL, Neil, BAPAT, Charusheel, MCLACHLAN, Ken. *Stable periodic vibroimpacts of an oscillator*. Journal of Sound and Vibration 87(1):41-59, 1983. [hal-01510833].

[34] PUN, D, LAU, SL, LAW, SS, CAO, DQ. *Forced vibration analysis of a multidegree impact vibrator*. Journal of Sound and Vibration 213(3):447-466, 1998. [hal-01509022].
[35] SAMI, Masri. Steady-state response of a multidegree system with an impact damper. Journal of Applied Mechanics 40(1):127-132, 1973. [hal-01511958].

[36] SHAW, Steven Wayne, HOLMES, PJ. A periodically forced piecewise linear oscillator. Journal of Sound and Vibration 90(1):129-155, 1983. [hal-01509017].

[37] SHAW, Steven, PIERRE, Christophe. Non-linear normal modes and invariant manifolds. Journal of Sound and Vibration 150(1):170-173, 1991. [hal-01310674].

[38] STUDER, Christian, LEINE, Remco, GLOCKER, Christoph. Step size adjustment and extrapolation for time-stepping schemes in non-smooth dynamics. International journal for numerical methods in engineering 76(11):1747-1781, 2008. [hal-01316098].

[39] THORIN, Anders, DELEZOIDE, Pierre, LEGRAND, Mathias. Non-smooth modal analysis of piecewise-linear impact oscillators. [Hyperlink for the animations of the current paper], 2016. [hal-01298983].

[40] VAKAKIS, Alexander, MANEVITCH, Leonid, MIKHLIN, Yuri, PILIPCHUK, Valery, ZEVIN, Alexander. Normal modes and localization in nonlinear systems. Kluwer Academic Publishers, 1996. [ISBN: 9780792370109].

[41] VEDENOVA, EG, MANEVICH, Leonid, PILIPCHUK, Valery. Normal oscillations of a string with concentrated masses on non-linearly elastic supports. Journal of Applied Mathematics and Mechanics 49(2):153-159, 1985. [hal-01510835].

[42] VON GROLL, Götz, EWINS, David. The harmonic balance method with arc-length continuation in rotor/stator contact problems. Journal of sound and vibration 241(2):223-233, 2001. [hal-01333704].

[43] WAGG, David. Multiple Non-Smooth Events in Multi-Degree-of-Freedom Vibro-Impact Systems. Nonlinear Dynamics 43(1-2):137-148, 2006. [hal-01403572].

[44] WHISTON, GS. Global dynamics of a vibro-impacting linear oscillator. Journal of Sound and Vibration 118(3):395-424, 1987. [hal-01508924].

A. Dimensions of modes in the symmetric case The transformation $\tau$ is an involution since $\tau(\tau(i)) = a - (a - i - 1) - 1 = i$ and if $\mathbf{P}_{\tau}$ is the associated permutation matrix, $\mathbf{\Pi}$ and $\mathbf{\Sigma}$ are invariant by the transformation $\mathcal{T}$ defined for a generic matrix $\mathbf{A}$ by $\mathcal{T}(\mathbf{A}) = (A_{\tau(j),\tau(i)})_{i,j} = \mathbf{P}_{\tau}^{-1}\mathbf{A}^{\top}\mathbf{P}_{\tau}$. Let $\lambda$ be such that $\mathbf{\Sigma}\mathbf{\lambda} = -g_0\mathbf{j}$. $\mathbf{\Sigma}$ is symmetric so $\mathbf{\Sigma} = \mathbf{P}_{\tau}^{-1}\mathbf{\Sigma}^{\top}\mathbf{P}_{\tau} = \mathbf{P}_{\tau}^{-1}\mathbf{\Sigma}\mathbf{P}_{\tau}$ and $\mathbf{P}_{\tau}\mathbf{\Sigma} = \mathbf{\Sigma}\mathbf{P}_{\tau}$. $\mathbf{P}_{\tau}$ is a permutation matrix so $\mathbf{P}_{\tau}\mathbf{j} = \mathbf{j}$ and

$$-g_0\mathbf{j} = -g_0\mathbf{P}_{\tau}\mathbf{j} = \mathbf{P}_{\tau}\mathbf{\Sigma}\mathbf{\lambda} = \mathbf{\Sigma}\mathbf{P}_{\tau}\mathbf{\lambda} \quad (A.1)$$

hence if $g_0 \neq 0$, $\mathbf{P}_{\tau}\mathbf{\lambda} = -g_0\mathbf{\Sigma}^{-1}\mathbf{j} = \mathbf{\lambda}$. Additionally, $\mathbf{\Pi}\mathbf{\lambda} = \mathbf{0} := \mathbf{E}$ and $\mathbf{\Pi} = \mathbf{P}_{\tau}^{-1}\mathbf{\Pi}^{\top}\mathbf{P}_{\tau} = -\mathbf{P}_{\tau}^{-1}\mathbf{\Pi}\mathbf{P}_{\tau}$ so that

$$\mathbf{P}_{\tau}\mathbf{E} = \mathbf{P}_{\tau}\mathbf{\Pi}\mathbf{\lambda} = -\mathbf{\Pi}\mathbf{P}_{\tau}\mathbf{\lambda} = -\mathbf{\Pi}\mathbf{\lambda} = -\mathbf{E} \quad (A.2)$$

which gives a priori $k$ equations. The equations corresponding to invariant elements of $\tau$ are always satisfied (for some $i$, $\mathbf{E}_i = -\mathbf{E}_i \Rightarrow \mathbf{E}_i = 0$). The others are in even number and half of them can be dropped ($\mathbf{E}_i = -\mathbf{E}_j \Rightarrow \mathbf{E}_j = -\mathbf{E}_i$). The number of remaining equations is the number of transpositions in the decomposition of the permutation $\tau$.

From Theorem 3.3, $\mathbf{\Pi}\mathbf{\Sigma}^{-1}\mathbf{j} = \mathbf{0}$ yields $k-1$ equations for the $k$ unknowns $t_1, \dots, t_k$ corresponding, generically, to a solution space of dimension 1. However, when the sequence $(\sigma_i)_i$ is invariant by a transformation $\tau: i \mapsto a - i - 1$, the number of independent unknowns is reduced. It was just shown that this reduction comes together with the reduction of the number of independent equations. The following proposition indicates that, even in the subgeneric case where symmetries in the flight durations sequence reduce the number of independent unknowns, the solution space dimension is 1.

We now proceed with the proof of Proposition 6.6.

*Proof.* The number of unknowns is determined by the involution $i \mapsto a - i$; the number of equations is determined by the involution $i \mapsto a - 1 - i$.

* If $k = 2r - 1$, both transformations have one invariant element (see Prop. 6.2) and are composed of $r - 1$ transpositions, hence the number of unknowns is $(r - 1) + 1$ and the number of equations is $r - 1$. The space of solutions is of dimension 1.

* If $k = 2r$ and $a$ is even, $i \mapsto a - i - 1$ has no invariant elements and is made of $r$ transpositions, so the number of equations is $r$. The involution $i \mapsto a - i$ has two invariant elements and is made of $r - 1$ transpositions, so the number of unknowns is $r - 1 + 2 = r + 1$. The space of solutions is of dimension 1.

* If $k = 2r$ and $a$ is odd, $i \mapsto a - i$ has two invariant elements and is made of $r - 1$ transpositions, which corresponds to $r - 1$ equations. $i \mapsto a + i - 1$ has no invariant element and is made of $r$ transpositions, yielding $r$ unknowns so, again, a space of solutions of dimension 1.

<figure><img src="image_23.png" /></figure>