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

HAL
open science

# First return time to the contact hyperplane for n-degree-of-freedom vibro-impact systems

Huong Le Thi, Stéphane Junca, Mathias Legrand

► To cite this version:

Huong Le Thi, Stéphane Junca, Mathias Legrand. First return time to the contact hyperplane for n-degree-of-freedom vibro-impact systems. Discrete and Continuous Dynamical Systems - Series B, 2022, 27 (1), pp.1-44. 10.3934/dcdsb.2021031. hal-01957546v2

HAL Id: hal-01957546

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

Submitted on 11 May 2020

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
FIRST RETURN TIME TO THE CONTACT HYPERPLANE FOR
N-DEGREE-OF-FREEDOM VIBRO-IMPACT SYSTEMS

HUONG LE THI

*Department of Mathematics and Informatics, Thang Long University, Hanoi, Vietnam*

STÉPHANE JUNCA

*LJAD Mathematics, Inria & CNRS, Université Côte d'Azur, Nice, France*

MATHIAS LEGRAND

*Department of Mechanical Engineering, McGill University, Montréal, Québec, Canada*

**ABSTRACT.** The paper deals with the dynamics of conservative N-degree-of-freedom vibro-impact systems involving one unilateral contact condition and a linear free flow. Among all possible trajectories, grazing orbits exhibit a contact occurrence with vanishing incoming velocity which generates mathematical difficulties. Such problems are commonly tackled through the definition of a Poincaré section and the attendant First Return Map. It is known that the First Return Time to the Poincaré section features a square-root singularity near grazing. In this work, a non-orthodox yet natural and intrinsic Poincaré section is chosen to revisit the square-root singularity. It is based on the unilateral condition but is not transverse to the grazing orbits. A detailed investigation of the proposed Poincaré section is provided. Higher-order singularities in the First Return Time are exhibited. Also, activation coefficients of the square-root singularity for the First Return Map are defined. For the linear and periodic grazing orbits from which bifurcate nonlinear modes, one of these coefficients is necessarily non-vanishing. The present work is a step towards the stability analysis of grazing orbits, which still stands as an open problem.

CONTENTS

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Main results</td><td>3</td></tr><tr><td>2.1. Poincaré section</td><td>3</td></tr><tr><td>2.2. Square-root singularity</td><td>5</td></tr><tr><td>2.3. Square-root singularity near a grazing linear mode</td><td>7</td></tr><tr><td>3. Domain of definition of the First Return Map</td><td>8</td></tr><tr><td>3.1. Contact times</td><td>8</td></tr><tr><td>3.2. Zero, one or infinitely many?</td><td>10</td></tr><tr><td>3.3. Poincaré section</td><td>12</td></tr><tr><td>4. Implicit function theorem and power-root singularity</td><td>14</td></tr><tr><td>4.1. Square-root singularity</td><td>15</td></tr><tr><td>4.2. Power-root singularity</td><td>18</td></tr><tr><td>4.3. Discontinuous First Return Time</td><td>24</td></tr><tr><td>5. Square-root instability</td><td>25</td></tr><tr><td>5.1. A nonlinear n-dimensional map</td><td>26</td></tr><tr><td>5.2. Two-dimensional maps with critical instability</td><td>28</td></tr><tr><td>6. In the vicinity of a grazing orbit</td><td>30</td></tr><tr><td>6.1. First Return Time</td><td>30</td></tr><tr><td>6.2. Square-root singularity in the vicinity of a grazing linear mode</td><td>32</td></tr><tr><td>7. Conclusion</td><td>34</td></tr><tr><td>Appendix A. One-dof system: on the square-root singularity</td><td>35</td></tr></table>

*Key words and phrases.* non-smooth analysis, grazing, unilateral contact, Poincaré section, First Return Time.
# 1. INTRODUCTION

In the present work, the dynamics of a mechanical system with *N* degrees-of-freedom (dof), one of which being unilaterally constrained, is of interest. An example of such a system is

<figure><img src="image_3.png" /><figcaption>FIGURE 1. A unilaterally constrained *N*-degree-of-freedom chain with *d* > 0.</figcaption></figure>

depicted in Figure 1. The governing equations considered in this work read<fnref n="1" />:

$$ \mathbf{M}\ddot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{r} \qquad (1.1a) $$

$$ \mathbf{u}(0) = \mathbf{u}_0, \quad \dot{\mathbf{u}}(0) = \dot{\mathbf{u}}_0 \qquad (1.1b) $$

$$ u_N(t) \le d, \quad R(t) \le 0, \quad (u_N(t) - d)R(t) = 0 \qquad (1.1c) $$

$$ \dot{\mathbf{u}}^{\top}(t)^{\top} \mathbf{M} \dot{\mathbf{u}}^{\top}(t) + \mathbf{u}^{\top}(t) \mathbf{K} \mathbf{u}(t) = \mathbf{E}(\mathbf{u}(t), \dot{\mathbf{u}}^{\top}(t)) = \mathbf{E}(\mathbf{u}(0), \dot{\mathbf{u}}(0)), \quad \forall t \ge 0 \qquad (1.1d) $$

with

$$ \mathbf{M} = \mathrm{diag}(m_j)_{j=1}^N; \quad \mathbf{K} = (k_{ij})_{i,j=1}^N; \quad \mathbf{u}(t) = (u_j)_{j=1}^N; \quad \mathbf{r}(t) = (0, \dots, 0, R(t)); $$

where $\dot{u}_j$, $\dot{u}_j^+$, $\dot{u}_j^-$ and $\ddot{u}_j$ represent the velocity, the right limit of the velocity, the left limit of the velocity and the acceleration of mass $j$, $j = 1, \dots, N$, respectively. Matrices **M** and **K** are assumed to be symmetric constant positive definite.<fnref n="2" /> Hence, there exists a matrix **P** of **M**-orthogonal eigenmodes which diagonalizes both **M** and **K**, that is $\mathbf{P}^{\top}\mathbf{MP} = \mathbf{I}$ and $\mathbf{P}^{\top}\mathbf{KP} = \mathbf{\Omega}^2 = \mathrm{diag}(\omega_j^2)|_{j=1,...,N}$ where **I** is the $N \times N$ identity matrix, $\omega_j^2$, the eigenfrequencies and $T_j$, the linear periods with $\omega_j T_j = 2\pi$, $j = 1, \dots, N$. Also, condition (1.1c) says that mass *N* is constrained on the right side by a rigid obstacle at a distance *d* > 0 from its equilibrium. There is only one constraint on mass *N*. The other masses are not constrained in any way. The quantity $R(t)$ is the reaction force induced by the obstacle on mass *N* at the time of gap closure. Generally, $R(t)$ is a measure.

The Initial Value Problem (1.1) without (1.1d) is not well-posed: it is known that uniqueness is not guaranteed. To overcome this issue, an impact law is usually incorporated into the formulation. The present work targets non-dissipative dynamics and condition (1.1d) is enforced: the total energy of the system is preserved during the motion. This implies the existence of a perfectly elastic impact law of the form $\dot{u}_N^+ = -e\dot{u}_N^-$ with $e = 1$ where $\dot{u}_N^-$ and $\dot{u}_N^+$ respectively stand for the pre- and post-impact velocities of mass *N*. For the well-posedness of the initial-value problem with constant energy, see [1, 13].

With the above impact law, two types of closing contacts should be addressed: contacts with non-zero pre-impact velocity (or simply “impacts”) and contacts with zero pre-impact velocity. This second category can itself be divided as follows: “grazing contact” if the mass leaves the obstacle immediately after contact is closed, or “sticking contact”<fnref n="3" /> if the mass stays in contact with the obstacle for a finite time interval. These types of contacts are illustrated in Figure 2.

<footnote n="1">¹The fact that $u_N$ is chosen to be unilaterally constrained is not limiting. It could be any other degree-of-freedom through a permutation of indices.</footnote>

<footnote n="2">²A positive definite mass matrix slightly more general than a diagonal one which has a single non-zero entry *NN* on row *N* and column *N* can be considered as well.</footnote>

<footnote n="3">³Note that the term “sticking” does not mean that there is glue on the wall. The contact force $R(t)$ of the wall on mass *N* is still necessarily negative as stated in Equation (1.1c).</footnote>
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Possible types of contacts in unilaterally constrained discrete dynamics.</figcaption></figure>

Recent results on nonlinear modal analysis of discrete vibro-impact systems suggest that the First Return Map (FRM) could be defined on the unilateral constraint hyperplane {$u_N = d$} in the phase-space [9, 18, 8]. This hyperplane gives access to special periodic solutions with closed-form expressions [9, 17] and to reduced-order systems of nonlinear equations together with the companion stability analysis [16, 15]. Yet this Poincaré section has drawbacks as already reported [2, p. 261-262]: it is not transverse to the grazing orbits. This is the main reason why this section is not chosen in [2]. Instead, the authors proposed a First Return Map, called the *discontinuity map*, transverse to the flow but which requires two times, first to reach the contact interface and then, the Poincaré section.

First, we prove that all orbits with a nonlinear behavior intersect the chosen Poincaré section infinitely many times. Second, the square-root singularity is here revisited in a detailed mathematical framework.

The paper is organized as follows: in Section 2, the main results are stated including the definition of the Poincaré section and the square-root singularity of the First Return Time (FRT). Section 3 deals with the exact subset of initial data at the contact constraint $u_N = d$ where the trajectory will come back to the constraint, that is the domain of the FRM. The square-root singularity and higher-order singularities are thoroughly explored in Section 4 in a mathematical framework where the implicit function theorem is used in a degenerate case. The dynamics induced by the square-root singularity is explored in Section 5. Coefficients of the asymptotic expansion of the FRM activating the square-root singularity are identified in Section 6. These coefficients are explicitly computed for grazing linear modes (GLM). The main results of this article are then summarized in the conclusion.

## 2. MAIN RESULTS

The Poincaré section is first defined in Section 2.1. The theorems on the square-root singularity which emerges in the vicinity of grazing orbits are provided in Section 2.2. The last section addresses a result on the square-root singularity near a periodic solution with a sole grazing occurrence for all times.

2.1. **Poincaré section.** The Poincaré section corresponds to the domain of definition of the First Return Map. In nonsmooth analysis, where the vector field governing the dynamics is commonly piecewise smooth only, the hyperplane $\mathcal{H} = \{[\mathbf{u}^T, \dot{\mathbf{u}}^{-T}]^T \in \mathbb{R}^{2N}, u_N = d\}$ of the phase-space defined by $u_N = d$ is a natural choice for the Poincaré section when targeting trajectories with non-vanishing pre-impact velocities since they hit the section transversally [9]. However, for grazing trajectories, the transversality condition is lost. A discontinuity mapping on another suitable section for which the transversality condition is retrieved shall be used instead [2, 7, 11]. To study the stability of grazing periodic orbits, the approach proposed in [2, 7, 11] is quite natural and efficient. To describe the behavior near grazing contact, our approach is more direct than the previous one. We obtain new insights on the square-root singularity and on the complex behavior near a grazing contact in the state space. In particular, we take advantage of the simple linear free dynamics outside the contact.
In the present work, the Poincaré section is still a subset of the hyperplane $\mathcal{H}$. It is the simplest cross-section to describe the dynamics with only two phases: contact dynamics and free-flight dynamics. This does not have the adverse effect of introducing a second free-flight dynamics as the discontinuity mapping does. Nevertheless, our critical choice necessitates a very careful description of the domain of definition of the FRM. This is stated in the following theorem that categorizes the initial data generating orbits which will always come back to the section.

**Proposition 2.1 — Finite sticking duration.** When $d$ is positive, the sticking phase of a solution to (1.1) is of finite duration.

This is proven in [15] for a two-dof vibro-impact system. A general proof for the $N$-dof system is given in Section 3.1 of the present work. This property is used to show that a solution to (1.1) has zero, one or infinitely many closing contacts with the hyperplane $\mathcal{H}$. Before stating this result, the following assumption is needed.

**Assumption 2.1 — No internal resonances.** The linear frequencies of system (1.1) are $\mathbb{Z}$ independent which means

$$ \sum_{i=1}^{N} k_i \omega_i = 0 \text{ and } k_i \in \mathbb{Z} \quad \Rightarrow \quad k_i = 0, \quad \forall i = 1, \dots, N. \qquad (2.1) $$

Through this assumption, internal resonances are discarded in the current work.

**Theorem 2.2 — Zero, one or infinitely many closing contacts with $\mathcal{H}$.** Let $\mathbf{u}(t)$ be a solution to (1.1). Under Assumption 2.1 of no internal resonance, the solution is such that

**Case 1 — linear solution:** mass $N$ never hits $\mathcal{H}$, i.e. $u_N(t) < d$ for all $t$.

**Case 2 — linear solution:** mass $N$ experiences only one closing contact, i.e. there exists $t_0$ such that $u_N(t_0) = d$ and $u_N(t) < d$ for all $t \neq t_0$.

**Case 3 — nonlinear solution:** mass $N$ experiences a countably infinite number of isolated closing contacts with $\mathcal{H}$.

This theorem is proven in Section 3.2. The affine space $\mathcal{H}$ is of dimension $2N - 1$. It is divided into three disjoint subsets:

$$ \mathcal{H}^{-} = \{[\mathbf{u}^{\mathrm{T}}, \dot{\mathbf{u}}^{-\mathrm{T}}]^{\mathrm{T}} \in \mathbb{R}^{2N}, u_{N} = d \text{ and } \dot{u}_{N}^{-} > 0\}, \qquad (2.2) $$

$$ \mathcal{H}^{+} = \{[\mathbf{u}^{\mathrm{T}}, \dot{\mathbf{u}}^{+\mathrm{T}}]^{\mathrm{T}} \in \mathbb{R}^{2N}, u_{N} = d \text{ and } \dot{u}_{N}^{+} < 0\}, \qquad (2.3) $$

$$ \mathcal{H}^{0} = \{[\mathbf{u}^{\mathrm{T}}, \dot{\mathbf{u}}^{\mathrm{T}}]^{\mathrm{T}} \in \mathbb{R}^{2N}, u_{N} = d \text{ and } \dot{u}_{N} = 0\}. \qquad (2.4) $$

By a careful investigation of the free dynamics, it can be seen that the solution which features a contact with non-zero pre-impact velocity will always experience a later (and a previous in negative time) closing contact. Accordingly, the problem of the existence of a subsequent (or previous) closing contact emerges only on $\mathcal{H}^0$. Theorem 2.2 implies that $\mathcal{H}^0$ is the union of the two subsets

$$ \mathcal{H}_{\infty}^{0} = \{[\mathbf{u}^{\mathrm{T}}, \dot{\mathbf{u}}^{\mathrm{T}}]^{\mathrm{T}} \in \mathcal{H}^{0} \text{ with an infinite number of contacts}\}, \qquad (2.5) $$

$$ \mathcal{H}_{1}^{0} = \{[\mathbf{u}^{\mathrm{T}}, \dot{\mathbf{u}}^{\mathrm{T}}]^{\mathrm{T}} \in \mathcal{H}^{0} \text{ with a single grazing contact}\}. \qquad (2.6) $$

Since $\mathcal{H}_1^0$ contains solutions with only one grazing contact, it does not belong to the Poincaré section.

**Definition 2.1 — Poincaré section.** The Poincaré section $\mathcal{H}_\rho \subset \mathcal{H}$ is formed by the union of the set of the states with non-zero velocity contacts and the set of the states with zero velocity contact that gives rise to an infinite number of closing contacts, that is

$$ \mathcal{H}_{\rho} = \mathcal{H}^{-} \cup \mathcal{H}_{\infty}^{0}. \qquad (2.7) $$
**Remark 2.1.** There are two options for the choice of the Poincaré section in $\mathcal{H}$ depending on whether one wants to start right before or right after the contact occurrence. The former gives the Poincaré section defined as above, and the latter yields $\mathcal{H}_{\rho}^{+} = \mathcal{H}^{+} \cup \mathcal{H}_{\infty}^{0}$ as the Poincaré section.

The set $\mathcal{H}_{\infty}^{0}$ can also be split into two subsets: $\mathcal{H}_{S}^{0}$ including all the initial data belonging to $\mathcal{H}_{\infty}^{0}$ such that the solution starts by a sticking contact and $\mathcal{H}_{G}^{0}$ of initial data such that the solution starts by a grazing contact:

$$ \mathcal{H}_{\infty}^{0} = \mathcal{H}_{S}^{0} \cup \mathcal{H}_{G}^{0}. \qquad (2.8) $$

**2.2. Square-root singularity.** The First Return Time is known to be generically an analytic function of the initial data. Let $W_0 = [u_0, \dot{u}_0]^\top \in \mathcal{H}^-$. If the first return to $\mathcal{H}$, named $W_1$, belongs to $\mathcal{H}^-$ then the FRM is analytic near $W_0$ [2, 9]. However, if $W_1 \in \mathcal{H}^0$, then there is a grazing contact and it is known that a square-root singularity appears [7, 11].

By definition of $\mathcal{H}_{\rho}$ and through Theorem 2.2, there exists a time, called “First Return Time”, at which the orbit emanating from $W \in \mathcal{H}_{\rho}$ comes back to $\mathcal{H}_{\rho}$.

**Definition 2.2 — First Return Time.** Let $\mathbf{u}(t)$ be a solution to (1.1) with the initial data $W = (W_i)_{i=1}^{2N} \in \mathcal{H}_{\rho}$ at the initial time $t=0$, i.e. $W_N = u_N(0) = \mathbf{e}_N^\top \mathbf{u}(0) = d$. The First Return Time $T = T(W) > 0$ is defined by

$$ T(W) = \begin{cases} \min\{t > 0 : u_N(t) = d\} & \text{if there is no sticking phase at } t = 0 \\ \min\{t > \tau(W) : u_N(t) = d\} & \text{if there is sticking phase at } t = 0 \end{cases} \qquad (2.9a) $$

$$ T(W) = \begin{cases} \min\{t > \tau(W) : u_N(t) = d\} & \text{if there is sticking phase at } t = 0 \\ \min\{t > 0 : u_N(t) = d\} & \text{if there is no sticking phase at } t = 0 \end{cases} \qquad (2.9b) $$

where $\tau(W)$ is the sticking duration.

This definition can be made more concise by saying that $T(W) = \min\{t > \tau(W) : u_N(t) = d\}$ with the convention that $\tau(W) = 0$ if there is no sticking at $t = 0$. When a sticking phase occurs at $t = 0$, Proposition 2.1 ensures that the First Return Time is well defined since the duration of the sticking phase is finite.

Let $\mathbf{u}(t, W)$ be the solution associated with the initial data $W \in \mathcal{H}_{\rho}$. By the definition of $\mathcal{H}_{\rho}$, $W_N = d$, hence $W$ is viewed as a vector of $2N-1$ variables. The displacement of mass $N$ is $u_N(t, W) = \mathbf{e}_N^\top \mathbf{u}(t, W)$ which is a function of time variable $t$ and the initial data $W$. Consider the smooth function

$$ \Phi(t, W) = \mathbf{e}_{2N}^{\top} \mathbf{R}(t) SW \qquad (2.10) $$

defined for all $t \in \mathbb{R}$ and for all $W \in \mathbb{R}^{2N}$ where $(1) \mathbf{S} = \text{diag}(1, \dots, 1, -1)$ is a $2N \times 2N$ diagonal matrix with last entry $-1$ to reflect the impact law and (2) the operator

$$ \mathbf{R}(t) = \begin{bmatrix} \mathbf{P} \cos(t\mathbf{\Omega})\mathbf{P}^{-1} & \mathbf{P}\mathbf{\Omega}^{-1} \sin(t\mathbf{\Omega})\mathbf{P}^{-1} \\ -\mathbf{P}\mathbf{\Omega} \sin(t\mathbf{\Omega})\mathbf{P}^{-1} & \mathbf{P} \cos(t\mathbf{\Omega})\mathbf{P}^{-1} \end{bmatrix} \qquad (2.11) $$

describes the free-flight dynamics. The dynamics strictly between two successive closing contacts is smooth. The function $\Phi(t, W)$ coincides with $u_N(t, W)$ as long as $u_N(t, W) < d$.

**Assumption 2.3 — Non-vanishing acceleration.** Is considered $W_0$ the initial data leading to an orbit which has a grazing contact at the first time $T_0$ with the condition $\partial_t^2 \Phi(T_0, W_0) \neq 0$.<fnref n="4" />

This is an important assumption which activates the square-root singularity near a grazing orbit. Otherwise, a stronger singularity will emerge.

**Assumption 2.3.** $\nabla_W \Phi(T_0, W_0) \neq 0_{\mathbb{R}^{2N}}$.

This assumption is generally satisfied and means that there exists $1 \le k \le 2N$, $k \ne N$ such that $\partial_{W_k} \Phi(T_0, W_0) \ne 0$. Without it, the square root singularity in the First Return Time is not expected to arise. Also, this assumption 2.3 is not verified for $N = 1$, as explained in Appendix A.

<footnote n="4">⁴In this paper, $\partial_t^n = (\partial_t)^n$ is the n-th order partial derivative with respect to $\bullet$.</footnote>
However, for $N \ge 2$, it always holds true for an initial data $W_0$ corresponding to a grazing linear mode, as shown in Lemma 6.1.

With Assumptions 2.2 and 2.3, the scalar

$$ \gamma_k = - \frac{\partial_t^2 \Phi(T_0, W_0)}{2 \partial_{W_k} \Phi(T_0, W_0)} \qquad (2.12) $$

is well defined and does not vanish. The expression of the First Return Time near a grazing contact is then given in Theorem 2.3. Throughout the paper, the expression “near a grazing contact” means “with initial data $W$ near $W_0$”, and the corresponding solution is investigated near time $T_0$. The implicit function theorem is applied on the smooth function $\Phi$ in place of the nonsmooth function $u_N$. It is easier to differentiate the smooth function $\Phi$ than the non-smooth function $u_N$.

**Theorem 2.3 — Square-root singularity near a grazing contact.** Let $W_0 = (W_{0i})_{i=1}^{2N} \in \mathcal{H}^-$ be the initial data generating an orbit with a grazing contact at the First Return Time $T_0 = T(W_0)$. If Assumption 2.3 applies, then there exists a component $W_k$ of $W$ such that $\partial_{W_k} \Phi(T_0, W_0) \ne 0$. Let $\underline{W} \in \mathbb{R}^{2N-1}$ be the reduced vector obtained from $W$ by removing $W_k$. If Assumptions 2.1 and 2.2 hold then:

(1) There exist two neighborhoods $V_{T_0}$ and $V_{W_0} \subset \mathcal{H}_\phi$ of $T_0$ and $W_0$, respectively as well as two smooth scalar functions $\eta$ and $\alpha$ defined on $V_{\underline{W}_0} = \{\underline{W}, W \in V_{\underline{W}_0}\}$ such that the set

$$ S_c = \{(t, W) \in \mathbb{R} \times \mathbb{R}^{2N}, \Phi(t, W) = d \text{ and } \partial_t \Phi(t, W) = 0\} \qquad (2.13) $$

where the square-root singularity will emerge, is locally parametrized as follows:

$$ S_c \cap \{V_{T_0} \times V_{W_0}\} = \{(\eta(\underline{W}), \alpha(\underline{W}), \underline{W}), \underline{W} \in V_{\underline{W}_0}\}, \qquad (2.14) $$

where $\eta(\underline{W}_0) = T_0$ and $\alpha(\underline{W}_0) = W_{0k}$.

(2) Let $s_k = \text{sign } \gamma_k$ and the set

$$ \mathcal{B}_k = \{W \in V_{\underline{W}_0}, s_k(W_k - \alpha(\underline{W})) \ge 0\} \qquad (2.15) $$

Let $\sigma = \text{sign}(\dot{u}_N^-(T_0, \underline{W}_0))$. There exists a smooth function $\psi$ such that $\psi(0, \underline{W}_0) = 0$, $\partial_1 \psi(0, \underline{W}_0) = 1/\sqrt{|\gamma_k|}$ and for all $W \in \mathcal{B}_k$, the First Return Time $T$ is given by

$$ T(W) = \eta(\underline{W}) + \psi(\sigma \sqrt{s_k(W_k - \alpha(\underline{W})), \underline{W}}). \qquad (2.16) $$

Moreover, the set $\{W \in \mathcal{H}_\phi, (T(W), W) \in V_{T_0} \times V_{\underline{W}_0}\}$ is exactly $\mathcal{B}_k$.

A mathematical proof of this theorem is provided in Section 6.1. Take note that the square-root dependence on the initial data only appears on the set $\mathcal{B}_k$ which is simply the region above or below the hypersurface $W_k = \alpha(\underline{W})$. The issue induced by this square-root singularity is already known for a Poincaré section transverse to $\mathcal{H}$ [7, 2, 11, 19]. These existing results are compared to our FRT in Section 4.3.

The presence of the square-root singularity is the consequence of the vanishing incoming velocity $\dot{u}_N^-(T_0) = 0$ and the non-vanishing acceleration $\ddot{u}_N^-$ at $T_0$. Contact occurrences with both vanishing velocity and vanishing acceleration, that is $\dot{u}_N^-(T_0) = \ddot{u}_N^-(T_0) = 0$, but with $\ddot{u}_N^-(T_0) \ne 0$ are expected to generate a cube-root singularity. More generally, if $u_N^{(\ell)-}(T_0) = 0, 0 < \ell < n$, and $u_N^{(n)-}(T_0) \ne 0$, then a $n^{\text{th}}$-root singularity is expected, where $u_N^{(\ell)-}$ refers to the $\ell^{\text{th}}$ left time-derivative of $u_N$. This is proved in Section 4.2 under generic assumptions and $n$ is shown to be bounded by $2N$ where $N$ is the number of degrees-of-freedom.

**Remark 2.2 — Discontinuous First Return Time.** Since $\mathcal{B}_k$ is a one-sided set bounded by an hypersurface, the set of initial data with their First Return Time near $T_0$ is not a neighborhood of $W_0$. This means that in the vicinity of $W_0$, there are initial data with a First Return Time far from $T_0$. In other words, the First Return Time $T(W)$ is discontinuous at $W_0$, see Section 4.3.
This geometric singularity is induced by our choice of a non-transverse Poincaré section. The dynamical interpretation of this discontinuity is discussed in Section 4.3.

Theorem 2.3 also has to be generalized for initial data $W_0 \in \mathcal{H}_G^0$, for instance the initial data of a grazing linear mode. The neighborhood $V_{W_0}$ has to be replaced by the half neighborhood $V_{W_0}^+ = \{W \in V_{W_0}, \mathbf{e}_{2N}^\top W \ge 0\}$ where the initial velocity is non-positive and the set $\mathcal{B}_k$ is replaced by the smaller set

$$ \mathcal{B}_k^+ = \mathcal{B}_k \cap V_{W_0}^+ = \{W \in V_{W_0}, \gamma_k(W_k - \alpha(\underline{W})) \ge 0, W_{2N} \ge 0\}. \quad (2.17) $$

The following theorem is stated with the notations of Theorem 2.3.

**Theorem 2.4** — Square-root singularity for two successive grazing contacts. Assume that $W_0 \in \mathcal{H}_G^0$, $\ddot{u}_N^+(0, W_0) < 0$ and $W(T_0) \in \mathcal{H}^0$ where $T_0 = T(W_0)$. If Assumption 2.3 applies, there exists $k \in \{1, \dots, 2N\}$ and $k \ne N$ such that $\partial_{W_k} \Phi(T_0, W_0) \ne 0$ and $s_k = \text{sign}(\gamma_k)$ together with $\sigma = \text{sign}(\ddot{u}_N^-(T_0, W_0))$. If Assumptions 2.1 and 2.2 hold, then there exist two neighborhoods $V_{T_0}$ and $V_{W_0} \subset \mathcal{H}_\rho$ of $T_0$ and $W_0$, respectively as well as two smooth scalar functions $\eta$ and $\alpha$ defined on $V_{W_0}$ containing $\underline{W}_0$ where $\eta(\underline{W}_0) = T_0$ and $\alpha(\underline{W}_0) = W_{0k}$. There also exists a smooth function $\psi$ such that $\psi(0, \underline{W}_0) = 0$, $\partial_1 \psi(0, \underline{W}_0) = |\gamma_k|^{-1/2}$ and for all $W \in \mathcal{B}_k^+$, the First Return Time is

$$ T(W) = \eta(\underline{W}) + \psi(\sigma \sqrt{s_k(W_k - \alpha(\underline{W})), \underline{W}}). \quad (2.18) $$

Moreover, the set $\{W \in \mathcal{H}_\rho, (T(W), W) \in V_{T_0} \times V_{W_0}^+\}$ is exactly $\mathcal{B}_k^+$.

The difference between Theorems 2.3 and 2.4 is that formula (2.16) of the First Return Time in the former is defined on $\mathcal{B}_k$ whereas in the latter it is defined on a smaller subset $\mathcal{B}_k^+$. This is due to the fact that $\mathcal{H}_G^0$ lies on the boundary of the cross-section $\mathcal{H}_\rho$, and not all perturbations of $W_0$ are thus admissible. The condition $\ddot{u}_N^+(0, W_0) < 0$ guarantees that solutions with initial data near $W_0$ and with a non-negative velocity for the last mass have a First Return Time near $T_0$. The case $W_0 \in \mathcal{H}_S^0$ can also be considered and can add an additional singularity due to the sticking phase [15].

### 2.3. Square-root singularity near a grazing linear mode.
The square-root singularity near the periodic solutions with one grazing contact per period is addressed. For a $N$-dof system without internal resonances, there are $N$ such periodic solutions which are called *grazing linear modes* (GLM). It is known that many invariant submanifolds of periodic solutions with $k$ Impact-Per-Period ($k$-IPP) [9, 18, 16] might emerge in the vicinity of GLM.

We define the $j^{th}$ grazing linear mode as a periodic trajectory $\mathbf{u}$ associated to the $j^{th}$ linear mode which satisfies $\max_{t \in \mathbb{R}} u_N(t) = d$, i.e. the contacts are at most of grazing type. An essential tool to investigate the dynamics near a grazing linear mode is the First Return Map [2]. This map is well defined on $\mathcal{H}_\rho$.

**Definition 2.3** — First Return Map. Suppose $W \in \mathcal{H}_\rho$ and $T = T(W) > 0$ is the First Return Time to $\mathcal{H}_\rho$ of the orbit emanating from $W$. The map which associates points in $\mathcal{H}_\rho$ to their first return images to $\mathcal{H}_\rho$ is called the First Return Map $\mathcal{F}$. To be more precise $\mathcal{F}: \mathcal{H}_\rho \to \mathcal{H}_\rho$ with

$$ \mathcal{F}(W) = \begin{cases} \mathbf{R}(T(W))\mathbf{S}\mathbf{W} & \text{if } W \in \mathcal{H}^- \cup \mathcal{H}_G^0 \\ \mathbf{R}(s(\mathbf{U}(\tau(W)))\mathbf{U}(\tau(W)) & \text{if } W \in \mathcal{H}_S^0 \end{cases} \quad (2.19a) $$

$$ \mathcal{F}(W) = \begin{cases} \mathbf{R}(s(\mathbf{U}(\tau(W)))\mathbf{U}(\tau(W)) & \text{if } W \in \mathcal{H}_S^0 \\ \mathbf{R}(T(W))\mathbf{S}\mathbf{W} & \text{if } W \in \mathcal{H}^- \cup \mathcal{H}_G^0 \end{cases} \quad (2.19b) $$

where the matrix $\mathbf{S}$ describes the impact law; $\mathbf{U}(t) = [\mathbf{u}(t), \dot{\mathbf{u}}(t)]$ is the state of the system at time $t$; $\tau$ and $s$ are the duration of the sticking phase and of the free-flight phase, respectively.

Let us investigate formula (2.19b). If there is a sticking phase, i.e. $W \in \mathcal{H}_S^0$, then $\tau = \tau(W)$ is the sticking duration and $s = T(W) - \tau(W)$ denotes the duration of the free-flight after the sticking phase until the next contact. The state of the system at the end of the sticking phase is called $\mathbf{U}(\tau(W))$. If $W \in \mathcal{H}^- \cup \mathcal{H}_G^0$, i.e. there is no sticking at $t = 0$, then $\tau$ is assumed to be 0,
hence $s \equiv T$ and $U(\tau(W)) \equiv W$. In other words, formula (2.19b) is valid for all the cases, with or without sticking. An explicit formula of $\tau$ for a two-dof system is exposed in [15].

In (2.16), the square-root singularity shows that periodic orbits with grazing contacts yield intricate dynamics such as instability of the periodic orbits or grazing bifurcations [7, 11]. This square-root term may produce the so-called square-root singularity of the First Return Map via the coefficients $C_k$, as defined below. In the particular framework of Section 5, instability of the grazing linear modes is expected. However, this is not true in the one-dof case, see Appendix A.

**Definition 2.4** — **Square-root singularity coefficients.** Suppose that $W_0 \in \mathcal{H}_\rho$ generates an orbit with the first contact at $T_0$ of grazing type, i.e. $u_N(T_0) = d$ and $\dot{u}_N(T_0) = 0$. Under Assumptions 2.1 of no internal resonance and 2.3 for each $k \in \{1, \dots, 2N\}$ such that $\partial_{W_k} \Phi(T_0, W_0) \neq 0$, the square-root singularity coefficient $C_k$ is defined by

$$C_k = \begin{cases} \sqrt{|\gamma_k|} \mathbf{e}_k^\top \mathbf{P} \dot{\mathbf{q}}(T_0) & \text{if } 1 \le k < N, \\ -\sqrt{|\gamma_k|} \mathbf{e}_{k-N}^\top \mathbf{M} \mathbf{P} \mathbf{\Omega}^2 \mathbf{q}(T_0) & \text{if } N < k \le 2N \end{cases} \quad (2.20a)$$

$$C_k = \begin{cases} -\sqrt{|\gamma_k|} \mathbf{e}_{k-N}^\top \mathbf{M} \mathbf{P} \mathbf{\Omega}^2 \mathbf{q}(T_0) & \text{if } N < k \le 2N \end{cases} \quad (2.20b)$$

where $[\mathbf{q}, \dot{\mathbf{q}}]$ are the modal coordinates defined by the transformation $\mathbf{u} = \mathbf{P}\mathbf{q}$.

The square-root singularity for the First Return Map is then defined as follows.

**Definition 2.5** — **Square-root singularity for the First Return Map.** System (1.1) is said to feature square-root singularity for the First Return Map near a grazing periodic solution associated to the initial condition $W_0$ if there exists at least a coefficient $C_k$, $1 \le k \le 2N$, $k \ne N$ which does not vanish.

Near a grazing linear mode, the square-root singularity for the First Return Map is shown to exist under a generic condition as stated in Theorem 2.5 proven in Section 6.2.

**Theorem 2.5** — **Square-root singularity for the First Return Map near grazing linear modes.** Consider the $j^{th}$ grazing linear mode of (1.1) associated to the initial state $W_0 \in \mathcal{H}_\rho$ and $T_j$, its period. Assumptions 2.1, 2.2 and 2.3 hold. If there exists an index $N < i \le 2N$ such that $\partial_{W_i} u_N(T_j, W_0) \neq 0$ and $P_{ij} \neq 0$, then there exists a square-root singularity for the First Return Map near the $j^{th}$ grazing linear mode.

Instead, if $C_k = 0$ for all $k$ then the square-root singularity for the First Return Map is not activated and the dynamics near the periodic orbit is expected to be similar to that of the smooth dynamics. Roughly speaking, the FRM annihilates the singularity in the same way as squaring removes a square-root singularity.

Theorem 2.5 gives a hint on how to explore the instability of the grazing linear modes. For instance, if the dynamics starts in $B_k$ and stays in $B_k$, then the framework of Section 5 follows where the instability of the associated First Return Map fixed-point is explained.

The condition $P_{ij} \neq 0$ comes from the definition of $C_i$ with $N < i \le 2N$. Despite the fact that $P_{ij} \neq 0$ is a generic property, it might be violated for a chain of masses with $N > 2$. Hence, in order to have such a $C_i \neq 0$, a condition on $P_{ij}$ is required. When $N = 2$ [15], it is shown that $P_{ij} \neq 0$ for all $i, j = 1, 2$ and there is always a square-root singularity for the First Return Map near the grazing linear modes.

### 3. DOMAIN OF DEFINITION OF THE FIRST RETURN MAP

This section details the Poincaré section in a comprehensive manner.

3.1. **Contact times.** The times at which contact occurs are defined and then categorized. In particular, if a sticking phase starts, then it is of finite duration. Moreover, the total energy of the system is conserved during the sticking phase. These features will be used to show that there is a countably infinite number of closing contacts after a sticking phase, further details are found in Section 3.2.
**Definition 3.1** — Contact time. Assume that $\mathbf{u}(t)$ is the solution to system (1.1); $T$ is a contact time if $u_N(T) = d$ and there exists $\delta > 0$ such that $u_N(T - t) < d$ for all $0 < t < \delta$.

In other words, a contact time $T$ is the time at which mass $N$ touches the rigid obstacle after a free-flight phase. Invoking [1], the contact time is well defined in the conservative case which is without chattering.

Contact times are classified into three categories:

(1) It is a contact with non-zero pre-velocity if $u_N(T) = d$ and $\dot{u}_N(T) > 0$.

(2) It is a contact with zero pre-velocity if $u_N(T) = d$ and $\dot{u}_N(T) = 0$, with two possibilities:

(a) a grazing contact if the mass leaves the obstacle right after the contact time;

(b) a sticking contact if the mass stays in contact with the obstacle.

In this paper, the term “closing contact” indifferently refers to either an impulsive impact or a grazing contact or the beginning of a sticking contact phase. The sticking system dictates the dynamics during the sticking phase [2]. Since the last mass rests against the obstacle, the system of interest “loses” one degree-of-freedom. The sticking system complemented by the initial data at the beginning of a sticking phase reads

$$ \underline{\mathbf{M}}\ddot{\mathbf{u}} + \underline{\mathbf{K}}\mathbf{u} = -d \underline{\mathbf{I}}_N^\top \qquad (3.1) $$

$$ m_N \ddot{u}_N = 0, \quad u_N(0) = d, \quad \dot{u}_N(0) = 0 \qquad (3.2) $$

where $\underline{\mathbf{M}}$ and $\underline{\mathbf{K}}$ are the mass and stiffness matrices after removing the last row and last column, $\underline{\mathbf{u}}$ is a $N - 1$ vector solution to the sticking dynamics, $\underline{\mathbf{I}}_N^\top$ is the last column of $\underline{\mathbf{K}}$ where the last entry $k_{NN}$ has been removed.

Proposition 2.1 can now be proven.

*Proof.* Let us proceed by contradiction and assume that the sticking phase never ends. The new equilibrium $\underline{\mathbf{u}}_e$ of the sticking system satisfies

$$ \underline{\mathbf{K}}\underline{\mathbf{u}}_e = -d \underline{\mathbf{I}}_N^\top \text{ that is } \underline{\mathbf{u}}_e = -d \underline{\mathbf{K}}^{-1}\underline{\mathbf{I}}_N^\top. \qquad (3.3) $$

During the sticking phase, the last equation reduces to $\ddot{u}_N = 0$. Indeed, the last equation says $R(t) = F(t)$ since $m_N\ddot{u}_N + \underline{\mathbf{I}}_N\underline{\mathbf{u}} + k_{NN}d = R(t)$, or $m_N\ddot{u}_N = -F(t) + R(t)$, with

$$ F(t) = \underline{\mathbf{I}}_N \underline{\mathbf{u}}(t) + k_{NN} d \le 0. \qquad (3.4) $$

It should be understood that $F(t)$ cannot be positive during the sticking phase. Otherwise, there exists $t_0$ during the sticking phase such that $F(t_0) > 0$. By continuity of the function $F$, it is strictly positive on an open interval including $t_0$. By integrating the acceleration, $u_N(t_0) < d$ which contradicts that $u_N(t) = d$ during the sticking phase.

The solution $\underline{\mathbf{u}}$ of the sticking system is quasi-periodic and continuous and its mean value $\langle \underline{\mathbf{u}} \rangle$ is the equilibrium $\underline{\mathbf{u}}_e$. Therefore, the mean value of the scalar quasi-periodic function $F$ is $\langle F \rangle = \underline{\mathbf{I}}_N \langle \underline{\mathbf{u}} \rangle + k_{NN}d = \underline{\mathbf{I}}_N \underline{\mathbf{u}}_e + k_{NN}d$. Since $F$ is continuous, there exists $t_0 > 0$ such that $F(t_0) = \langle F \rangle = -d \underline{\mathbf{I}}_N \underline{\mathbf{K}}^{-1} \underline{\mathbf{I}}_N^\top + k_{NN}d = d \underline{\mathbf{X}}^\top \underline{\mathbf{K}} \underline{\mathbf{X}}$ with $\underline{\mathbf{X}} = [\underline{\mathbf{K}}^{-1} \underline{\mathbf{I}}_N^\top, -1]^\top$. $F(t_0)$ is positive because of the condition $d > 0$ and the positive definiteness of $\underline{\mathbf{K}}$, which contradicts (3.4). Therefore, the sticking duration is finite. $\square$

The conservation of energy during the sticking phase is not necessarily obvious since the sticking and non-sticking systems are different.

**Lemma 3.1** — Energy during sticking phase. The solution to (3.1)-(3.2) preserves the total energy $E$ defined in (1.1d).

*Proof.* Assume that $t = 0$ is the beginning of a sticking phase and $t = \tau$, the end. During this sticking phase on the interval $[0; \tau]$, the governing equations become $\underline{\mathbf{M}}\ddot{\underline{\mathbf{u}}} + \underline{\mathbf{K}}\underline{\mathbf{u}} = \underline{\mathbf{0}}$, with
$\bar{\mathbf{u}} = \mathbf{u} - \mathbf{u}_e$ where $\mathbf{u}_e = -d \mathbf{K}^{-1} \mathbf{I}_N^\top$ is the new equilibrium of the sticking system. This sticking system conserves the following energy:

$$ \bar{\mathbf{E}}(t) = \dot{\bar{\mathbf{u}}}^\top(t) \mathbf{M} \dot{\bar{\mathbf{u}}}(t) + \bar{\mathbf{u}}^\top(t) \mathbf{K} \bar{\mathbf{u}}(t) = \bar{\mathbf{E}}(0). \quad (3.5) $$

Moreover, since $\mathbf{u} = \bar{\mathbf{u}} + \mathbf{u}_e$, an easy manipulation yields

$$ \underline{\mathbf{E}}(t) = \dot{\underline{\mathbf{u}}}^\top(t) \underline{\mathbf{M}} \dot{\underline{\mathbf{u}}}(t) + \underline{\mathbf{u}}^\top(t) \underline{\mathbf{K}} \underline{\mathbf{u}}(t) = \bar{\mathbf{E}}(t) + \underline{\mathbf{u}}_e^\top \underline{\mathbf{K}} \underline{\mathbf{u}}_e + 2 \bar{\mathbf{u}}^\top(t) \underline{\mathbf{K}} \underline{\mathbf{u}}_e. \quad (3.6) $$

In particular,

$$ \underline{\mathbf{E}}(0) = \bar{\mathbf{E}}(0) + \underline{\mathbf{u}}_e^\top \underline{\mathbf{K}} \underline{\mathbf{u}}_e + 2 \bar{\mathbf{u}}^\top(0) \underline{\mathbf{K}} \underline{\mathbf{u}}_e \quad (3.7) $$

and the total energy of the system can now be calculated. Since $\dot{\mathbf{u}}$ is continuous during a sticking phase including its beginning and end, the exponents $\pm$ are dropped, and

$$ \begin{aligned} \mathbf{E}(t) &= \dot{\mathbf{u}}^\top(t) \mathbf{M} \dot{\mathbf{u}}(t) + \mathbf{u}^\top(t) \mathbf{K} \mathbf{u}(t) \\ &= \dot{\mathbf{u}}^\top(t) \underline{\mathbf{M}} \dot{\mathbf{u}}(t) + \underline{\mathbf{u}}^\top(t) \underline{\mathbf{K}} \underline{\mathbf{u}}(t) + k_{NN} u_N^2(t) + 2u_N(t) \underline{\mathbf{I}}_N \underline{\mathbf{u}}(t) \\ &= \underline{\mathbf{E}}(t) + k_{NN} d^2 + 2d \underline{\mathbf{I}}_N \underline{\mathbf{u}}(t) \\ &= \bar{\mathbf{E}}(t) + \underline{\mathbf{u}}_e^\top \underline{\mathbf{K}} \underline{\mathbf{u}}_e + 2 \bar{\mathbf{u}}^\top(t) \underline{\mathbf{K}} \underline{\mathbf{u}}_e + k_{NN} d^2 + 2d \underline{\mathbf{I}}_N (\bar{\mathbf{u}}(t) + \underline{\mathbf{u}}_e) \\ &= \bar{\mathbf{E}}(t) + \underline{\mathbf{u}}_e^\top \underline{\mathbf{K}} \underline{\mathbf{u}}_e - 2d \bar{\mathbf{u}}^\top(t) \underline{\mathbf{I}}_N^\top + k_{NN} d^2 + 2d \underline{\mathbf{I}}_N (\bar{\mathbf{u}}(t) + \underline{\mathbf{u}}_e) \\ &= \bar{\mathbf{E}}(t) + \underline{\mathbf{u}}_e^\top \underline{\mathbf{K}} \underline{\mathbf{u}}_e + k_{NN} d^2 + 2d \underline{\mathbf{I}}_N \underline{\mathbf{u}}_e. \end{aligned} \quad (3.8) $$

Similarly,

$$ \begin{aligned} \mathbf{E}(0) &= \underline{\mathbf{E}}(0) + k_{NN} d^2 + 2d \underline{\mathbf{I}}_N \underline{\mathbf{u}}(0) \\ &= \bar{\mathbf{E}}(0) + \underline{\mathbf{u}}_e^\top \underline{\mathbf{K}} \underline{\mathbf{u}}_e + 2 \bar{\mathbf{u}}^\top(0) \underline{\mathbf{K}} \underline{\mathbf{u}}_e + k_{NN} d^2 + 2d \underline{\mathbf{I}}_N \underline{\mathbf{u}}(0) \\ &= \bar{\mathbf{E}}(0) + \underline{\mathbf{u}}_e^\top \underline{\mathbf{K}} \underline{\mathbf{u}}_e + k_{NN} d^2 + 2d \underline{\mathbf{I}}_N \underline{\mathbf{u}}_e \end{aligned} \quad (3.9) $$

It follows that $\mathbf{E}(t) = \mathbf{E}(0)$, i.e. the total energy of the system is constant during the sticking phase. $\square$

**3.2. Zero, one or infinitely many?** This section discusses the number of closing contacts of any trajectory. The main results are stated in Theorem 2.2. For the system of interest, the number of closing contacts can only be either 0, 1 or countably infinitely many. As a consequence, if there are at least two closing contacts, then mass $N$ will interact with the obstacle infinitely many times. Accordingly, it is strictly impossible for a solution to enjoy exactly 2, 3, ..., $m$ closing contacts.

The proof of Theorem 2.2 is mainly grounded on the relative position of mass $N$ with the wall. In case 1 of the theorem, system (1.1) is linear. In case 2 with only one grazing contact, the nonlinear contact condition with zero velocity does not affect the solution of the linear system. Now we discuss the case with at least two contacts. Between every two consecutive closing contacts, say $]t_j; t_{j+1}]$, system (1.1a) with $R(t) = 0$ is linear. There is a unique solution associated with the initial data at time $t_j$. Hence, $u_N$ is a quasi-periodic trigonometric polynomial on every interval $]t_j; t_{j+1}]$. The proof of Theorem 2.2 has two steps. We first consider the quasi-periodic function $\varphi$ defined on $\mathbb{R}$ which coincides with $u_N$ on a free-flight interval. The properties of almost periodic functions dictate the behavior of $\varphi$ [5]. Then, the results obtained for $\varphi$ are extended to its restriction $\mathbf{u}(t)$. The first step is contained in the following lemma.

**Lemma 3.2** — Maximum of a quasi-periodic function. Let $\varphi$ be a quasi-periodic function defined on $\mathbb{R}$ such that

$$ \varphi(t) = \sum_{j=1}^{N} (c_j \cos \omega_j t + s_j \sin \omega_j t) \quad \text{and} \quad \varphi(0) = \sum_{j=1}^{N} c_j = d > 0 \quad (3.10) $$
If $(\omega_1, \dots, \omega_N)$ are $\mathbb{Z}$-independent, then two possibilities arise:

(1) If $s_j = 0$ and $c_j \ge 0$ for all $j = 1, \dots, N$ then $\sup_{\mathbb{R}} \varphi = d$. Moreover, if $\varphi$ is not periodic, i.e there exist at least two coefficients $c_j$ and $c_k$ with $j \ne k$ such that $c_j > 0$ and $c_k > 0$, then $\varphi(t) < d$ for all $t \ne 0$ and $\varphi(t) = d$ for only $t = 0$.

(2) Otherwise, if there exists at least one $\ell \in \{1, \dots, N\}$ such that $s_\ell \ne 0$ or $c_\ell < 0$ then $\sup_{\mathbb{R}} \varphi > d$.

The latter is equivalent to saying that the converse of the first case is true. In the first case, $\varphi(t) = \sum_{j=1}^N c_j \cos \omega_j t$ and the $\mathbb{Z}$-independence of $\{\omega_j\}$ implies that $\varphi$ is periodic if and only if there exists a unique $c_j > 0$ with $c_k = 0$ for all $k \ne j$. Then, the set $\{t : \varphi(t) = d\}$ is the infinite set $\{kT_j, k \in \mathbb{Z}\}$ where $T_j = 2\pi/\omega_j$. This corresponds to a solution with infinitely many grazing contacts. Discarding the periodic case, the value of $\varphi$ can be very close to $d$ but will never equal it again in the past and in the future. This argument proves the case when the response of system (1.1) has only one grazing contact and never reaches the obstacle again.

Later, it is proven that the function $\varphi$ with a supremum strictly greater than $d$ corresponds to a solution with infinitely many closing contacts. Now, Lemma 3.2 is proven.

*Proof.*

(1) Consider $s_j = 0$ and $c_j \ge 0$ for all $j = 1, \dots, N$, then $\varphi(t) = \sum_{j=1}^N c_j \cos \omega_j t$ and $\sup_{\mathbb{R}} \varphi = \sum_{j=1}^N |c_j| = \sum_{j=1}^N c_j = d$. If $\exists t > 0$ such that $\sum_{j=1}^N c_j \cos \omega_j t = \varphi(t) = d$, then $\sum_{j=1}^N c_j (1 - \cos \omega_j t) = 0$ where $c_j (1 - \cos \omega_j t) \ge 0$ for all $j = 1, \dots, N$. Hence, $\cos \omega_j t = 1$ for all $j = 1, \dots, N$. Thus, $\omega_j t = k_j 2\pi$, $k_j \in \mathbb{Z}$, which contradicts the $\mathbb{Z}$-independence assumption. Hence, $\varphi$ is always smaller than $d$ for $t > 0$.

(2) Otherwise, if there exists $\ell \in \{1, \dots, N\}$ such that $s_\ell \ne 0$, or $c_\ell < 0$, then [5]

$$ \sup_{\mathbb{R}} \varphi = \sum_{j=1}^{N} \sqrt{c_j^2 + s_j^2} \ge \max\left(\sum_{j=1}^{N} \sqrt{c_j^2 + s_j^2}, \sum_{j=1}^{N} |c_j|\right) > \sum_{j=1}^{N} c_j = d. $$

Hence, there exists $t > 0$ such that $\varphi(t) = d$. $\square$

Now, using Lemma 3.2, Theorem 2.2 is proven.

*Proof.* When $R(t) = 0$, the solution to (1.1a) for all $t \in \mathbb{R}$ reads

$$ \begin{aligned} \mathbf{\Phi}(t) &= \mathbf{P} \cos (t\mathbf{\Omega}) \mathbf{P}^{-1} \mathbf{u}(0) + \mathbf{P}\mathbf{\Omega}^{-1} \sin (t\mathbf{\Omega}) \mathbf{P}^{-1} \dot{\mathbf{u}}(0) \\ &= \sum_{j=1}^{N} (c_j \cos \omega_j t + s_j \sin \omega_j t) \mathbf{P} \mathbf{e}_j \end{aligned} \quad (3.11) $$

where $\mathbf{P}\mathbf{e}_j$, $j = 1, \dots, N$ are eigenvectors corresponding to eigenvalues $\omega_j^2$ of $\mathbf{M}^{-1}\mathbf{K}$; also $\cos(t\mathbf{\Omega}) \equiv \mathrm{diag}(\cos \omega_j t)_{j=1,\dots,N}$ and $\sin(t\mathbf{\Omega}) \equiv \mathrm{diag}(\sin \omega_j t)_{j=1,\dots,N}$ are used as notations. We can choose $\mathbf{P}\mathbf{e}_j$ such that $P_{Nj} = 1$, $j = 1, \dots, N$. It is clear that $\mathbf{u}(t) = \mathbf{\Phi}(t)$ as long as $u_N(t) < d$.

Henceforth, consider the solution $\mathbf{u}(t)$ to the vibro-impact system (1.1). If $u_N(t) < d$, for all $t \ge 0$, then the solution does not impact the wall. Otherwise, since the system is autonomous, the impacting time is chosen to be $t = 0$, i.e. $u_N(0) = d$. Based on the velocity of the $N^{th}$ mass at the contact time $t = 0$, two possibilities emerge

**Strictly positive pre-contact velocity:** By assumption, $\dot{u}_N^-(0) > 0$ which implies $\dot{u}_N^+(0) < 0$. Let $\varphi$ be a function defined on $\mathbb{R}$ such that $\varphi(t) = u_N(t)$ on $\{t, s \in ]0; t]\,|\, u_N(s) < d\}$. It follows that $\varphi(0) = d$ and $\dot{\varphi}(0) < 0$, thus there exists $\tau > 0$ such that $\varphi(t) > d$ for $t \in ]-\tau; 0[$. Therefore, $\sup_{\mathbb{R}^-} \varphi \ge \sup_{]-\tau; 0[} \varphi > d$. For an almost periodic function $\varphi$, the supremum taken on $\mathbb{R}^-$ is also the one taken on $\mathbb{R}^+$ [5], this yields $\sup_{\mathbb{R}^+} \varphi > d$. By Lemma 3.2, the first instant $t_1 > 0$ such that $\varphi(t_1) = d$ exists, i.e. mass $N$ will come back to the obstacle at time $t_1$.
**Vanishing pre-contact velocity:** By assumption, $\dot{u}_N^-(0) = 0$ which implies $\dot{u}_N^+(0) = 0$ as well as $\sum_{j=1}^N c_j = d$ and $\sum_{j=1}^N s_j\omega_j = 0$. Let $\varphi$ be the function defined on $\mathbb{R}$ such that $\varphi(t) = u_N(t)$ on $\{t, s \in ]0; t]\,|\, u_N(s) < d\}$. If $s_j = 0$ and $c_j \ge 0$ for all $j = 1, \dots, N$, by Lemma 3.2, it follows that $\sup_{\mathbb{R}} \varphi = d$ and $\varphi(t) < d$ for all $t > 0$. Thus, the solution has only one grazing contact: it is a linear solution.

Otherwise, Lemma 3.2 shows that $\sup_{\mathbb{R}} \varphi > d$. It follows that there is a sticking phase from 0 to $\tau$, i.e. $u_N(t) = d$ for all $t \in [0; \tau]$ and $u_N(t) < d$ for $\tau < t < \tau + \delta$, where $\delta > 0$. Assume that $\mathbf{w}$ is the solution to (1.1a) after the sticking interval with new initial data at $t = \tau$. Thus, $w_N$ reads

$$w_N(t) = \sum_{j=1}^{N} (c_j \cos(\omega_j(t-\tau)) + s_j \sin(\omega_j(t-\tau))), \quad t \ge \tau. \quad (3.12)$$

Let $\varphi$ be the function defined on $\mathbb{R}$ satisfying $\varphi(t) = w_N(t)$ on $\{t, s \in ]0; t]\,|\, w_N(s) < d\}$. Consider the solution $\mathbf{u}$ to the vibro-impact system (1.1) just before and just after the sticking phase. Its component $u_N(t)$ can be written as

$$u_N(t) = \begin{cases} \varphi(t) & \text{for } t \le 0 \\ d & \text{for } t \in [0; \tau] \\ \varphi(t) & \text{for } t \ge \tau \end{cases} \quad (3.13)$$

We then show that $\sup_{[\tau; \infty[} \varphi > d$. Otherwise, $\sup_{\mathbb{R}} \varphi = d$, then $\varphi(t) < d$, for all $t > \tau$ and $\varphi(\tau) = d$. By reversibility ($e = 1$) and uniqueness arguments, the solution to (1.1) never activates the contact condition and there is only one grazing contact at $\tau$. This is in contradiction with the fact that the solution involves a sticking phase. As a consequence, $\sup_{\mathbb{R}} \varphi > d$. Hence, there exists $t > \tau$ such that $\varphi(t) = d$.

The same procedure can be repeated infinitely many times in both cases, thus the set of closing contacts is countably infinite. $\square$

**3.3. Poincaré section.** This section of the paper deals with the construction of the Poincaré section on which the First Return Map is well-defined. Consider an orbit $[\mathbf{u}, \dot{\mathbf{u}}] \subset \mathbb{R}^{2N}$ satisfying (1.1). The First Return Map (also called Poincaré map) is used to study the dynamics in a neighborhood of such an orbit. This map is defined on a Poincaré section which is classically a $(2N-1)$-dimensional manifold in $\mathbb{R}^{2N}$ that contains a point $\mathbf{U}(t) = [\mathbf{u}(t), \dot{\mathbf{u}}(t)]$ of the previous orbit and is transverse to the orbit at $\mathbf{U}(t)$. In the current work, the transversality is lost on $\mathcal{H}^0$ since an orbit starting on $\mathcal{H}^0$ may not intersect $\mathcal{H}$ again. Hence, an important task is to eliminate the set of initial data such that the associated orbit does not intersect $\mathcal{H}_\varphi$ again. This is achieved by investigating the structure of $\mathcal{H}^0$. Via Theorem 2.2, there are two possibilities after the grazing contact: the orbit never comes in contact again ($\mathcal{H}_1^0$) or the orbit comes in contact infinitely many times ($\mathcal{H}_\infty^0$). The explicit description of $\mathcal{H}_1^0$ is given in the proposition below.

**Proposition 3.3** — $\mathcal{H}_1^0$ and solutions with only one grazing contact. The set $\mathcal{H}_1^0$ of initial data on $\mathcal{H}$ such that the associated orbits have only one contact is the convex subset of $\mathcal{H}^0$,

$$\{[\mathbf{u}^\top, \dot{\mathbf{u}}^\top]^\top \in \mathbb{R}^{2N}, u_N = d, \dot{u}_k = 0, k = 1, \dots, N\} \quad (3.14)$$

of a $(N-1)$-dimensional affine subspace such that the $N-1$ components $u_1, u_2, \dots, u_{N-1}$ satisfy the inequalities:

$$P_{Nk} \sum_{j=1}^{N-1} B_{kj} u_j \ge -P_{Nk} B_{kN} d, \quad \forall k = 1, \dots, N \quad (3.15)$$

with two of which, at least, are strict inequalities and where $\mathbf{B} = (B_{ij})_{i,j=1}^N$ denotes $\mathbf{P}^{-1}$.
Note that $\mathcal{H}_1^0$ is a $N-1$-dimensional convex set or an empty set.

*Proof.* The solution to (1.1) with a unique grazing contact is analytic and $u_N(t)$ has a closed form-expression:

$$ 
\begin{aligned}
u_N(t) &= \mathbf{e}_N^\top (\mathbf{P} \cos(t\Omega) \mathbf{P}^{-1} \mathbf{u}(0) + \mathbf{P}\Omega^{-1} \sin(t\Omega) \mathbf{P}^{-1} \dot{\mathbf{u}}(0)) \\
&= \sum_{k=1}^{N} (c_k \cos \omega_k t + s_k \sin \omega_k t)
\end{aligned} 
\quad (3.16) $$

where $c_k = P_{Nk} \mathbf{e}_k^\top \mathbf{P}^{-1} \mathbf{u}(0)$ and $s_k = P_{Nk} \omega_k^{-1} \mathbf{e}_k^\top \mathbf{P}^{-1} \dot{\mathbf{u}}(0)$, $k = 1, \dots, N$. Since $[\mathbf{u}^\top, \dot{\mathbf{u}}^\top]^\top$ belongs to $\mathcal{H}_1^0$, this corresponds to case 1 of Lemma 3.2, thus $u_N(0) = d$ and $\dot{u}_N(0) = 0$, and the coefficients $(c_k, s_k)$ satisfy

$$ 
\left\{ 
\begin{array}{ll}
c_k \ge 0, & \forall k = 1, \dots, N, \quad \exists \ell \neq m : c_\ell > 0, c_m > 0 \\
s_k = 0, & \forall k = 1, \dots, N.
\end{array}
\right. 
\quad (3.17) $$

This gives

$$ 
\left\{ 
\begin{array}{ll}
P_{Nk} \mathbf{e}_k^\top \mathbf{P}^{-1} \mathbf{u}(0) \ge 0, & \forall k = 1, \dots, N, \quad \exists \ell \neq m : c_\ell > 0, c_m > 0 \\
P_{Nk} \omega_k^{-1} \mathbf{e}_k^\top \mathbf{P}^{-1} \dot{\mathbf{u}}(0) = 0, & \forall k = 1, \dots, N.
\end{array}
\right. 
\quad (3.18) $$

The second condition yields a linear system $\mathbf{P}^{-1}\dot{\mathbf{u}}(0) = \mathbf{0}$. Since $\det(\mathbf{P}^{-1}) \neq 0$, it follows that $\dot{\mathbf{u}}(0) = \mathbf{0}$ and $u_1, \dots, u_{N-1}$ satisfy the inequalities (3.18). This gives the explicit formula of $\mathcal{H}_1^0$ stated in Proposition 3.3. $\square$

An immediate consequence of Proposition 3.3 is that most orbits which belong to $\mathcal{H}^0$ are in $\mathcal{H}_\infty^0$. To be precise, the set of initial data such that the associated orbits belong to $\mathcal{H}_1^0$ is empty or is a $(N-1)$-dimensional subset of the $2(N-1)$-dimensional set $\mathcal{H}^0$. As a consequence, the Poincaré section chosen in Definition 2.1 is reasonable.

**Corollary 3.1** — Domain of definition of the First Return Map. The maximal subset of $\mathcal{H}$ where the First Return Map $\mathcal{F}$ is well-defined is $\mathcal{H}_\mathcal{P} = \mathcal{H}^- \cup \mathcal{H}_\infty^0$.

Consequences of Proposition 3.3 are now stated.

**Corollary 3.2.** Orbits with a sticking phase exhibit infinitely many countable closing contacts.

*Proof.* An orbit including a sticking phase intersects $\mathcal{H}^0$. As proven in Proposition 3.3, the initial data must be in $\mathcal{H}_\infty^0$ since $\mathcal{H}_1^0$ only involves the data on $\mathcal{H}$ such that the orbits have a unique grazing contact but no sticking phase. Hence, it belongs to the set with an infinite number of closing contacts. $\square$

**Corollary 3.3** — Infinite number of closing contacts. If the orbit intersects $\mathcal{H}^-$ then it intersects $\mathcal{H}^-$ infinitely many times in the future and in the past.

*Proof.* By Theorem 2.2, an orbit corresponding to the initial data in $\mathcal{H}^-$ experiences at least one impact at $t = 0$. This eliminates the possibilities of case 1 and case 2 in the theorem. Hence, the orbit belongs to the third category, which means that there will be an infinite number of closing contacts. However, if the set of closing contacts involves only grazing contacts in the future, the system becomes linear with $u_N(t) \le d$ for all $t > 0$. Let $\varphi$ be a function defined on $\mathbb{R}$ with $\varphi(t) = u_N(t)$, it follows that $\sup_{\mathbb{R}} \varphi = \sup_{t>0} \varphi = d$. This contradicts the fact that the supremum of $\varphi$ must be greater than $d$ since there is at least one impact at $t = 0$. The process can be repeated to get another closing contact and so on. $\square$

The subset $\mathcal{H}^--$ of $\mathcal{H}^-$ containing all the initial data generating orbits with impacts only seems to be dense in $\mathcal{H}^-$. Indeed, this is an interesting open problem.
About the set $\mathcal{H}_1^0$ for two-dof systems.

**Corollary 3.4** — $\mathcal{H}_1^0$ for a two-dof system. Consider $N = 2$ in (1.1), then, from Proposition 3.3, $\mathcal{H}_1^0$ is defined as

$$ \begin{aligned} \mathcal{H}_1^0 &= \{[\mathbf{u}^\top, \dot{\mathbf{u}}^\top]^\top \in \mathbb{R}^4 : P_{2k}B_{k1}u_1 > -P_{2k}B_{k2}d, k = 1, 2, u_2 = d, \dot{u}_1 = \dot{u}_2 = 0\} \\ &= D_{u_1} \times \{d\} \times \{0\} \times \{0\} \end{aligned} \quad (3.19) $$

where $D_{u_1}$ is the subset of $\mathbb{R}$ in which $u_1$ satisfies the two strict inequalities $P_{21}B_{11}u_1 > -P_{21}B_{12}d$ and $P_{22}B_{21}u_1 > -P_{22}B_{22}d$.

Therefore, the set $D_{u_1}$ can be either empty or an open interval of the form $(b, \infty)$, $(-\infty, a)$, or $(a, b)$, with $a = \min\{\alpha_1, \alpha_2\}$, $b = \max\{\alpha_1, \alpha_2\}$, where $\alpha_1 = -dB_{22}/B_{21}$ and $\alpha_2 = -dB_{12}/B_{11}$.

**Remark 3.1** — Linear gazing modes of a two-dof system and boundary of $\mathcal{H}_1^0$. The two scalars $\alpha_1$ and $\alpha_2$ are the distinct initial values of $u_1$ corresponding to the first and the second grazing linear mode.

The set $\mathcal{H}_1^0$ is an interval in the two-dimensional space $\mathcal{H}^0$. If $\mathcal{H}_1^0 \neq \emptyset$ then

$$ \mathcal{H}_1^0 = ]a, b[ \times \{d\} \times \{0\} \times \{0\}, \quad (3.20) $$

in the system of coordinates $(u_1, u_2, v_1, v_2)$. This means that the set of initial data generating orbits that have only one grazing contact is a very small subset of all the initial data of orbits which contain zero velocity closing contacts.

For the 2-dof system, the Poincaré section $\mathcal{H}_\rho$ reads

$$ \mathcal{H}_\rho = \mathcal{H}^- \cup \mathcal{H}_\infty^0 = \{(u_1, u_2 = d, v_1, v_2), u_1 \in \mathbb{R}, v_1 \in \mathbb{R}, v_2 > 0\} \cup \mathcal{H}_\infty^0. \quad (3.21) $$

It is displayed in figure 3. The only complicated part of $\mathcal{H}_\rho$ is $\mathcal{H}_\infty^0 = \mathcal{H}^0 - \mathcal{H}_1^0$. It means that the interval $\mathcal{H}_1^0$ must be subtracted from the 2-dimensional set $\mathcal{H}^0$. The set $\mathcal{H}^0$ is studied and depicted in figure 6 at the end of Section 4.2. In figure 6, the data generating a grazing or a sticking phase is also completely explicit.

<figure><img src="image_5.png" /><figcaption>FIGURE 3. The 3D Poincaré section for the 2 dof case on $\{u_2 = d\}$. It consists in the upper part $\{\dot{u}_2 > 0\}$ and the part $\{\dot{u}_2 = 0\}$ minus the interval $]a, b[$.</figcaption></figure>

## 4. IMPLICIT FUNCTION THEOREM AND POWER-ROOT SINGULARITY

This section is divided into three parts. First, the square-root singularity of maps is investigated in a general framework in order to apply it on the First Return Time. This singularity generically induces the square-root singularity for the First Return Map near the grazing linear modes as stated later in Section 6. Second, a general power-root singularity of maps is then unveiled. Finally, the discontinuity of the First Return Time is shortly discussed in the third part.
**4.1. Square-root singularity.** The First Return Time $T$ is implicitly defined in the equation $f(T, W) = 0$ reflecting the equality $u_N(T) = d$. Unfortunately, if the contact at $T_0$ is of grazing type, $\partial_t f(T_0, W_0)$ vanishes and hence the implicit function theorem does not apply in a straightforward fashion. Indeed, a square-root singularity is expected at the intersection of the hypersurfaces $f = 0$ and $\partial_t f = 0$. The implicit function theorem is then exploited on a variable other than $t$. We state the results for the function $f$ defined in a two-dimensional space as well as a $m+2$-dimensional space with $m \ge 1$.

The presentation of the square-root singularity is elementary and self-contained. It corresponds to the local geometry of a folded map, well known in the singularity and bifurcation theory, see for instance the pioneering work [14].

**4.1.1. In two dimensions.** Let $f$ be a function of two variables $x$ and $y$, where the relation of $x$ with respect to $y$ is implicitly given by $f(x, y) = 0$. Let us write $x$ locally as a function of $y$ when $f$ satisfies some unusual conditions as follows.

**Theorem 4.1—Square-root singularity in two dimensions.** Suppose that $f(x, y) \in C^3(\mathbb{R}^2, \mathbb{R})$ satisfies the following conditions at $(x_0, y_0)$:

(1) $f(x_0, y_0) = 0$,

(2) $\partial_x f(x_0, y_0) = 0$,

(3) $\partial_x^2 f(x_0, y_0) \neq 0$,

(4) $\partial_y f(x_0, y_0) \neq 0$.

Denote the ratio $2\gamma = -\partial_x^2 f(x_0, y_0)/\partial_y f(x_0, y_0) \neq 0$ and $s_\gamma = \text{sign} \gamma$. There exist two neighborhoods $I$ and $J$ of $y_0$ and $x_0$, respectively, as well as a function $\psi \in C^3(I_\gamma, J)$ defined on the subinterval

$$I_\gamma = \{y \in I, s_\gamma(y - y_0) \ge 0\} \qquad (4.1)$$

such that the set $\{(x, y) \in I \times J, f(x, y) = 0\}$ is formed by the following two branches near $(x_0, y_0)$:

$$x = x_0 + \psi(\sqrt{s_\gamma(y - y_0)}) \quad \text{or} \quad x = x_0 + \psi(-\sqrt{s_\gamma(y - y_0)}). \qquad (4.2)$$

Moreover, $\psi$ satisfies $\psi(0) = 0$ and $\dot{\psi}(0) = 1/\sqrt{|\gamma|} \neq 0$.

An illustration of this theorem is given in Figure 4. Note that $\psi$ is not defined on the whole open interval $I$ but only on half of it, determined by the sign of $\gamma$.

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Two branches on the right of the line $y = y_0$ when $\gamma > 0$.</figcaption></figure>
*Proof.* Without loss of generality, assume that $(x_0, y_0)$ is at the origin. If not, consider the new variables $x^* = x - x_0$ and $y^* = y - y_0$ so that at the point $(0,0)$, the function $f^*(x^*, y^*) = f(x_0 + x^*, y_0 + y^*)$ satisfies the conditions of the theorem. The conclusions are then modified accordingly.

The implicit function theorem guarantees the existence of the neighborhoods $V_x$ of $x_0 = 0$ and $V_y$ of $y_0 = 0$, and a unique function $y = \varphi(x)$ for $x \in V_x$ and $y \in V_y$ such that $f(x, \varphi(x)) = 0$ for $x \in V_x$, $\varphi(0) = 0$, $\dot{\varphi}(0) = -\partial_x f(0, 0)/\partial_y f(0, 0) = 0$ and $\ddot{\varphi}(0) = -\partial_x^2 f(0, 0)/\partial_y f(0, 0) = 2\gamma \neq 0$. Moreover, $\varphi$ is as smooth as $f$. The Taylor expansion of $\varphi$ with an integral remainder at 0 is

$$y = \varphi(x) = x^2 r(x) \quad \text{with} \quad r(x) = \int_0^1 (1-s) \ddot{\varphi}(sx) \, ds \qquad (4.3)$$

where $r \in C^1(\mathbb{R}, \mathbb{R})$ and $r(0) = \ddot{\varphi}(0)/2 = \gamma \neq 0$. The expression of $x$ with respect to $y$ depends on the sign of $\ddot{\varphi}(0)$ (or the sign of $\gamma$) as follows:

**If $\gamma > 0$:** then $r(0) > 0$ and the continuity of $r$ implies that $r(x) > 0$ in some neighborhood of 0. On that neighborhood, let us define the continuously differentiable function $\phi(x) = x\sqrt{r(x)}$: at $x = 0$, $\phi$ has a nonzero derivative since $r(0) \neq 0$. Through the inverse function theorem, there exists an inverse function $\phi^{-1}$ in some neighborhood $I$ of $\phi(0) = 0$. From $y = \phi^2(x)$, it is required that $y \ge 0$. Accordingly, $I_\gamma$ is simply a right neighborhood of 0. Hence, for $y \in I_\gamma$, $x = \phi^{-1}(\sqrt{y})$ or $x = \phi^{-1}(-\sqrt{y})$.

**If $\gamma < 0$:** then $y = -\phi^2(x)$ with $\phi(x) = x\sqrt{-r(x)}$. A similar proof holds for $y$ defined on $I_\gamma$ which is now a left neighborhood of 0. The conclusion is obtained by replacing $\sqrt{y}$ by $\sqrt{-y}$ and the function $\psi$ satisfies $\dot{\psi}(0) = \dot{\phi}^{-1}(0) = 1/\sqrt{-\gamma}$.

In both cases, by denoting $\psi = \phi^{-1}$, we have shown that there exists a function $\psi$ defined in a neighborhood $I_\gamma$ of 0 satisfying $\psi(0) = 0$ and $\dot{\psi}(0) = \dot{\phi}^{-1}(0) = 1/\sqrt{|\gamma|}$ such that $x = \psi(\sqrt{s_\gamma y})$ or $x = \psi(-\sqrt{s_\gamma y})$. This concludes the proof. $\square$

**Remark 4.1.** Another way to express the square-root singularity of $x$ is $x = \sqrt{\beta(y)}$. However, with the present expression, say $x = \lambda(\sqrt{y})$, we have a better regularity of $\lambda$. For instance, if $\lambda(y) = y+y^2$, $y \in \mathbb{R}^+$, then $\lambda \in C^\infty(\mathbb{R}^+)$. However, $\beta(y) = \lambda^2(\sqrt{y}) = (\sqrt{y} + y)^2 = y + 2y\sqrt{y} + y^2$ is not smooth: $\beta \notin C^\infty$ but $\beta \in C^{1.5}$.

4.1.2. In $m+2$ dimensions ($m \ge 1$). Consider a general function $f$ of $m+2$ variables with a non-zero gradient.

**Theorem 4.2** — Square-root singularity in $m+2$ dimensions. Consider the smooth function $f(x, y, \mathbf{z}) : \mathbb{R} \times \mathbb{R} \times \mathbb{R}^m \to \mathbb{R}$ such that at point $\mathbf{X}_0 = (x_0, y_0, \mathbf{z}_0)$, it satisfies the following conditions:

(1) $f(\mathbf{X}_0) = 0,$

(2) $\partial_x f(\mathbf{X}_0) = 0,$

(3) $\partial_x^2 f(\mathbf{X}_0) \neq 0,$

(4) $\partial_y f(\mathbf{X}_0) \neq 0.$

Denote $2\gamma = -\partial_x^2 f(\mathbf{X}_0)/\partial_y f(\mathbf{X}_0) \neq 0$ and $s_\gamma = \text{sign} \, \gamma$. Then, there exist three neighborhoods $V_{x_0}, V_{y_0}$, and $V_{z_0}$ of $x_0, y_0$ and $z_0$ respectively, and two smooth scalar functions $\eta : V_{z_0} \to V_{x_0}$ and $\alpha : V_{z_0} \to V_{y_0}$ satisfying $\eta(\mathbf{z}_0) = x_0$ and $\alpha(\mathbf{z}_0) = y_0$ such that the set

$$S_c = \{(x, y, \mathbf{z}) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}^m, f(x, y, \mathbf{z}) = 0 \text{ and } \partial_x f(x, y, \mathbf{z}) = 0\} \qquad (4.4)$$

is parameterized by $S_c \cap \Omega = \{(\eta(\mathbf{z}), \alpha(\mathbf{z}), \mathbf{z}), \mathbf{z} \in V_{z_0}\}$ where $\Omega = V_{x_0} \times V_{y_0} \times V_{z_0}$. Let $B_\gamma$ be a subset of $V_{y_0} \times V_{z_0}$ where $B_\gamma = \{(y, \mathbf{z}) \in V_{y_0} \times V_{z_0}, s_\gamma(y - \alpha(\mathbf{z})) \ge 0\}$. There exists a smooth scalar function $\psi$ defined on $B_\gamma$ such that there are two graphs:

$$x = \eta(\mathbf{z}) + \psi(\sqrt{s_{\gamma}(y - \alpha(\mathbf{z}))}, \mathbf{z}), \quad (y, \mathbf{z}) \in B_{\gamma}, \qquad (4.5)$$
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Square-root singularity in dimension 3. Two graphs of $x(y, \mathbf{z})$ are shown in red and blue. The square-root singularity arises along the curve $y = \alpha(\mathbf{z})$ in purple.</figcaption></figure>

and

$$x = \eta(\mathbf{z}) + \psi(-\sqrt{s_{\gamma}(y - \alpha(\mathbf{z}))}), \quad (y, \mathbf{z}) \in B_{\gamma}. \qquad (4.6)$$

In particular, $\psi(0, \mathbf{z}_0) = 0$ and $\partial_y\psi(0, \mathbf{z}_0) = 1/|\gamma|$.

An illustration of the square-root singularity in three dimensions is provided in Figure 5. The function $\psi$ is not defined on the whole neighborhood $V_{y_0} \times V_{z_0}$ but only on $B_{\gamma}$. The subset $B_{\gamma}$ in the theorem gives the one-sided condition from which the square-root singularity arises. When the sign of $\gamma$ is positive or negative, respectively, $B_{\gamma}$ is simply a region above or below the hypersurface $y = \alpha(\mathbf{z})$ and includes it. A more general result when condition (3) in Theorem 4.2 does not hold is discussed in Remark 4.2.

*Proof.* Let $S_c$ be the manifold of $m$ dimensions which is the intersection of the two hypersurfaces of dimension $m+1$ near $(x_0, y_0, \mathbf{z}_0)$, $S_c = S^0 \cap S^1$ where $S^0 = \{(x, y, \mathbf{z}) \in \Omega, f(x, y, \mathbf{z}) = 0\}$ and $S^1 = \{(x, y, \mathbf{z}) \in \Omega, G(x, y, \mathbf{z}) = 0\}$. Moreover, it can be parameterized in a neighborhood of $\mathbf{z}$ by two parametric functions $\eta$ and $\alpha$ as follows.

From (1) and (4) (in Theorem 4.2), by the implicit function theorem, there exist two neighborhoods $V$ of $(x_0, \mathbf{z}_0)$ and $W$ of $y_0$, and a function $\varphi: V \to W$, $(x, \mathbf{z}) \mapsto y = \varphi(x, \mathbf{z})$ such that $f(x, \varphi(x, \mathbf{z}), \mathbf{z}) = 0$ for $(x, \mathbf{z}) \in V$, $\varphi(x_0, \mathbf{z}_0) = y_0$. In particular, $\partial_x\varphi(x_0, \mathbf{z}_0) = 0$ and $\partial_x^2\varphi(x_0, \mathbf{z}_0) = -\partial_x^2f(x_0, y_0, \mathbf{z}_0)/\partial_yf(x_0, y_0, \mathbf{z}_0) = 2\gamma \neq 0$. Moreover, $\varphi$ inherits the smoothness of $f$.

For $(x, \mathbf{z}) \in V$, once again we apply the implicit function theorem for the function $G(x, \mathbf{z}) = \partial_x f(x, \varphi(x, \mathbf{z}), \mathbf{z})$ which satisfies:

$$G(x_0, \mathbf{z}_0) = \partial_x f(x_0, y_0, \mathbf{z}_0) = 0$$

$$\partial_x G(x_0, \mathbf{z}_0) = \partial_x^2 f(x_0, y_0, \mathbf{z}_0) + \partial_y \partial_x f(x_0, y_0, \mathbf{z}_0) \partial_x \varphi(x_0, \mathbf{z}_0) = \partial_x^2 f(x_0, a_0, \mathbf{z}_0) \neq 0.$$

This guarantees the existence of the neighborhoods $B_{z_0}$ of $z_0$ and $B_{x_0}$ of $x_0$ such that there is a smooth function $\eta: B_{z_0} \to B_{x_0}$, $\mathbf{z} \mapsto x = \eta(\mathbf{z})$ satisfying $x_0 = \eta(\mathbf{z}_0)$, $G(\eta(\mathbf{z}), \mathbf{z}) = 0$ for all $\mathbf{z} \in B_{z_0}$. It follows that $\partial_x f(\eta(\mathbf{z}), \alpha(\mathbf{z}), \mathbf{z}) = 0$ for all $\mathbf{z} \in B_{z_0}$ where $\alpha(\mathbf{z}) := \varphi(\eta(\mathbf{z}), \mathbf{z})$ belongs to a neighborhood of $y_0$, denoted $V_{y_0}$. In particular, $\alpha(\mathbf{z}_0) = \varphi(x_0, \mathbf{z}_0) = y_0$. As a consequence,
the local parameterization

$$S_c = \{(x, y, \mathbf{z}) = (\eta(\mathbf{z}), \alpha(\mathbf{z}), \mathbf{z}), \mathbf{z} \in B_{z_0}\} \qquad (4.7)$$

exists. This parameterization is also found by using the implicit function theorem for the vector function $H = (f, \partial_x f)$ with the invertible matrix

$$D_{x,y}H(\mathbf{X}_0) = \begin{bmatrix} \partial_x f & \partial_y f \\ \partial_x^2 f & \partial_y \partial_x f \end{bmatrix} (\mathbf{X}_0) = \begin{bmatrix} 0 & \neq 0 \\ \neq 0 & \partial_y \partial_x f(\mathbf{X}_0) \end{bmatrix}. \qquad (4.8)$$

We will now show that the square-root singularity arises along the hypersurface $y = \alpha(\mathbf{z})$. For each fixed $\mathbf{z} \in B_{z_0}$, the Taylor expansion with an integral remainder of $\varphi$ with respect to $x$ near $\eta(\mathbf{z})$ reads

$$\varphi(x, \mathbf{z}) = \varphi(\eta(\mathbf{z}), \mathbf{z}) + \partial_x \varphi(\eta(\mathbf{z}), \mathbf{z})(x - \eta(\mathbf{z})) \\ + (x - \eta(\mathbf{z}))^2 \int_0^1 (1-s) \partial_x^2 \varphi(s(x - \eta(\mathbf{z}))) + \eta(\mathbf{z}), \mathbf{z} \, \mathrm{d}s. \qquad (4.9)$$

Since $y = \varphi(x, \mathbf{z})$, $\partial_x \varphi(\eta(\mathbf{z}), \mathbf{z}) = 0$, $\varphi(\eta(\mathbf{z}), \mathbf{z}) = \alpha(\mathbf{z})$, this yields

$$y - \alpha(\mathbf{z}) = (x - \eta(\mathbf{z}))^2 r(x - \eta(\mathbf{z}), \mathbf{z}) \qquad (4.10)$$

where $r(x - \eta(\mathbf{z}), \mathbf{z})$ is a smooth function with respect to $x$ and $\mathbf{z}$, satisfying $r(0, \mathbf{z}) = \partial_x^2 \varphi(\eta(\mathbf{z}), \mathbf{z})/2$.

The remainder $r(x - \eta(\mathbf{z}), \mathbf{z})$ has the sign of $\gamma$ for all $\mathbf{z}$ in a neighborhood of $z_0$. Consider the function $\gamma(\mathbf{z}) = r(0, \mathbf{z})$, $\gamma(\mathbf{z})$: it is as smooth as $r$ and satisfies $\gamma(z_0) = \partial_x^2 \varphi(x_0, \mathbf{z})/2 = \gamma \neq 0$. Without loss of generality, assume that $\gamma > 0$. By continuity, $\gamma(\mathbf{z}) > 0$ in some neighborhood of $z_0$, say $V_{z_0}$. It follows that $r(0, \mathbf{z}) > 0$ for $\mathbf{z} \in V_{z_0}$. Similarly, by the continuity of $r$ with respect to $x$, $r(x, \mathbf{z}) > 0$ for all $x$ in a neighborhood of $\eta(\mathbf{z})$.

Let $B_\gamma \subset V_{y_0} \times V_{z_0}$ be the region adjacent to and including the hypersurface $y = \alpha(\mathbf{z})$ such that $B_\gamma = \{(y, \mathbf{z}) \in V_{y_0} \times V_{z_0}, \gamma(y - \alpha(\mathbf{z})) \geq 0\}$ and denote by $\phi$ the function $\phi(x - \eta(\mathbf{z}), \mathbf{z}) = (x - \eta(\mathbf{z}))\sqrt{r(x - \eta(\mathbf{z}), \mathbf{z})}$. Equation (4.10) then becomes

$$y - \alpha(\mathbf{z}) = (\phi(x - \eta(\mathbf{z}), \mathbf{z}))^2. \qquad (4.11)$$

This follows that $y - \alpha(\mathbf{z}) \geq 0$. The subset $B_\gamma$ is then a region above and including the hypersurface $y = \alpha(\mathbf{z})$. Since $\phi$ satisfies $\phi(0, \mathbf{z}) = 0$, $\partial_x \phi(0, \mathbf{z}) = \sqrt{r(0, \mathbf{z})} = \sqrt{\gamma(\mathbf{z})} \neq 0$, there exists an inverse function $\phi^{-1}$ in some neighborhood of 0. The inverse function theorem is used uniformly with respect to the parameter $\mathbf{z}$ and eventually reduces the neighborhood of $z_0$ and $x_0$. Therefore, for $(y, \mathbf{z}) \in B_\gamma$, Equation (4.11) leads to $\phi(x - \eta(\mathbf{z}), \mathbf{z}) = \sqrt{y - \alpha(\mathbf{z})}$ or $\phi(x - \eta(\mathbf{z}), \mathbf{z}) = -\sqrt{y - \alpha(\mathbf{z})}$. This implies $x = \eta(\mathbf{z}) + \phi^{-1}(\sqrt{y - \alpha(\mathbf{z})}, \mathbf{z})$ or $x = \eta(\mathbf{z}) + \phi^{-1}(-\sqrt{y - \alpha(\mathbf{z})}, \mathbf{z})$. In particular, $\phi^{-1}(0, \mathbf{z}) = 0$ and $\partial_y \phi^{-1}(0, \mathbf{z}) = 1/\sqrt{\gamma(\mathbf{z})}$.

Similarly, for the case $\gamma < 0$, $r(x - \eta(\mathbf{z}), \mathbf{z}) < 0$ for $x$ in some neighborhood of $\eta(\mathbf{z})$, it follows that $y - \alpha(\mathbf{z}) = -\phi^2(x - \eta(\mathbf{z}), \mathbf{z})$, where $\phi(x - \eta(\mathbf{z}), \mathbf{z}) = (x - \eta(\mathbf{z}))\sqrt{-r(x - \eta(\mathbf{z}), \mathbf{z})}$. The subset $B_\gamma$ is now the region below and including the hypersurface $y = \alpha(\mathbf{z})$. The same conclusions hold with a change of $\sqrt{y - \alpha(\mathbf{z})}$ in $\sqrt{\alpha(\mathbf{z})} - y$ and $\phi^{-1}(0, \mathbf{z}) = 1/\sqrt{-\gamma(\mathbf{z})}$.

Consequently, by denoting $\psi = \phi^{-1}$, we have shown that $\psi$ is defined in $B_\gamma$ and

$$x = \eta(\mathbf{z}) + \psi(\sqrt{s_\gamma(y - \alpha(\mathbf{z}))}, \mathbf{z}) \quad \text{or} \quad x = \eta(\mathbf{z}) + \psi(-\sqrt{s_\gamma(y - \alpha(\mathbf{z}))}, \mathbf{z}). \qquad (4.12)$$

In particular, $\psi(0, z_0) = 0$, $\partial_y \psi(0, z_0) = 1/\sqrt{\gamma(z_0)} = 1/|\gamma|$. This ends the proof. $\square$

4.2. **Power-root singularity.** The square-root singularity shown in Theorem 4.2 stems from the non-vanishing second derivative $\partial_x^2 f(x_0, y_0, z_0)$. If it vanishes, then higher order derivatives of $f$ must be considered. The following remark gives analogous conditions from which we define power-root singularity.
**Remark 4.2.** Let $n \le m + 1$ and $f(x, y, \mathbf{z}) : \mathbb{R} \times \mathbb{R} \times \mathbb{R}^m \to \mathbb{R}$ be a smooth function, such that at point $\mathbf{X}_0 = (x_0, y_0, \mathbf{z}_0)$, it satisfies the following conditions:

(1) $f(\mathbf{X}_0) = \partial_x f(\mathbf{X}_0) = \partial_x^2 f(\mathbf{X}_0) = \dots = \partial_x^{n-1} f(\mathbf{X}_0) = 0$, $\partial_x^n f(\mathbf{X}_0) \ne 0$, $n \ge 3$,

(2) $\partial_y f(\mathbf{X}_0) \ne 0$,

(3) linear independence of vectors $\nabla_{(x,y,z)} \partial_x^\ell f(x_0, y_0, \mathbf{z}_0)$, $0 \le \ell < n$, then a $n^{th}$-root singularity emerges.

Let us explain this remark in details. Denote $\gamma = -f^{(n)}(\mathbf{X}_0)/n! \partial_y f(\mathbf{X}_0) \ne 0$ and $s_\gamma = \text{sign } \gamma$. There exist three neighborhoods $V_{x_0}, V_{y_0}$, and $V_{z_0}$ of $x_0, y_0$ and $z_0$ respectively, and two smooth scalar functions $\eta : V_{z_0} \to V_{x_0}, \alpha : V_{z_0} \to V_{y_0}$ such that the critical set

$$S_c = \{(x, y, \mathbf{z}) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}^m, \partial_x^\ell f(x, y, \mathbf{z}) = 0, \forall \ell = 0, \dots, n-1\} \quad (4.13)$$

can be parametrized under the classical assumption (3) of linear independence. The set $S_c$ has dimension $m + 2 - n$. Therefore the functions $\eta$ and $\alpha$ do not depend on all components of $\mathbf{z}$. Two cases arise when solving the equation $f(x, y, \mathbf{z}) = 0$ on a small box $\Omega = V_{x_0} \times V_{y_0} \times V_{z_0}$:

$$(1) \quad n = 2\ell, \ell \in \mathbb{N}: \text{ let } B_\gamma \text{ be subset of } V_{y_0} \times V_{z_0} \text{ defined by} \\ B_\gamma = \{(y, \mathbf{z}) \in V_{y_0} \times V_{z_0}, \gamma(y - \alpha(\mathbf{z})) \ge 0\}. \quad (4.14)$$

Then there exists a smooth real function $\psi$ such that $\forall \mathbf{z} \in B_\gamma$, there are only two branches in $\Omega$:

$$x = \eta(\mathbf{z}) + \psi(\pm(s_\gamma(y - \alpha(\mathbf{z})))^{\frac{1}{n}}, \mathbf{z}), \quad (4.15)$$

$$(2) \quad n = 2\ell + 1, \ell \in \mathbb{N}: \text{ there exists an interval } I \text{ containing } \alpha(\mathbf{z}), \text{ and a smooth function } \psi \text{ defined on } I \times V_{z_0} \text{ such that}$$

$$x = \eta(\mathbf{z}) + \psi((y - \alpha(\mathbf{z}))^{\frac{1}{n}}, \mathbf{z}). \quad (4.16)$$

In both cases, $\psi$ is not degenerate at $(0, \mathbf{z}_0)$ with respect to the first variable: $\psi(0, \mathbf{z}_0) = 0$ and $\partial_y \psi(0, \mathbf{z}_0) = 1/\gamma^{\frac{1}{n}} \ne 0$.

The existence of a cube-root or fourth-root singularity is investigated in [4]. In the present work, the so-called power-root singularity is defined for any positive integer $n$.

**Definition 4.1** — **Power-root singularity.** A function $F$ defined on a subset of $\mathbb{R}^m$ is said to have an $n^{th}$-root singularity at 0 if there exists $n > 0$ such that, by a change of variable if needed, $F$ can be written as

$$F(\mathbf{X}) = f(X_1^{\frac{1}{n}}, X_2, \dots, X_m) \quad (4.17)$$

where $f$ is a smooth function; $n$ is called the singularity exponent.

Using this definition, Theorem 4.2 states the conditions needed for a power-root singularity. This singularity is associated to the multiplicity of a root of a function whose definition is recalled below.

**Definition 4.2** — **Multiplicity.** Given a smooth real function $f$, a positive integer $n$ is said to be the multiplicity of a root of $f$, denoted by $\text{mult}(f)(r) = n$, if it satisfies

$$f(r) = \dot{f}(r) = \dots = f^{(n-1)}(r) = 0 \quad \text{and} \quad f^{(n)}(r) \ne 0 \quad (4.18)$$

By convention, we denote $\text{mult}(f)(r) = 0$ when $f(r) \ne 0$. Before stating the main results of this section, recall that $\underline{\mathbf{M}} = \text{diag}(m_j)_{j=1}^{N-1}$ and $\underline{\mathbf{K}} = (k_{ij})_{i,j=1}^{N-1}$ are the matrices obtained from $\mathbf{M}$ and $\mathbf{K}$ by removing their last row and column, with $\underline{\mathbf{I}}_N = \mathbf{e}_N^T \underline{\mathbf{K}}$ as the last row of $\mathbf{K}$ and $\underline{\mathbf{I}}_N$ is the last row of $\mathbf{K}$ removing the last entry $k_{NN}$. Lemma 4.3 requires the following assumptions on the matrix $\underline{\mathbf{D}} = \underline{\mathbf{M}}^{-1} \underline{\mathbf{K}}$ and the matrix $\underline{\mathbf{P}}$ of eigenvectors of $\underline{\mathbf{M}}^{-1} \underline{\mathbf{K}}$ (see Section 1).

**Assumption 4.1.** The rank $\text{rank}(\underline{\mathbf{I}}_N, \underline{\mathbf{I}}_N \underline{\mathbf{D}}, \dots, \underline{\mathbf{I}}_N \underline{\mathbf{D}}^{N-2}) = N - 1$ is maximal.
This assumption on the Krylov subspace states that the vectors $\underline{I}_N \underline{D}^i$ are linearly independent and therefore constitute an invertible matrix. Such an assumption is well-known in controllability, for instance, the so-called Kalman's criterion to ensure controllability of a linear system [6]. Another application in the form of Krylov subspaces is used in modern iterative methods for finding one (or a few) eigenvalue(s) of large sparse matrices or for solving large systems of linear equations. This assumption is also related to the Frobenius decomposition with the block companion matrix [6].

**Assumption 4.2.** All the eigenvalues of $\mathbf{K}$ are distinct and none of the last components of the eigenvectors $\mathbf{P}_{ek}$ shall vanish, i.e. $P_{Nk} \neq 0$, for all $k = 1, \dots, N$.

This assumption guarantees the existence of $N$ grazing linear modes [9]. The following lemma shows that for non-trivial solutions of a linear differential system to exist, the multiplicity of the last entry of the solution must be bounded and related to the dimension of the solution.

**Lemma 4.3.** Assumption 4.1 or Assumption 4.2 holds. Let $\mathbf{x}(t) \in \mathbb{R}^N$ be a solution to the linear differential system

$$ \mathbf{M}\ddot{\mathbf{x}} + \mathbf{K}\mathbf{x} = \mathbf{0}. \qquad (4.19) $$

Then, either $\mathrm{mult}(x_N)(0) \le 2N - 1$ or $\mathbf{x} \equiv \mathbf{0}$.

In other words, a non-trivial solution to (4.19) exists only if $\mathrm{mult}(x_N)(t) < 2N$, for all $t$. We first prove this Lemma by using Assumption 4.1; the proof under Assumption 4.2 is provided after.

*Proof.* [With Assumption 4.1] The statement of Lemma 4.3 is equivalent to saying that if $\mathrm{mult}(x_N)(0) > 2N - 1$, that is

$$ x_N(0) = \dot{x}_N(0) = \dots = x_N^{(2N-1)}(0) = 0, \qquad (4.20) $$

then $\mathbf{x} \equiv \mathbf{0}$. From the $N^{\mathrm{th}}$ equation of (4.19), $m_N \ddot{x}_N + \underline{I}_N \mathbf{x} = 0$. Differentiating this equation $2N - 3$ times and using (4.20) implies that the following relations hold

$$ (A_k) : \underline{I}_N \underline{\mathbf{x}}^{(k)}(0) = 0, \quad k = 1, \dots, 2N - 3. \qquad (4.21) $$

The other $N - 1$ equations form a reduced system $\underline{\mathbf{M}}\ddot{\mathbf{x}} + \underline{\mathbf{K}}\mathbf{x} = \underline{\mathbf{c}}$ where $\underline{\mathbf{x}} = (x_j)_{j=1}^{N-1}$ and $\underline{\mathbf{c}} = -x_N \underline{I}_N^T$. Differentiating this system with respect to $t$ gives $\underline{\mathbf{M}}\mathbf{x}^{(3)}(0) + \underline{\mathbf{K}}\dot{\mathbf{x}}(0) = \mathbf{0}$, or $\underline{\mathbf{x}}^{(3)}(0) + \underline{\mathbf{D}}\dot{\mathbf{x}}(0) = \mathbf{0}$. Similarly, differentiating this system $2N - 5$ times yields

$$ (B_k) : \underline{\mathbf{x}}^{(k+2)}(0) + \underline{\mathbf{D}}\underline{\mathbf{x}}^{(k)}(0) = \mathbf{0}, \quad k = 1, \dots, 2N - 5. \qquad (4.22) $$

Since $x_N(0) = \dot{x}_N(0) = 0$, it is sufficient to prove that $\underline{\mathbf{x}}(0) = \mathbf{0}$ and $\dot{\underline{\mathbf{x}}}(0) = \mathbf{0}$. To show the former, a linear system with $\underline{\mathbf{x}}(0)$ as the unknown is constructed as follows. Relation (A$_1$) from (4.21) gives the first equation of the linear system: $\underline{I}_N \underline{\mathbf{x}}(0) = 0$. Multiplying (B$_1$) in (4.22) with $\underline{I}_N$ and using the relation (A$_3$) to eliminate the term $\underline{I}_N \underline{\mathbf{x}}^{(3)}(0)$ results in the second equation: $\underline{I}_N \underline{\mathbf{D}}\underline{\mathbf{x}}(0) = 0$. The third equation is obtained after the following steps. First, multiply (B$_3$) with $\underline{I}_N$, then use relation (A$_5$) to get $\underline{I}_N \underline{\mathbf{D}}\underline{\mathbf{x}}^{(3)}(0) = 0$. This vanishing term appears when (B$_1$) is multiplied by $\underline{I}_N \underline{\mathbf{D}}$, which gives rise to the third equation: $\underline{I}_N \underline{\mathbf{D}}^2 \underline{\mathbf{x}}(0) = 0$. The same recursive process can be used for each $k = 1, \dots, N - 1$: the use of the relations of odd indices $(A_{2k-1}, B_{2k-3}, B_{2k-5}, \dots, B_1)$ gives the $k^{\mathrm{th}}$ equation of the system

$$ \underline{I}_N \underline{\mathbf{D}}^{k-1} \underline{\mathbf{x}}(0) = 0. \qquad (4.23) $$

Combining the $N - 1$ equations (4.23), $k = 1, \dots, N - 1$, it follows that $\underline{\mathbf{x}}(0)$ satisfies $\underline{I}_N \underline{\mathbf{x}}(0) = 0$, $\underline{I}_N \underline{\mathbf{D}}\underline{\mathbf{x}}(0) = 0, \dots, \underline{I}_N \underline{\mathbf{D}}^{N-2}\underline{\mathbf{x}}(0) = 0$. By Assumption 4.1, the unique solution of this linear system is $\underline{\mathbf{x}}(0) = \mathbf{0}$. Together with the hypothesis that $x_N(0) = 0$, then $\underline{\mathbf{x}}(0) = \mathbf{0}$.

Similarly, $\dot{\underline{\mathbf{x}}}(0) = \mathbf{0}$ is shown by constructing another linear system with $\dot{\underline{\mathbf{x}}}(0)$ as the variable. The $k^{\mathrm{th}}$ equation $\underline{I}_N \underline{\mathbf{D}}^{k-1} \dot{\underline{\mathbf{x}}}(0) = 0, k \in [1, \dots, N - 1]$, of the system can be derived by using a similar recursive process now involving the even-indexed relations $(A_{2k}, B_{2k-2}, B_{2k-4}, \dots, B_2)$,
whence, $\underline{I}_N\dot{\underline{x}}(0) = 0$, $\underline{I}_N\underline{D}\dot{\underline{x}}(0) = 0, \dots, \underline{I}_N\underline{D}^{N-2}\dot{\underline{x}}(0) = 0$. It then follows that $\dot{\underline{x}}(0) = \underline{0}$, and together with the assumption $\dot{x}_N(0) = 0$, eventually gives $\dot{x}(0) = \underline{0}$.

It was shown that $\underline{x}(0) = \underline{0}$ and $\dot{\underline{x}}(0) = \underline{0}$. A solution of (4.19) associated to this initial data is then identically zero: $\underline{x}(t) \equiv \underline{0}$. This ends the proof by using Assumption 4.1.

[With Assumption 4.2] Consider the $N^{\text{th}}$ component of $\underline{x}(t)$:

$$ 
\begin{aligned}
x_N(t) &= \underline{e}_N^\top \underline{x}(t) = \underline{e}_N^\top (\underline{\mathbf{P}} \cos (t \underline{\Omega}) \underline{\mathbf{P}}^{-1} \underline{x}(0) + \underline{\mathbf{P}} \underline{\Omega}^{-1} \sin (t \underline{\Omega}) \underline{\mathbf{P}}^{-1} \dot{\underline{x}}(0)) \\
&= \sum_{k=1}^{N} (\alpha_k \cos \omega_k t + \beta_k \sin \omega_k t) \underline{e}_N^\top \underline{\mathbf{P}} \underline{e}_k = \sum_{k=1}^{N} (\alpha_k \cos \omega_k t + \beta_k \sin \omega_k t) \nu_k
\end{aligned} 
\quad (4.24) $$

where $\alpha_k$ and $\beta_k$ are coefficients depending on the initial data $[\underline{x}(0), \dot{\underline{x}}(0)]$, $\nu_k = \underline{e}_N^\top \underline{\mathbf{P}} \underline{e}_k$ are the components of the last row of the matrix $\underline{\mathbf{P}}$.

Assume that $\text{mult}(x_N)(0) > 2N - 1$, using (4.20) for $x_N$ and $N - 1$ first even-order derivatives of $x_N$ yields $\sum_{k=1}^N \alpha_k \nu_k = 0$, $\sum_{k=1}^N \alpha_k \omega_k^2 \nu_k = 0, \dots, \sum_{k=1}^N \alpha_k \omega_k^{2N-2} \nu_k = 0$. Denote $\lambda_k = \omega_k^2$, this can be rewritten as a linear system where $\alpha_k$ are the unknowns:

$$ 
\begin{bmatrix} 1 & 1 & \dots & 1 \\ \lambda_1 & \lambda_2 & \dots & \lambda_N \\ \vdots & \vdots & \ddots & \vdots \\ \lambda_1^{N-1} & \lambda_2^{N-1} & \dots & \lambda_N^{N-1} \end{bmatrix} \begin{bmatrix} \nu_1 & 0 & \dots & 0 \\ 0 & \nu_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \nu_N \end{bmatrix} \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_N \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix}. 
\quad (4.25) $$

Through Assumption 4.2, it follows that $\lambda_k, k = 1, \dots, N$ are distinct, and $\nu_k = P_{Nk} \neq 0$ for all $k = 1, \dots, N$, then there is a unique solution $\alpha_k = 0$, for all $k = 1, \dots, N$.

Similarly, using (4.20) for $N$ odd-order derivatives of $x_N$ yields

$$ \sum_{k=1}^{N} \beta_k \omega_k \nu_k = 0, \quad \sum_{k=1}^{N} \beta_k \omega_k^3 \nu_k = 0, \quad \dots, \quad \sum_{k=1}^{N} \beta_k \omega_k^{2N-1} \nu_k = 0. 
\quad (4.26) $$

This can be recast in a linear system where the $\beta_k$ are the unknowns:

$$ 
\begin{bmatrix} 1 & 1 & \dots & 1 \\ \lambda_1 & \lambda_2 & \dots & \lambda_N \\ \vdots & \vdots & \ddots & \vdots \\ \lambda_1^{N-1} & \lambda_2^{N-1} & \dots & \lambda_N^{N-1} \end{bmatrix} \begin{bmatrix} \omega_1 \nu_1 & 0 & \dots & 0 \\ 0 & \omega_2 \nu_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \omega_N \nu_N \end{bmatrix} \begin{bmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_N \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix}. 
\quad (4.27) $$

With the same argument, the unique solution of this system is $\beta_k = 0$, for all $k = 1, \dots, N$. As a consequence, all the components of $\underline{x}(t)$ vanish, hence, $\underline{x}(t) \equiv \underline{0}$ for all $t$. Lemma 4.3 is then proven. $\square$

This result can be proven to be true for all components $x_k$ of $\underline{x}$, $k = 1, \dots, N$ with suitable Assumption 4.1 or 4.2. Let us apply this general result to the solution to (1.1) in order to show that, at the contact time, $u_N$ has at most $2N - 1$ derivatives which vanish. Between the successive closing contacts, the system is linear and hence the solution $\underline{u}(t)$ is analytic. Therefore, for the sake of simplicity, assume that all the derivatives are taken on the left of 0; $u_N^{(n)-}$ is denoted for the left $n^{\text{th}}$-derivative of $u_N$ at 0 with respect to $t$.

**Proposition 4.4.** Assume that $\underline{u}(t)$ is a solution to (1.1) which has a closing contact at $t = 0$, i.e. $u_N(0) = d$. Then, under Assumption 4.1 or 4.2,

$$ 0 \le \text{mult}(\underline{u}_N^{-})(0) \le 2N - 1. 
\quad (4.28) $$

Moreover, if there is a sticking phase after $t = 0$ of duration $\tau$ and with an assumption similar to Assumption 4.1 or 4.2 for the $(N-1) \times (N-1)$ sticking system rewritten in a suitable basis where $\underline{I}_N\dot{\underline{u}}$ corresponds to the last entry of the solution, then

$$ 0 \le \text{mult}(\underline{I}_N\dot{\underline{u}})(\tau) \le 2N - 3. 
\quad (4.29) $$
*Proof.* The first part of the proposition is obtained by applying Lemma 4.3 to $\mathbf{x} = \dot{\mathbf{u}}^{-}$. Note that the dynamics is linear away from the closing contacts, i.e.

$$ \mathbf{M}\ddot{\mathbf{u}}^{-} + \mathbf{K}\mathbf{u} = \mathbf{0}. \qquad (4.30) $$

Hence, $\dot{\mathbf{u}}$ is the solution to $\mathbf{M}\ddot{\mathbf{u}}^{-} + \mathbf{K}\dot{\mathbf{u}} = \mathbf{0}$. By Lemma 4.3, it follows that $\mathrm{mult}(\dot{u}_{N}^{-})(0) \le 2N - 1$ or $\dot{\mathbf{u}}(t) = \mathbf{0}$. However, if $\dot{\mathbf{u}}(t) = \mathbf{0}$, then substituting this into the linear system (4.30), it follows that $\mathbf{u}(t) = \mathbf{0}$. This contradicts the assumption that $u_{N}(0) = d > 0$.

If there is a sticking phase at $t = 0$, then $\dot{u}_{N}^{-}(0) = 0$ and $\dot{u}_{N}^{-}(\tau) = 0$ where $\tau$ is the end of sticking phase. During the sticking phase from $t = 0$ to $t = \tau$, the $N^{th}$ mass lies on the obstacle, that is

$$ m_N \ddot{u}_N^{-}(t) + \mathbf{l}_N \mathbf{u}(t) = R(t) \qquad (4.31) $$

where $R(t) \le 0$, and $F(t) = \mathbf{l}_N \mathbf{u}(t)$. The duration of the sticking phase $\tau$ is also the end of the sticking phase, it satisfies:

$$ F(\tau) = \mathbf{l}_N \mathbf{u}(\tau) = \mathbf{l}_N \mathbf{u}(\tau) + k_{NN} d = 0 \qquad (4.32) $$

$$ F(\tau + \delta) > 0, \quad \forall \delta \in ]0; \delta_0[, \quad \delta_0 > 0 \qquad (4.33) $$

where $\mathbf{u}$ is the solution of the sticking system $\mathbf{M}\ddot{\mathbf{u}} + \mathbf{K}\dot{\mathbf{u}} = \mathbf{c}$, with $\mathbf{c} = -d\mathbf{l}_N^{\top}$. Differentiating this system with respect to $t$ yields $\mathbf{M}\ddot{\mathbf{u}} + \mathbf{K}\dot{\mathbf{u}} = \mathbf{0}$. A change of variables $\dot{\mathbf{u}} = \mathbf{Q}\mathbf{v}$ using an $\mathbf{M}$-orthogonal matrix $\mathbf{Q}$, i.e. $\mathbf{Q}^{\top}\mathbf{M}\mathbf{Q} = \mathbf{I}$, such that $\mathbf{v}_{N-1} = c\mathbf{l}_N\dot{\mathbf{u}}$, with $c \ne 0$ yields

$$ \ddot{\mathbf{v}} + \mathbf{Q}^{\top}\mathbf{K}\mathbf{Q}\mathbf{v} = \mathbf{0}. \qquad (4.34) $$

Invoking Theorem 4.3 for this $(N - 1) \times (N - 1)$ reduced system, it follows that $\mathrm{mult}(\mathbf{l}_N\dot{\mathbf{u}})(\tau) \le 2N - 3$. $\square$

In the following proposition, the singularity exponents of both the free-flight duration and the sticking phase duration are bounded. These durations have a relationship with the First Return Time. More precisely, if $W \in \mathcal{H}^{-} \cup \mathcal{H}_{G}^{0}$, then the First Return Time of an orbit stemming from W coincides with the duration of the free-flight $s(W)$. Otherwise, if $W \in \mathcal{H}_{S}^{0}$, then the First Return Time involves the duration of sticking phase and the duration of the free-flight: $T(W) = \tau(W) + s(W)$. The bounded multiplicity shown in Proposition 4.4 is the main ingredient to prove the singularity.

**Proposition 4.5** — Power-root singularity. Suppose that $W_0 \in \mathcal{H}_{\rho}$ generates an orbit $[\mathbf{u}, \dot{\mathbf{u}}]$ of (1.1) which has a closing contact at $t = T_0$. Assumption 4.1 or 4.2 and an equivalent assumption for the sticking system holds.

(1) If $W_0 \in \mathcal{H}^{-}$ and $W(T_0) \in \mathcal{H}^{0}$, then the duration of the free-flight of the orbit generated from a nearby W until the next closing contact, denoted by $s(W)$, has an $n^{th}$-root singularity at $W_0$ with

$$ n = 1 + \mathrm{mult}(\dot{u}_{N}^{-})(T_0), \quad 1 \le n \le 2N. \qquad (4.35) $$

(2) Denote $\mathcal{H}_{S}^{0}$ the interior of $\mathcal{H}_{S}^{0}$ in the topology of $\mathcal{H}^{0}$. If $W_0 \in \mathcal{H}_{S}^{0}$, i.e. there is a sticking phase right after $t = 0$ of duration $\tau_0 = \tau(W_0)$, then for all $W \in \mathcal{H}_{S}^{0}$ near $W_0$, the duration of the sticking phase $\tau = \tau(W)$ has an $m^{th}$-root singularity at $W_0$ with

$$ m = 1 + \mathrm{mult}(\mathbf{l}_N \dot{\mathbf{u}})(\tau(W_0)), \quad 1 \le m \le 2N - 2. \qquad (4.36) $$

**Remark 4.3.** The First Return Time can involve two power-root singularities: one from the duration of the sticking phase, and one from the duration of the free-flight phase.

**Remark 4.4.** When $W_0 \in \mathcal{H}_{G}^{0} \cap \bar{\mathcal{H}}_{S}^{0}$, the behaviour of the First Return Time is more complicated. It depends on whether W belongs to $\mathcal{H}_{S}^{0}$ or not. In the former, there are two power-root singularities while in the latter, there is a single power-root singularity.
*Proof.* The idea of the proof is to use the implicit function theorem for the duration of the free-flight or the duration of the sticking phase.

(1) Consider the duration $T_0$ of the free-flight between the two closing contacts at $t = 0$ and $t = T_0$. For $(s, W)$ in the neighborhood of $(T_0, W_0)$, consider the function $f(s, W) = e_N^T R(s) SW$. It is smooth and defined for all arguments and corresponds to the nonsmooth function (due to the closing contacts) $u_N(s, W)$. A key point is to apply the implicit function theorem to $f$ and then to interpret the result for the nonsmooth function $u_N$. By Proposition 4.4, it follows that $0 \le \text{mult}(\dot{u}_N^*)(T_0) \le 2N - 1$. Using Remark 4.2 for the function $f$ in the neighborhood of $(T_0, W_0)$, it follows that there exists a subset containing $W_0$ such that $s$ has an $n^{th}$-root singularity at $W_0$ with $n = 1 + \text{mult}(\dot{u}_N^*)(T_0)$, and thus, $1 \le n \le 2N$.

(2) Recall that $U(\cdot, W) = [u(\cdot, W), \dot{u}(\cdot, W)]$ is the solution to (1.1) associated to the initial data $W$. During the sticking phase, $u_N(t, W) = d$ and $\dot{u}_N(t, W) = 0$ for all $t \in [0; \tau]$. By denoting $F(t, W) = I_N u(t, W) = I_N \underline{u}(t, W) + k_{NN} d$, the last equation of (1.1) yields $m_N \ddot{u}_N(t, W) = -F(t, W) + R(t, W)$ where $R(t, W) \le 0, \forall t \in [0; \tau]$. Hence, the sticking phase starts at $t = 0$ when $F(0, W) = 0$ and $F(\delta, W) < 0, \forall \delta \in ]0; \delta_0[$, $\delta_0 > 0$. The second condition ensures that mass $N$ does not leave the wall immediately. The sticking phase holds as long as $F(t, W) \le 0$, it ends at $t = \tau$ when $F(\tau, W) = 0$ and $F(\tau+\delta, W) > 0, \forall \delta \in ]0; \delta_0[$, $\delta_0 > 0$. The latter condition makes sure that the total force acting on the $N^{th}$ mass becomes strictly negative right after $t = \tau$ and therefore the mass leaves the obstacle. Similarly, denote by $F^*(t) = I_N \underline{u}(t) + k_{NN} d$, where $\underline{u}$ is the solution of the smooth system $M\ddot{\underline{u}} + K\underline{u} = c$, with $c = -dI_N^T$. By Remark 4.2 for the smooth function $F^*$ in the neighborhood of $(\tau_0, W_0)$, it follows that there is a subset of $H_S^0$ containing $W_0$ such that $\tau$ has an $m^{th}$-root singularity at $W_0$ with $m = 1 + \text{mult}(\dot{F})(\tau(W_0))$. Proposition 4.4 implies that $0 \le \text{mult}(I_N \dot{\underline{u}})(\tau(W_0)) \le 2N - 3$, therefore $1 \le m \le 2N - 2$. $\square$

**Remark 4.5.** The lower bound of this multiplicity is optimal for a chain, i.e. $1/2N$ is the greatest lower bound of the power-root singularity of the duration $s$. Let us verify this by showing that, if $u_N(T_0) = d$, $\dot{u}_N(T_0) = \dots = u_N^{(2N-1)}(T_0) = 0$ then there is a unique solution to (1.1) with the maximal $2N$-root singularity of $s$.

Notice that, in this case, the symmetric matrix $K$ has the tridiagonal form

$$ K = \begin{bmatrix} k_{11} & k_{21} & \dots & 0 & 0 \\ k_{21} & k_{22} & \dots & 0 & 0 \\ & & \ddots & & \\ 0 & 0 & \dots & k_{N-1,N-1} & k_{N,N-1} \\ 0 & 0 & \dots & k_{N,N-1} & k_{NN} \end{bmatrix} \quad (4.37) $$

with $k_{i+1,i} \ne 0$ for all $i = 1, \dots, N-1$. Away from closing contacts, $u$ is the solution of the linear system

$$ M\ddot{u} + Ku = 0. \quad (4.38) $$

From the last equation of (4.38), that is $m_N \ddot{u}_N + k_{N,N-1} u_{N-1} + k_{NN} u_N = 0$, together with $k_{N,N-1} \ne 0$, it follows that $u_{N-1}(0) = C_{N-1}d$, $\dot{u}_{N-1}(0) = \dots = u_{N-1}^{(2N-3)}(0) = 0$, and $C_{N-1} = -k_{NN}/k_{N,N-1}$. Similarly, this recursive process can be used until the $(k+1)^{th}$ equation, $k = 2, \dots, N-1$, which, together with the assumption that $k_{k,k-1} \ne 0$, gives

$$ u_k(0) = C_k d \quad \text{and} \quad \dot{u}_k(0) = \dots = u_k^{2k-1}(0) = 0 \quad (4.39) $$

where $C_k = -(k_{kk}C_k + k_{k+1,k}C_{k+1})/k_{k,k-1}$. As a consequence, the initial data is:

$$ u(0) = [C_1 d, \dots, C_{N-1} d, d], \quad \dot{u}(0) = 0. \quad (4.40) $$

Hence, there exists a unique solution to (1.1) associated to this initial data.
However, it is not sure that $1/(2N-2)$ is the greatest lower bound of the singularity of $\tau$ since it depends on the admissibility condition of the end of the sticking phase which is $\mathbf{I}_N\mathbf{u}^{(2N-2)}(\tau) > 0$.

**Remark 4.6** — Power-root singularity for the two-dof chain. Let us clarify the power-root singularity when $N = 2$. In this case, the duration of the free-flight may have a square-root, cube-root, or at most fourth-root singularity, while the duration of the sticking phase is analytic.

An illustration of these power-root singularities is given in Figure 6. Let the state of the system at $t = 0$ be $W_0 = [u_1(0), d, \dot{u}_1(0), \dot{u}_2^-(0)] \in \mathcal{H}_p$. The interesting criteria happens when the first contact $(u_2(T_0) = d)$ is with zero velocity, $\dot{u}_2^-(T_0) = 0$. That means $W(T_0) \in \mathcal{H}^0$.

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Power-root singularity (1/2, 1/3 and 1/4) in the plane $(u_1(T_0), \dot{u}_1(T_0))$ which is isomorphic to the set $\mathcal{H}^0$ since $u_2(T_0) = d$ and $\dot{u}_2(T_0) = 0$. The blue area $u_1 < d$ corresponds to grazing. The red branch $u_1 = d$ and $\dot{u}_1 > 0$ corresponds to the beginning of a sticking phase. The green area $u_1 > d$ corresponds to states within a sticking phase. The solid green line $u_1 = d$ and $\dot{u}_1 < 0$ corresponds to the end of a sticking phase. The dark red dot corresponds to a unique orbit the worst power-root singularity.</figcaption></figure>

(1) If there is no sticking phase after $t = 0$, that is $W_0 \in \mathcal{H}^- \cup \mathcal{H}_G^0$, three configurations exist:

(a) $u_1(T_0) < d$: this implies $m_2\ddot{u}_2^-(T_0) = k_2(u_1(T_0) - u_2(T_0)) < 0$. The duration of the free-flight $s = s(W)$ has a square-root singularity at $W_0$.

(b) $u_1(T_0) = d$ and $\dot{u}_1(T_0) > 0$: this is the beginning of a sticking phase, $\ddot{u}_2(T_0) = 0$ and $m_2\ddot{u}_2^-(T_0) = k_2(\dot{u}_1(T_0) - \dot{u}_2(T_0)) > 0$. Hence, $s = s(W)$ has a cube-root singularity at $W_0$.

(c) $u_1(T_0) = d$ and $\dot{u}_1(T_0) = 0$: this gives $\ddot{u}_2^-(T_0) = \ddot{u}_2^-(T_0) = 0$. However, $m_2u_2^{(4)-}(T_0) = k_2(\ddot{u}_1(T_0) - \ddot{u}_2^-(T_0)) = k_1k_2d > 0$. This case corresponds to a grazing contact and the fourth-root singularity of $s(W)$ arises.

(2) If there is a sticking phase of duration $\tau_0$ at $t = 0$, that is $W_0 \in \mathcal{H}_S^0$, then there are two possibilities: $u_1(0) > d$ or $(u_1(0) = d$ and $\dot{u}_1(0) > 0)$. During a sticking phase, $u_1$ is a solution of $m_1\ddot{u}_1 + (k_1 + k_2)u_1 = k_2d$. The end of the sticking phase is at $t = \tau_0$, when $u_1(\tau_0) = d$ and $\dot{u}_1(\tau_0) < 0$. This shows that $\mathrm{mult}(\dot{u}_1)(\tau_0) \le 2N - 3 = 1$. This is a particular case where the duration of the sticking phase $\tau(W)$ is analytic for $W \in \mathcal{H}_S^0$.

An example of a cube-root singularity is depicted in Figure 7.

4.3. **Discontinuous First Return Time.** The discontinuity of the return time near grazing orbits is a consequence of Theorem 2.3. The simplest case is stated in Corollary 4.1. This result yields a possibly discontinuous First Return Map, higher singularity degree, or instantaneous instability. It is a geometric discontinuity induced by the choice of the Poincaré section. Another choice for the Poincaré section could annihilate this discontinuity.
<figure><img src="image_9.png" /><figcaption>FIGURE 7. First Return Time $T$ with respect to $v_1 = \dot{u}_1$ (near a periodic solution with one sticking phase per period [15]). A cube-root singularity appears near $v_1(0) = \dot{u}_1(0) = 5.86$.</figcaption></figure>

**Corollary 4.1** — **Discontinuous First Return Time.** If $W_0 \in \mathcal{H}^{-}$ and $W(T_0) \in \mathcal{H}^{0}$ where $T_0 = T(W_0) > 0$ and $W(T_0)$ satisfies Assumption 2.2 which says that the grazing contact is not degenerate, then the First Return Time $T(W)$ is discontinuous at $W_0$.

*Proof.* Theorem 2.3 states that $T(W)$ is near $T_0 = T(W_0)$ only on $B_k = \{W \in V_{W_0}, s_k(W_k - \alpha(W)) \ge 0\}$ which is one half of a neighborhood of $W_0$. Accordingly, $T(W)$ is not near $T_0$ for $s_k(W_k - \alpha(W)) < 0$. $\square$

This kind of discontinuity is already reported by Nordmark in [12, p. 282]. His clear explanation is reproduced here: “In the case of mapping between impacts, there is a more fundamental problem due to tangential flow. At $v = 0$, the flow is tangent to the Poincaré section, so that nearby points may impact with low velocity, or may miss and impact at a later time. This is an inherent problem when choosing a section where the flow is not everywhere transversal, and it will cause discontinuities to appear at points being mapped to the line $u = 0$.” Moreover, when Nordmark chose a Poincaré section orthogonal to the contact surface, which corresponds to our Poincaré section $\mathcal{H}_\rho$, he discovered the square-root singularity for a 1-dof model with a source term [12]. This Nordmark map is well summarized in [10, p. 2]. The square-root singularity appears in exactly a half neighborhood which corresponds to our half neighborhood $B_k$. For the other half part of the neighborhood, the Nordmark map is smooth. It corresponds to orbits coming back to the contact surface with a First Return Time which is not almost the time $T_0$ but a later time. It is the reason why our FRT is not continuous.

## 5. SQUARE-ROOT INSTABILITY

In this section, the so-called square-root instability is introduced for a fixed-point of the map. Loosely speaking, the square-root singularity may affect the dynamics and the fixed-point may become unstable.

Let $\mathbf{F}: \mathbb{R}^n \to \mathbb{R}^n$ be a map defined by $\mathbf{F}(\mathbf{X}) = \mathbf{G}(\sqrt{|x_1|}, x_2, \dots, x_n)$ where $\mathbf{G}: \mathbb{R}^n \to \mathbb{R}^n$ is at least a $C^2$ function. $\mathbf{F}$ has a fixed-point $\mathbf{0} \in \mathbb{R}^n$. Consider the dynamical system obtained by iterating $\mathbf{F}$:

$$ \mathbf{X}^{m+1} = \mathbf{F}(\mathbf{X}^m), \quad \mathbf{X}^m = (x_1^m, x_2^m, \dots, x_n^m) \in \mathbb{R}^n, \quad m = 1, 2, \dots \qquad (5.1) $$

With the presence of the square-root term $\sqrt{|x_1|}$, a question on the stability of the fixed-point $\mathbf{0}$ arises. This section provides a generic condition for which the fixed-point $\mathbf{0}$ of the map $\mathbf{F}$ is unstable. Besides, examples show that this condition is only a sufficient condition of the instability of the fixed point even in small dimensions.
**5.1. A nonlinear n-dimensional map.** In the coming theorem, it is shown under a specific condition that the square-root term acting on the component $x_1$ of $\mathbf{X}$ yields the instability of the fixed-point $\mathbf{0}$. Moreover, this instability occurs along the direction of $x_1$.

**Theorem 5.1 — Unstable fixed-point.** Consider two functions $\mathbf{F}$ and $\mathbf{G}: \mathbb{R}^n \to \mathbb{R}^n$ such that $\mathbf{F}(\mathbf{X}) = \mathbf{G}(\mathbf{X})$, where $\mathbf{G}(\mathbf{X}) = (g_1, g_2, \dots, g_n)(\mathbf{X})$, $\mathbf{X} = (x_1, \dots, x_n)$ and $\mathbf{X} = (\sqrt{|x_1|}, x_2, \dots, x_n)$. If $\mathbf{G}$ belongs to $C^2(\mathbb{R}^n, \mathbb{R}^n)$, $\mathbf{G}(\mathbf{0}) = \mathbf{0}$, and the Jacobian $\mathbf{D}\mathbf{G}(\mathbf{0}) = (a_{ij})_{i,j=1}^n$ satisfies

$$a_{11} \neq 0 \qquad (5.2)$$

then $\mathbf{0}$ is an unstable fixed-point of $\mathbf{F}$.

**Remark 5.1.** More generally, if the square-root is activated on the component $x_k$, then the instability of the fixed-point $\mathbf{0}$ occurs if $a_{kk} = \partial_{x_k} g_k(\mathbf{0}) \neq 0$.

*Proof.* Our goal is to show that there is a neighborhood of $\mathbf{0}$ such that many points arbitrarily close to $\mathbf{0}$ will go out of that neighborhood.

The Taylor expansion with an integral remainder of each $g_i$ near $\mathbf{0}$, $i = 1, \dots, n$ is

$$g_i(\mathbf{X}) = a_{i1}\sqrt{|x_1|} + \sum_{j=2}^{n} a_{ij}x_j + |x_1|r_{i,11}(\mathbf{X}) + \\ + 2\sum_{j=2}^{n} \sqrt{|x_1|}x_j r_{i,1j}(\mathbf{X}) + 2\sum_{n \ge k > \ell > 1} x_k x_\ell r_{i,k\ell}(\mathbf{X}) + \sum_{j=2}^{n} x_j^2 r_{i,jj}(\mathbf{X}) \quad (5.3)$$

where $r_{i,k\ell}(\mathbf{X}) = \int_0^1 (1-s)\partial_k\partial_\ell g_i(s\mathbf{X}) \mathrm{d}s$, for all $k, \ell = 1, \dots, n$. On a compact set which will be chosen later, there exists a constant $M > 0$ such that $|r_{i,k\ell}(\mathbf{X})| \le M$ for all $k, \ell = 1, \dots, n$ and for each $i = 1, \dots, n$. Hence,

$$|a_{i1}|\sqrt{|x_1|} - \sum_{j=2}^{n} |a_{ij}||x_j| - M\theta(\mathbf{X}) \le |g_i(\mathbf{X})| \le |a_{i1}|\sqrt{|x_1|} + \sum_{j=2}^{n} |a_{ij}||x_j| + M\theta(\mathbf{X}) \quad (5.4)$$

where $\theta(\mathbf{X}) = |x_1| + 2\sqrt{|x_1|} \sum_{j=2}^n |x_j| + 2\sum_{n \ge k > \ell > 1} |x_k||x_\ell| + \sum_{j=2}^n x_j^2$. A suitable neighborhood, denoted $D_\epsilon$, is constructed so that a sequence starting from any point in $D_\epsilon$ will eventually go away from $\mathbf{0}$ in the direction of $x_1$. To define $D_\epsilon$, the following notations are needed. Let $\alpha_i$, $i = 2, \dots, n$, be:

$$\alpha_i = \begin{cases} \frac{|a_{11}|}{2n|a_{1i}|} & \text{if } a_{1i} \neq 0 \\ 1 & \text{if } a_{1i} = 0. \end{cases} \qquad (5.5a) \\ (5.5b)$$

It follows from this definition that $\alpha_j|a_{1j}| \le |a_{11}|/2n$ while the equality occurs if $a_{1j} \ne 0$; otherwise it is a strict inequality. For $\mathbf{X} \in \mathbb{R}^n$ such that $|x_j| \le \alpha_j\sqrt{|x_1|}$, $\forall j = 2, \dots, n$, we have $|a_{1j}||x_j| \le |a_{1j}|\alpha_j\sqrt{|x_1|}$, and hence,

$$\sum_{j=2}^{n} |a_{1j}||x_j| \le \sum_{j=2}^{n} |a_{1j}|\alpha_j\sqrt{|x_1|} \le (n-1)\frac{|a_{11}|}{2n}\sqrt{|x_1|}. \qquad (5.6)$$

Moreover, the following inequality holds:

$$\theta(\mathbf{X}) \le |x_1| \left( 1 + 2 \sum_{j=2}^{n} \alpha_j + 2 \sum_{n \ge k > \ell > 1} \alpha_k \alpha_\ell + \sum_{j=2}^{n} \alpha_j^2 \right). \qquad (5.7)$$

Denote $\alpha = 1 + 2 \sum_{j=2}^n \alpha_j + 2 \sum_{n \ge k > \ell > 1} \alpha_k \alpha_\ell + \sum_{j=2}^n \alpha_j^2$, the above inequality becomes $\theta(\mathbf{X}) \le \alpha|x_1|$. Denote also $\gamma_0 = a_{11}^2/(2n\alpha M)^2$. Then, for $\mathbf{X} \in \mathbb{R}^n$ such that $|x_1| \le \gamma_0$, we have

$$M\alpha|x_1| \le \frac{|a_{11}|}{2n}\sqrt{|x_1|}. \qquad (5.8)$$
Since $M\theta(\mathbf{X}) \le M\alpha|x_1|$, it follows that

$$M\theta(\mathbf{X}) \le \frac{|a_{11}|}{2n} \sqrt{|x_1|}. \qquad (5.9)$$

Let $C > 0$ be a constant such that

$$C \ge \frac{2}{|a_{11}|} \left( |a_{i1}| + \sum_{j=2}^{n} |a_{ij}|\alpha_j + M\alpha \right), \quad \forall i = 2, \dots, n \qquad (5.10)$$

Let $\gamma > 0$ such that $\gamma = \min_{i=2,...,n}\{\gamma_i > 0 | C\gamma_i \le \alpha_i\sqrt{\gamma_i}\}$. Therefore, for any $\mathbf{X} \in \mathbb{R}^n$ such that $|x_1| \le \gamma$, we have $C|x_1| \le \alpha_i\sqrt{|x_1|}$, for all $i = 2, ..., n$. We can now define $D_\epsilon$ as follows. Let $\epsilon = \min\{\gamma_0, \gamma, a_{11}^2/8\} > 0$. Consider

$$D = \{ \mathbf{X} = (x_1, \dots, x_n) \in \mathbb{R}^n, |x_i| \le C|x_1|, \forall i = 2, \dots, n \} \qquad (5.11)$$

$$D_{\epsilon} = \{ \mathbf{X} = (x_1, \dots, x_n) \in \mathbb{R}^n, |x_1| \le \epsilon \text{ and } |x_i| \le C|x_1|, \forall i = 2, \dots, n \}. \qquad (5.12)$$

This choice of $D_\epsilon$ avoids the criteria in which the instability of $0$ is hidden by starting at a point

<figure><img src="image_10.png" /><figcaption>FIGURE 8. Neighborhood $D_\epsilon$</figcaption></figure>

near $0$ but the sequence comes back at $0$ after one step. See Example 5.1 for more information. An illustration of $D_\epsilon$ is provided in Figure 8. Any $\mathbf{X} \in D_\epsilon$ satisfies the important inequalities as stated next.

**Lemma 5.2.** For any $\mathbf{X} \in D_\epsilon$, the following inequalities hold:

$$\frac{|a_{11}|}{2} \sqrt{|x_1|} \le |g_1(\mathbf{X})| \le \frac{3|a_{11}|}{2} \sqrt{|x_1|}, \qquad (5.13)$$

$$|g_i(\mathbf{X})| \le C_i|g_1(\mathbf{X})|, \quad \forall i = 2, \dots, n. \qquad (5.14)$$

*Proof.* Suppose that $\mathbf{X} \in D_\epsilon$ then $|x_j| \le C|x_1|$ for all $j = 2, ..., n$, and $|x_1| \le \epsilon \le \gamma$. By the definition of $\gamma$, it follows that $C|x_1| \le \alpha_j\sqrt{|x_1|}$. Hence, $|x_j| \le C|x_1| \le \alpha_j\sqrt{|x_1|}$. Thus (5.6) holds. Substituting (5.6) and (5.9) into (5.4) when $i = 1$, the right hand side (RHS) of (5.4) becomes

$$\text{RHS}(5.4) \le |a_{11}|\sqrt{|x_1|} + (n-1)\frac{|a_{11}|}{2n}\sqrt{|x_1|} + \frac{|a_{11}|}{2n}\sqrt{|x_1|} \le \frac{3|a_{11}|}{2}\sqrt{|x_1|}. \qquad (5.15)$$

Similarly,

$$\text{LHS}(5.4) \ge |a_{11}|\sqrt{|x_1|} - (n-1)\frac{|a_{11}|}{2n}\sqrt{|x_1|} - \frac{|a_{11}|}{2n}\sqrt{|x_1|} \ge \frac{|a_{11}|}{2}\sqrt{|x_1|}. \qquad (5.16)$$
Therefore, $\mathbf{X}$ satisfies (5.13).

Let us prove inequality (5.14). For $i = 2, \dots, n$, expressions (5.4) and (5.9) imply

$$|g_i(\underline{\mathbf{X}})| \le \left( |a_{i1}| + \sum_{j=2}^{n} |a_{ij}|\alpha_j + M\alpha \right) \sqrt{|x_1|}. \quad (5.17)$$

Via inequality (5.13), it is shown that $|a_{11}|\sqrt{|x_1|} \le 2|g_1(\tilde{\mathbf{X}})|$. Hence, to prove (5.14), it is sufficient to show that

$$\frac{2}{|a_{11}|} \left( |a_{i1}| + \sum_{j=2}^{n} |a_{ij}|\alpha_j + M\alpha \right) \le C. \quad (5.18)$$

This is true by the choice of $C$ given in (5.10). Inequality (5.14) and Lemma 5.2 are then proven. $\square$

Back to the proof of Theorem 5.1, the idea is to show that the recurrence goes away from $\mathbf{0}$ in the direction of the first component. In other words, the square-root singularity acting on the first component plays an important role via the inequalities (5.13) and (5.14).

Let us show that, for any $0 < \delta < \epsilon$, if the sequence $(\mathbf{X}^m)_{m \ge 1}$ in $D_\epsilon$, where $\mathbf{X}^m = (x_1^m, \dots, x_n^m)$, is defined by

$$\mathbf{X}^0 = (\delta/2, 0, \dots, 0) \in D_\epsilon, \quad \mathbf{X}^{m+1} = \mathbf{F}(\mathbf{X}^m), \quad m \ge 1, \quad (5.19)$$

then there exists a $N_0 > 0$ such that $\mathbf{X}^{N_0} \notin D_\epsilon$. From (5.14), $\mathbf{X}^1 = \mathbf{F}(\mathbf{X}^0) = \mathbf{G}(\underline{\mathbf{X}}^0) \in D$. If $|x_1^1| > \epsilon$ then $\mathbf{X}^1 \notin D_\epsilon$. This shows the instability of $\mathbf{0}$. Otherwise, if $|x_1^1| \le \epsilon$ and $|x_1^1| = |g_i(\mathbf{X}^0)| \le C|x_1^1|$ by inequality (5.14), then $\mathbf{X}^1 \in D_\epsilon$, and $\mathbf{X}^2$ is considered. If there exists $\mathbf{X}^m \in D_\epsilon$ for all $m \ge 1$ then

$$\frac{a_{11}}{2} \sqrt{|x_1^m|} \le |x_1^{m+1}| \le \frac{3a_{11}}{2} \sqrt{|x_1^m|}. \quad (5.20)$$

Consider the sequence $z_m$ defined by $z_{m+1} = a_{11}/2\sqrt{|z_m|}$, $z_0 = x_1^0 = \delta/2 \le \epsilon \le a_{11}^2/8$. Then, inequality (5.20) yields $|x_1^m| \ge z_m$. It is known that $z_m$ increasingly converges to $a_{11}^2/4$ in an interval $(0, a_{11}^2/4)$. Hence, there exists $N_0 > 0$ such that $|z_{N_0}| > a_{11}^2/8 \ge \epsilon$, and thus $|x_1^{N_0}| \ge |z_{N_0}| > \epsilon$. That means, $\mathbf{X}^{N_0} \notin D_\epsilon$ and hence $\mathbf{0}$ is unstable. $\square$

**Example 5.1.** Consider $\mathbf{F}(x, y) = \mathbf{G}(\sqrt{|x|}, y)$ when $\mathbf{G}$ is a linear map

$$\mathbf{G}(\mathbf{X}) = \mathbf{G}(x, y) = \mathbf{A}\mathbf{X} = \begin{bmatrix} a & b \\ c & d \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix}, \quad a \neq 0. \quad (5.21)$$

Figure 9(a) shows that the fixed-point $(0, 0)$ of $\mathbf{F}$ is unstable and the recurrence goes away from this point along the line $y = cx/a$. An interesting case is when $c = \alpha a$, $d = \alpha b$ with $\alpha \neq 0$. The instability of $(0, 0)$ is hidden if the sequence starts at a point belonging to the curve (C): $y = -a\sqrt{|x|}/b$ because the sequence stays at $(0, 0)$ as soon as the second step, as illustrated in Figure 9(b). To show the instability, it is required to start at a point $(x_0, y_0)$ which does not lie on C. This is the reason why the set $D_\epsilon$ is chosen as specified in the previous proof.

5.2. **Two-dimensional maps with critical instability.** As proven in Theorem 5.1, the square-root instability of the fixed-point appears when $a_{11} \neq 0$. When $a_{11}$ vanishes, the square-root singularity may appear under additional assumptions on the other components of $\mathbf{A}$.

**Proposition 5.3.** Consider two functions $\mathbf{F}$ and $\mathbf{G}: \mathbb{R}^2 \to \mathbb{R}^2$ such that $\mathbf{F}(x, y) = \mathbf{G}(\sqrt{|x|}, y)$ where $\mathbf{G}$ is a linear map satisfying $\mathbf{G}(x, y) = (by, cx + dy)$ with $bc \neq 0$. Then, $(0, 0)$ is an unstable fixed-point of $\mathbf{F}$.
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Instability of the fixed-point (0,0) for the map (5.21). (a) $a = 1, c = 2, b = d = 0$: the recurrence goes away from (0,0) along the line $y = cx/a$; (b) $c = \alpha a, d = \alpha b$: the instability does not realize if the starting point sits on the curve $C: y = -a\sqrt{|x|}/b$ since the next iterate is (0,0) (c) $b = c = 1, a = d = 0$: iterates leave (0,0). The gradient color scale [color bar] shows initial iterates in blue to final iterates in red irrespective of the magnitude.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 10. Instability of the fixed-point (0,0) for the map (5.21). (a) when $0 < d = 0.5 < 1$ and $c = 1$; (b) when $-1 < d = -0.5 < 0$ and $c = 1$; (c) when $d = 0$ and $c = 1$. The gradient color scale [color bar] shows initial iterates in blue to final iterates in red irrespective of the magnitude.</figcaption></figure>

*Proof.* By considering the map $(\mathbf{F} \circ \mathbf{F})(x, y) = (bc\sqrt{|x|} + bdy, cd\sqrt{|x|} + d^2y + c\sqrt{|by|})$, it is clear that $a_{11} = bc \neq 0$, then the technique in the proof of Theorem 5.1 is used to show the instability of the fixed-point (0,0) of $\mathbf{F} \circ \mathbf{F}$.

□

Let us denote

$$ \mathbf{DG}(0,0) = \begin{bmatrix} 0 & b \\ c & d \end{bmatrix}. \qquad (5.22) $$

Then, in the case of a linear map $\mathbf{G}$, the fixed-point (0,0) of $\mathbf{F}$ is stable if and only if $bc = 0$ and $|d| < 1$. The square-root term in $x$ disappears because of the condition $bc = 0$ and the stability analysis of the origin for the two-dimensional map reduces to the stability analysis of the origin for the one-dimensional linear map, and hence $|d| < 1$ is needed to guarantee stability.

**Example 5.2.** A simple example of the instability mechanism is shown in Fig. 9(c) where $b = c = 1$ and $a = d = 0$. Starting at a point closed to (0,0), the recurrence goes away from this point.

**Example 5.3.** Is also of interest the situation where (0,0) is linearly stable but nonlinearly unstable. Consider $\mathbf{F}_1(x, y) = (0, \sqrt{|x| + dy})$ where $|d| < 1$. Thus (0,0) is a stable fixed-point of $\mathbf{F}_1$. However, for the map $\mathbf{F}_2(x, y) = (\sqrt{|x|y}, \sqrt{|x|} + dy)$, the fixed-point (0,0) is numerically shown to be unstable.

More precisely, three cases are expected:
(1) for $0 < d < 1$, the recurrence oscillates during several steps then goes away from $(0,0)$ along the curve similar to $y = k\sqrt{|x|}$, see Figure 10(a).

(2) for $-1 < d < 0$, the recurrence seems to always oscillate around $(0,0)$ (see Figure 10(b)).

(3) For $d = 0$ and when $G$ is linear, the fixed-point $(0,0)$ for $F$ is asymptotically stable. when $G$ is nonlinear, the problem is open, see Figure 10(c).

## 6. IN THE VICINITY OF A GRAZING ORBIT

The aim of this section is first to prove the existence of the square-root singularity near a grazing contact, that is to prove Theorems 2.3 and 2.4. Then, the square-root singularity near periodic solutions with one grazing contact per period is investigated. Such periodic solutions of the $N$-dof system are limited to $N$ grazing linear modes and there is no other such linear solutions. Periodic solutions with sticking contacts [15, 17] are excluded. Take note that $H_p$ ensures that the First Return Time exists and is finite. A challenging point in building the FRT is to check the nonlocal admissibility condition $u_N(t) < d$ before the closing contact. Moreover, data may lead to a FRT which is not in the vicinity of $T_0$ but instead in the vicinity of $2T_0, 3T_0$ and so on. We pay attention to the class of initial data near $W_0$ in $H^-$ or $H_G^0$ which leads to orbits having their FRT near $T_0$.

**6.1. First Return Time.** This section deals with the proof of Theorem 2.3. It will be shown that, in a certain class of initial data, the FRT has a particular form containing a square-root term. First, a lemma is stated, then a long proof is proposed. Recall that $\Phi(t, W) = e_N^T R(t) SW$ is smooth. It coincides with $u_N$ in the neighborhood of $(T_0, W_0)$ as long as $u_N(t) \le d$.

**Lemma 6.1.** The set $\mathcal{K}$ of indices $\{i \in \{N+1, \dots, 2N\}, \partial_{W_i} \Phi(T_0, W_0) \neq 0\}$ is not empty for $N \ge 2$.

In other words, there exists at least one non-vanishing partial derivative of $\Phi$ with respect to at least one initial velocity. This guarantees the existence of a component $W_k$ stated in Theorem 2.3.

*Proof.* This is proved by contradiction. The smooth function $\Phi(t, W)$ can be written as

$$ \Phi(t, W) = e_N^T R(t) SW = e_N^T (P \cos(t\Omega)P^{-1}u + P\Omega^{-1} \sin(t\Omega)P^{-1}v) \quad (6.1) $$

If $\mathcal{K} = \emptyset$ then $\partial_{W_i} \Phi(T_0, W_0) = 0$ for all $i = N+1, \dots, 2N$. It follows that

$$ e_N^T P \Omega^{-1} \sin(T_0 \Omega) P^{-1} e_i = 0 \text{ for all } i = 1, \dots, N. \quad (6.2) $$

In other words, the $N^{th}$ row of $A = P\Omega^{-1} \sin(T_0\Omega)P^{-1}$ is identically zero. Thus, $\sin(T_0\Omega) = 0$, i.e. $\sin \omega_i T_0 = 0, \forall i = 1, \dots, N$ which implies $\omega_i T_0 = k_i \pi, k_i \in \mathbb{N}$: this is impossible since $N \ge 2$, there is no internal resonance through Assumption 2.1. Hence, $\mathcal{K} \neq \emptyset$. $\square$

Now, the implicit equation $\Phi(T, W) = d$ is locally solved via a function $T = \theta(W)$ defined on half a neighborhood of $W_0$. Ultimately, it is proven that $\theta(W)$ is the First Return Time $T(W)$.

*Proof.* Consider the initial condition $W_0 \in H^- \cup H_G^0$ at initial instant $t = 0$, the associated orbit shows a first grazing contact at $T_0 = T(W_0)$. The First Return Time $T$ is implicitly given by the equation $u_N(T, W) = d$. The nonsmooth function $u_N$ is replaced by the smooth function

$$ \Phi(T, W) = d \quad (6.3) $$

on which the implicit function theorem applies and where $\Phi(T_0, W_0) = u_N(T_0, W_0) = d$, $\partial_t \Phi(T_0, W_0) = \dot{u}_N(T_0, W_0) = 0$. As a consequence of Lemma 6.1, there exists $k \in \mathcal{K}$ where $\partial_{W_k} \Phi(T_0, W_0) \neq 0$. Denote $\underline{W}$ is the reduced vector obtained from $W$ by removing $W_k$. Together with Assumption 2.2, function $\Phi(t, W)$ can be seen as $\Phi(t, W_k, \underline{W})$ where $W_0$ corresponds to $(W_{0k}, \underline{W}_0)$, and $f$ satisfies the following conditions: $\Phi(T_0, W_{0k}, \underline{W}_0) = d$, $\partial_t \Phi(T_0, W_{0k}, \underline{W}_0) = 0$,
$\partial_{W_k} \Phi(T_0, W_{0k}, \underline{W}_0) \neq 0$, $\partial_t^2 \Phi(T_0, W_{0k}, \underline{W}_0) \neq 0$. Via Theorem 4.2, there exist neighborhoods $V_{\underline{W}_0}$, $V_{T_0}$ and $V_{W_{0k}}$ of $\underline{W}_0$, $T_0$ and $W_{0k}$ respectively and smooth scalar functions $\eta, \alpha$ such that

$$ \eta: V_{\underline{W}_0} \to V_{T_0}, \underline{W} \mapsto T = \eta(\underline{W}) \qquad (6.4) $$

$$ \alpha: V_{\underline{W}_0} \to V_{W_{0k}}, \underline{W} \mapsto W_k = \alpha(\underline{W}) \qquad (6.5) $$

satisfying $\eta(\underline{W}_0) = T_0$ and $\alpha(\underline{W}_0) = W_{0k}$. The set $S_c$ defined as the intersection of the two hypersurfaces $\Phi = d$ and $\partial_t \Phi = 0$ can be parameterized as

$$ \begin{aligned} S_c &= \{ (t, W_k, \underline{W}) \in V_{T_0} \times V_{W_{0k}} \times V_{\underline{W}_0}, \Phi(t, W_k, \underline{W}) = d \text{ and } \partial_t \Phi(t, W_k, \underline{W}) = 0 \} \\ &= \{ (\eta(\underline{W}), \alpha(\underline{W}), \underline{W}), \underline{W} \in V_{\underline{W}_0} \}. \end{aligned} \qquad (6.6) $$

Recall that $\mathcal{B}_k = \{W = (W_k, \underline{W}) \in V_{W_{0k}} \times V_{W_{0k}}, s_k(W_k - \alpha(\underline{W})) \geq 0\}$ where $s_k = \text{sign } \gamma_k$. It is the region adjacent to and including the hypersurface $W_k = \alpha(\underline{W})$. By applying Theorem 4.2, there exists a smooth function $\psi$ such that there are two graphs solving equation (6.3):

$$ T = \theta(W) = \eta(\underline{W}) + \psi(\pm\sqrt{s_k(W_k - \alpha(\underline{W})), \underline{W}}), \quad W \in \mathcal{B}_k \qquad (6.7) $$

where $\psi(0, \underline{W}_0) = 0$ and $\partial_{W_k} \psi(0, \underline{W}_0) = |\gamma_k|^{-1/2}$. We choose the branch of $\theta(W)$ corresponding to the admissibility condition for the velocity at the contact so that the mass comes in contact in the right side of the wall, which is

$$ \dot{\Phi}(\theta(W), W) \geq 0. \qquad (6.8) $$

Denote $F(W) = R(\theta(W))SW$. Using the asymptotic expansion up to the first order of $\sqrt{s_k(W_k - \alpha(\underline{W}_0))}$ in the direction of the $k^{th}$ component of W, and taking note that $\eta(\underline{W}_0) = T_0$, one has

$$ \theta(W) = T_0 \pm |\gamma_k|^{-1/2} \sqrt{s_k(W_k - \alpha(\underline{W}_0))} + O(h^2) = T_0 + h + O(h^2) \qquad (6.9) $$

where $h = \pm|\gamma_k|^{-1/2}\sqrt{s_k(W_k - \alpha(\underline{W}_0))}$ from which we deduce that $W_k = \alpha(\underline{W}_0) \pm \gamma_k h^2 = W_{0k} + O(h^2)$ and $W = W_0 + h^2 e_k$, $e_k^\top = [0, \dots, 1, \dots, 0] \in \mathbb{R}^{2N}$ follows, since $\alpha(\underline{W}_0) = W_{0k}$, see Theorem 4.2. Similarly,

$$ \cos(\theta(W) \Omega) = \cos((T_0 + h + O(h^2)) \Omega) = \cos(T_0 \Omega) - h \Omega \sin(T_0 \Omega) + O(h^2) $$

$$ \sin(\theta(W) \Omega) = \sin((T_0 + h + O(h^2)) \Omega) = \sin(T_0 \Omega) + h \Omega \cos(T_0 \Omega) + O(h^2) $$

and $R(\theta(W))$ can be written as $R(T_0 + h + O(h^2)) = R(T_0) + h\dot{R}(T_0) + O(h^2)$ where $\dot{R}(t)$ denotes the matrix whose elements are the derivatives of the elements of $R$:

$$ \dot{R}(t) = \begin{bmatrix} -\mathbf{P} \Omega \sin(t \Omega) \mathbf{P}^{-1} & \mathbf{P} \cos(t \Omega) \mathbf{P}^{-1} \\ -\mathbf{P} \Omega^2 \cos(t \Omega) \mathbf{P}^{-1} & -\mathbf{P} \Omega \sin(t \Omega) \mathbf{P}^{-1} \end{bmatrix}. \qquad (6.10) $$

Inserting these expressions into $F$ and neglecting higher order terms in $h$ leads to

$$ F(W) \approx (R(T_0) + h\dot{R}(T_0))SW_0 \approx \begin{bmatrix} \mathbf{u}(T_0, W_0) \\ \mathbf{v}(T_0, W_0) \end{bmatrix} + h \begin{bmatrix} \mathbf{v}(T_0, W_0) \\ -\mathbf{Ku}(T_0, W_0) \end{bmatrix}. \qquad (6.11) $$

Thus, $\dot{\Phi}(\theta(W), W) = e_N^\top F(W) \approx e_N^\top (\mathbf{v}(T_0, W_0) - h\mathbf{Ku}(T_0, W_0)) = -h e_N^\top \mathbf{Ku}(T_0, W_0)$ where $e_{2N}$ is the vector in $\mathbb{R}^{2N}$ with all zero components except the $2N^{th}$ one equal to unity, in a way similar to the vector $e_N$ in $\mathbb{R}^N$. From the $N^{th}$ equation of (1.1) taken at $t = T_0$, that is $m_N \ddot{u}_N^-(T_0) + e_N^\top \mathbf{Ku}(T_0) = 0$, it follows that $m_N \ddot{u}_N^-(T_0, W_0) = -e_N^\top \mathbf{Ku}(T_0, W_0)$ and thus

$$ \dot{\Phi}(\theta(W), W) \approx hm_N \ddot{u}_N^-(T_0, W_0). \qquad (6.12) $$

Since $\ddot{u}_N^-(T_0, W_0) < 0$, $h$ is chosen with the negative sign in order to have $\dot{\Phi}(\theta(W)) \geq 0$. Hence, $\theta(W)$ satisfies (6.8). $\square$

The proof of Theorem 2.3 ends by showing that $\theta(W) = T(W)$. This means that $\Phi(t, W) < d$ for all $t \in ]0; \theta(W)[$.
*Proof.* First, $\ddot{\Phi}(T_0, W_0) < 0$ and $\dot{\Phi}(T_0, W_0) = 0$ imply that $\Phi(t, W_0) < d$ for $t < T_0$ and $t$ near $T_0$. By the smoothness of $\Phi$ with respect to $(t, W)$, there exists a neighborhood $[T_-; T_+] \times V_1$ of $(T_0, W_0)$ and $\delta_1 > 0$ such that $\ddot{\Phi}(T, W) < -\delta_1 < 0$ in this neighborhood and $\theta(W) \in [T_-; T_+]$ for $W \in V_1$. This will be the crucial ingredient to conclude at the end that $\Phi(t, W) < d$ for all $t \in ]T_-; \theta(W)]$.

Second, $\dot{\Phi}(0, W_0) < 0$ because $W_0$ belongs to $\mathcal{H}^-$. Thus, there exists a neighborhood $[0; T_2] \times V_2$ of $(0, W_0)$ with $V_2 \subset V_1$ and $\delta_2 > 0$ such that $\dot{\Phi}(T, W) < -\delta_2 < 0$ in this neighborhood. This implies that $\Phi(t, W) < d - \delta_2 t < d$ for all $t \in ]0; T_2]$. Let $\delta_3 = 2^{-1}\delta_2 T_2 > 0$, there exists a neighborhood of $W_0$ denoted by $V_3 \subset V_2$ such that $|\Phi(t, W) - \Phi(t, W_0)| < \delta_3$ for all $(t, W) \in [T_2; T_-] \times V_3$. This yields $\Phi(t, W) < d - \delta_3 < d$ for all $(t, W) \in [T_2; T_-] \times V_3$.

Finally, on $[T_-; \theta(W)]$ with $W \in V_3$, the function $\Phi(t, W)$ is concave and the velocity at time $\theta(W)$ is nonnegative from (6.8) and $\Phi(t, W) < d$ on $[T_-; \theta(W)] \times V_3$ which concludes the proof for $\theta(W) = T(W)$. This also ends the proof of Theorem 2.3. $\square$

The proof of Theorem 2.4 is similar in many aspects except two points: the neighborhood of $W_0$ is smaller and the proof that $\Phi(t, W) < d$ for $0 < t$ small enough requires a careful inspection. The difference is that the initial velocity of the last mass is zero. However, the acceleration is negative $\ddot{u}_N^+(0, W_0) < 0$ which is enough to ensure that $\Phi(t, W) < d$ for $0 < t$ small enough. This concludes the proof of Theorem 2.4.

**6.2. Square-root singularity in the vicinity of a grazing linear mode.** This subsection explores the possible square-root singularity near any GLM. One feature of a GLM is that the sticking phase does not occur near such a motion. This property is proven in the coming Proposition. Then, the square-root singularity coefficients are computed. At least one of them, non-vanishing, will activate the square-root singularity near a GLM. If the First Return Map has a particular expression with the square-root term in a class of initial data then one may expect the instability of any GLM as in Section 5. Such an instability near periodic grazing solutions is expected [3] but is not yet mathematically proven.

Throughout this section, let us consider the $j^{th}$ GLM (the notation $\text{GLM}_j$ is also used below) with period $T_j = T(W_0)$ [9] associated to the initial data

$$W_0 = [\mathbf{u}_0^\top, \mathbf{v}_0^\top]^\top, \quad \mathbf{u}_0 = \frac{d}{P_{Nj}} \mathbf{P} \mathbf{e}_j, \quad \mathbf{v}_0 = \mathbf{0}. \qquad (6.13)$$

**Proposition 6.2.** *The sticking phase does not occur near the $j^{th}$ GLM, $j = 1, \dots, N$.*

*Proof.* From the $N^{th}$ equation in (1.1), $m_N \ddot{u}_N(t) + \mathbf{e}_N^\top \mathbf{K} \mathbf{u}(t) = R(t)$, $R(t) \le 0$, the sticking phase does not occur if $u_N(t) = d$ and $F(t) = \mathbf{e}_N^\top \mathbf{K} \mathbf{u}(t) > 0$. By the periodicity of GLM, it is sufficient to show that $F(0) > 0$. The initial data of the GLM yields

$$F(0) = \frac{d}{P_{Nj}} \mathbf{e}_N^\top \mathbf{K} \mathbf{P} \mathbf{e}_j = \frac{d}{P_{Nj}} \mathbf{e}_N^\top \mathbf{M} \mathbf{P} \mathbf{\Omega}^2 \mathbf{P}^{-1} \mathbf{P} \mathbf{e}_j = d m_N \omega_j^2 > 0 \text{ for } d > 0. \qquad (6.14)$$

which ends the proof. $\square$

In some cases when $d > 0$ and for some initial data in $\mathcal{B}_k^+$ defined in Theorem 2.4, the associated orbit takes less time to come back to the Poincaré section. This is a consequence of Theorem 2.4 and a feature of GLM.

**Corollary 6.1.** *Consider a GLM with period $T_j$ and an initial perturbation of its initial data $W = W_0 + w \mathbf{e}_k$ with $k \ne 2N$ (or $k = 2N$ and $\sigma > 0$ where $\sigma = \text{sign } \gamma_k$), then there exists $\delta > 0$ such that $0 < \sigma w < \delta$ and $T(W) < T_j$.*

Such one-sided condition on the First Return Time is already known for nonlinear modes with one impact per period near a GLM [9]. It is expected that this inequality is valid on a larger set near $W_0$ but of course not in the whole neighborhood if $k \ne 2N$.
<figure><img src="image_13.png" /><figcaption>FIGURE 11. First Return Time (red lines) with respect to the initial displacement of the first mass: (a) near the first GLM, (b) near the second GLM.</figcaption></figure>

*Proof.* To apply Theorem 2.4, Assumption 2.2 must be verified. The second time derivative of $u_N$ at the grazing point is

$$ \ddot{u}_N(T_j, W_0) = -\frac{1}{m_N} \mathbf{e}_N^\top \mathbf{K} \mathbf{u}(T_j) = -\frac{1}{m_N} \mathbf{e}_N^\top \mathbf{K} \mathbf{u}(0) = -d\omega_j^2 < 0. \quad (6.15) $$

Hence, by Theorem 2.3, the First Return Time $T$ takes the form given by (2.16). Moreover, since $\ddot{u}_N(T_0, W_0) < 0$, it follows that in the direction of the $k^{th}$ component of $W$

$$ T = \eta(\underline{W}_0) + \psi(-\sqrt{\sigma w}, \underline{W}_0), \quad \eta(\underline{W}_0) = T_j, \quad \partial_{W_k} \psi(0, \underline{W}_0) = 1/|\gamma_k|. \quad (6.16) $$

Hence, $T = T_j - 1/|\gamma_k|\sqrt{\sigma w} + O(w) < T_j$. $\square$

Let us name $u_{1GLM_1}(0)$ and $u_{1GLM_2}(0)$, the initial displacement of the first mass along $GLM_1$ and $GLM_2$, respectively. The square-root singularity of the FRT with respect to the initial $u_1$, when $u_1(0) \ge u_{1GLM_1}(0)$ or $u_1(0) \le u_{1GLM_2}(0)$ is shown in Figure 11. The red lines are computed numerically, while the dashed lines illustrate the Taylor expansion of $\psi$. In both cases, $T < T_j, j = 1, 2$.

The following is the proof of Theorem 2.5 which gives the computation of the coefficients $C_k$ and the generic condition that causes the square-root singularity for the First Return Map near GLMs.

*Proof.* Via Proposition 6.2, the sticking phase does not occur near GLMs, thus the First Return Map takes the form (2.19a) $\mathcal{F}(W) = \mathbf{R}(T(W))\mathbf{S}\mathbf{W}$ where the FRT has the square-root dependence (2.16) in the subset $B_k^+$. Using the proof in Subsection 6.1, the First Return Map can be rewritten as $\mathcal{F}(W) = \mathbf{G}(\bar{W}) = (g_1, \dots, g_{2N})(\bar{W})$ where $\bar{W} = [\bar{W}_i]_{i=1}^{2N}$, with a change of variables:

$$ \bar{W}_k = \sqrt{s_k(W_k - \alpha(\underline{W}_0))} \quad \text{and} \quad \bar{W}_i = \underline{W}_i, \quad \forall i \ne k. \quad (6.17) $$

Note that $g_N(\bar{W}) = d$. From the proof in Subsection 6.1, it follows that $\bar{W}_k = \sqrt{|\gamma_k|}h$. Together with (6.11), $\mathcal{F}$ can then be written as

$$ \mathcal{F}(W) = \mathbf{G}(W_1, \dots, \sqrt{|\gamma_k|}h, \dots, W_{2N}) \approx \begin{bmatrix} \mathbf{u}(T_j, W_0) \\ \mathbf{v}(T_j, W_0) \end{bmatrix} + h \begin{bmatrix} \mathbf{v}(T_j, W_0) \\ -\mathbf{K} \mathbf{u}(T_j, W_0) \end{bmatrix}. \quad (6.18) $$

For $k \in \{1, \dots, 2N\}, k \ne N$ such that $\partial_{W_k} u_N(T_j, W_0) \ne 0$, the coefficient associated with the square-root term $h$ is $C_k = a_{kk} = \partial_{\bar{W}_k} g_k(T_j, W_0) = \sqrt{|\gamma_k|} \partial_h g_k(T_j, W_0)$. The expression of $C_k$ is
then obtained via (6.18) as

$$C_k = \begin{cases} \sqrt{|\gamma_k|} \mathbf{e}_k^\top \mathbf{v}(T_j, \mathbf{W}_0) & \text{if } 1 \le k < N, \\ -\sqrt{|\gamma_k|} \mathbf{e}_{k-N}^\top \mathbf{K} \mathbf{u}(T_j, \mathbf{W}_0) & \text{if } N < k \le 2N. \end{cases} \quad (6.19a)$$

$$C_k = \begin{cases} -\sqrt{|\gamma_k|} \mathbf{e}_{k-N}^\top \mathbf{K} \mathbf{u}(T_j, \mathbf{W}_0) & \text{if } N < k \le 2N. \end{cases} \quad (6.19b)$$

Since $\mathbf{v}(T_j) = \mathbf{v}_0 = \mathbf{0}$, it follows that $C_k = 0, \forall 1 \le k < N$. Another way to write $C_k$ is by using the modal coordinates $[\mathbf{q}, \dot{\mathbf{q}}]$ where $\mathbf{u} = \mathbf{P}\mathbf{q}$ as

$$C_k = \begin{cases} \sqrt{|\gamma_k|} \mathbf{e}_k^\top \mathbf{P} \dot{\mathbf{q}}(T_j, \mathbf{W}_0) & \text{if } 1 \le k < N, \\ -\sqrt{|\gamma_k|} \mathbf{e}_{k-N}^\top \mathbf{M} \mathbf{P} \mathbf{\Omega}^2 \mathbf{q}(T_j, \mathbf{W}_0) & \text{if } N < k \le 2N. \end{cases} \quad (6.20a)$$

$$C_k = \begin{cases} -\sqrt{|\gamma_k|} \mathbf{e}_{k-N}^\top \mathbf{M} \mathbf{P} \mathbf{\Omega}^2 \mathbf{q}(T_j, \mathbf{W}_0) & \text{if } N < k \le 2N. \end{cases} \quad (6.20b)$$

From Lemma 6.1, $\mathcal{K} \neq \emptyset$ and for every $k \in \mathcal{K}$, the coefficient $C_k$ is

$$C_k = -\sqrt{|\gamma_k|} \mathbf{e}_{k-N}^\top \mathbf{M} \mathbf{P} \mathbf{\Omega}^2 \mathbf{q}(T_j) = -\sqrt{|\gamma_i|} dm_i \omega_j^2 P_{kj} / P_{Nj}. \quad (6.21)$$

By the hypothesis that there exists $i \in \mathcal{K}$ where $P_{ij} \neq 0$, it follows that $C_i \neq 0$. This non-vanishing coefficient then facilitates the analysis of the square-root singularity for the First Return Map near GLMs. $\square$

The square-root singularity for the First Return Map is activated near $W_0$ and in a particular class of initial data. If the orbits stay in that regime of initial data, the dynamics will follow the framework of Section 5 and one may be able to determine the instability of GLMs.

**Remark 6.1.** In this section, we have only considered a First Return Time near the fundamental period $T_j$ of the GLM $j$. A similar study can be done for any multiple of this period $T_\ell = \ell T_j$, $\ell \in \mathbb{N}^*$. This yields exactly the same coefficient $C_k$.

## 7. CONCLUSION

A wealth of studies are available on one-dof vibro-impact oscillators with nonlinear ODEs and source terms [2]. In the present article, we consider $N$-dof autonomous vibro-impact systems with linear ODEs. The choice of the Poincaré section $\mathcal{H}_\mathcal{P}$ is critical for the qualitative exploration of the associated nonsmooth dynamics. It is a subset of the contact hyperplane $\mathcal{H} = \{u_N = d\} \subset \mathbb{R}^{2N}$. Indeed, the only nonlinear and nonsmooth behavior of orbits occurs only on this hyperplane and it seems natural to define it as the Poincaré section. However, it features a geometric drawback since it is not transverse to all orbits. It appears that this serious drawback is directly related to the richness of the nonsmooth dynamics of vibro-impact systems. The difficulty is with orbits involving a closing contact with a vanishing incoming velocity. Such orbits, initially investigated by Nordmark [12], can generate a grazing bifurcation and a square-root singularity for the First Return Time and the First Return Map. It is also known that many nonlinear modes emerge from the grazing linear modes. The present paper pays attention to the initial data on $\mathcal{H}$ corresponding to grazing orbits.

We define the Poincaré section $\mathcal{H}_\mathcal{P}$ as the largest subset of $\mathcal{H}$ where the orbits emitted from this subset return to $\mathcal{H}$ infinitely many times. Unfamiliar orbits hitting $\mathcal{H}$ only once are discovered. These very specific grazing orbits hit the subset $\mathcal{H}_1^0 \subset \mathcal{H}$. We also prove that an orbit with a sticking contact always has a finite sticking duration (when $d > 0$—no pre-stress—otherwise it can be false) and crosses $\mathcal{H}_\mathcal{P}$. This means that if an orbit has a sticking contact, then the orbit crosses $\mathcal{H}_\mathcal{P}$ infinitely many times.

Our work extensively analyzes the First Return Time fundamental function defined from $\mathcal{H}_\mathcal{P}$ to $]0; +\infty[$. By construction of $\mathcal{H}_\mathcal{P}$, the First Return Time is a well-defined function. We revisit the square-root singularity on the chosen unusual Poincaré section in a self-contained and elementary analytical framework. We also prove that stronger singularities exist. For two-dof systems, the data producing the cube-root or quartic-root singularity is explicitly obtained. Generically, for $N$-dof systems, it is shown that the order of the singularity is at most $2N$. Another striking result says that The First Return Time is not continuous over a whole neighborhood of data producing a
grazing contact. This strange new property is caused by orbits which do not return to $\mathcal{H}$ with almost the same time. In short, they can return to $\mathcal{H}$ after more of fewer cycles than the grazing orbit does.

The square-root singularity of the First Return Time may affect the First Return Map. We introduce the square root singularity coefficient $C_k$ for $k = 1, \dots, 2N$. When one $C_k \neq 0$ then the First Return Map also exhibits a square root singularity. We prove that for a vectorial map defined on an open set with such a square-root singularity at a fixed point, then the point is unstable. We cannot directly apply this result to the First Return Map because our study of the square-root singularity of the First Return Time is not valid over an entire neighborhood. Accordingly, the expected instability of the Grazing Linear Modes is still an open problem.

## APPENDIX A. ONE-DOF SYSTEM: ON THE SQUARE-ROOT SINGULARITY

Consider the one-dof dynamics depicted in Figure 12:

$$m\ddot{u} + ku = r \qquad (A.1)$$

$$u(0) = u_0, \quad \dot{u}(0^-) = \dot{u}_0^- \qquad (A.2)$$

$$u(t) \le d, \quad r(t) \le 0, \quad (d - u(t))r(t) = 0 \qquad (A.3)$$

$$\dot{u}^+ = -\dot{u}^- \quad \text{when} \quad u(t) = d \qquad (A.4)$$

The solution to (A.1) is always periodic with or without the obstacle and it takes the form

$$u(t) = u_0 \cos \omega t + \dot{u}_0^- / \omega \sin \omega t \qquad (A.5)$$

on $]0; T_0[$ where $\omega = \sqrt{k/m}$ is the linear frequency and $T_0$ is the period of the nonlinear solution and also the First Return Time. It can be expressed as

<figure><img src="image_14.png" /><figcaption>FIGURE 12. Simple unilaterally constrained one-dof system</figcaption></figure>

$$T_0 = T - 2\tau \qquad (A.6)$$

where $T = 2\pi/\omega$ and $2\tau$ is the complementary time for the mass to go from the state $(d, +\dot{u}_0^-)$ back to $(d, -\dot{u}_0^-)$ if there was no obstacle, see Figure 13. Since both linear and nonlinear

<figure><img src="image_15.png" /><figcaption>FIGURE 13. One-dof system orbit</figcaption></figure>

solutions are periodic, $T$ and $T_0$ are exactly the First Return Time to the Poincaré section $\{[u_0, \dot{u}_0^-] \in \mathbb{R}^2, u_0 = d, \dot{u}_0^- > 0\}$ of the linear/nonlinear orbit. In order to approximate $\tau$, a non-admissible solution is considered, see the gray branch in Figure 13. Suppose that the initial
condition is not on the Poincaré section: $u(0) = d(1 + \epsilon)$, $\epsilon > 0$ and $\dot{u}(0) = 0$. The solution is then $u(t) = d(1 + \epsilon) \cos \omega t$ and since at $t = \tau$, $u(\tau) = d$, it follows that $d(1 + \epsilon) \cos \omega \tau = d$. A Taylor expansion in the neighborhood of $\tau = 0$ and $\epsilon = 0$ reads

$$1 - \frac{\omega^2 \tau^2}{2} + O(\tau^4) = 1 - \epsilon + O(\epsilon^2) \qquad (A.7)$$

$$\omega^2 \tau^2 = 2(\epsilon + O(\epsilon^2) + O(\tau^4)). \qquad (A.8)$$

Accordingly, $\omega\tau = \sqrt{2\epsilon + O(\epsilon^2)} = \sqrt{2\epsilon} + O(\epsilon^{3/2})$ and the period of the nonlinear solution is

$$T_0 = T - 2\tau = T - \frac{2}{\omega}\sqrt{2\epsilon + O(\epsilon^2)} = T - \frac{2}{\omega}\sqrt{2\epsilon} + O(\epsilon^{3/2}). \qquad (A.9)$$

Thus, the square-root singularity emerges in the FRT. However, this is slightly misleading because the derivations proposed in the paper are expressed in terms of $\dot{u}^-$. In this one-dof example, the total energy reads $\omega^2 d^2 + (\dot{u}_0^-)^2 = \omega^2 d^2(1 + \epsilon)^2$ thus $\epsilon = (\dot{u}_0^-)^2/(2\omega^2 d^2) + O((\dot{u}_0^-)^2)$ and the FRT becomes:

$$T_0 = T - \frac{2\dot{u}_0^-}{d\omega^2} + O((\dot{u}_0^-)^2) \qquad (A.10)$$

keeping in mind that $W^T = (W_1, W_2) = (u_0, \dot{u}_0^-) = (d, \dot{u}_0^-)$. Another reason why there is no square-root singularity for the FRT is that Assumption 2.3 is not verified. The details are as follows: $\Phi(t, W) = \Phi(t, W_1, W_2) = \Phi(t, u_0, \dot{u}_0^-) = u_0 \cos \omega t + \dot{u}_0^- / \omega \sin \omega t$. However, $\partial_{W_2} \Phi = \dot{u}_0^- = 0$ for the sole grazing orbit. Accordingly, Assumption 2.3 which reduces to $\partial_{W_2} \Phi \neq 0$ is not satisfied.

Also, in this one-dof model, the FRM is the identity function $F(U) = U, \forall U \in H_p$ and hence it is smooth with respect to the initial data.

## APPENDIX B. POWER-ROOT SINGULARITY FOR A MASS-SPRING CHAIN

In the case of a chain, there is a simpler proof for Proposition 4.4 as stated below.

**Lemma B.1.** Suppose $u(t)$ is a solution to (1.1) which models a chain of $N$ masses and has a closing contact at $t = 0$. Under the assumption that $k_{j+1,j} \neq 0$ for all $j = 1, \dots, N-1$, then $\text{mult}(\dot{u}_N^-)(0) \le 2N-1$. Moreover, if there is sticking phase at $t = 0$, then $\text{mult}(\dot{u}_{N-1})(\tau) \le 2N-3$.

*Proof.* The proof includes two parts. The first part shows that $\text{mult}(\dot{u}_N^-)(0) \le 2N-1$. Otherwise, $\text{mult}(\dot{u}_N^-)(0) \ge 2N$, i.e. $\dot{u}_N^- (0) = \dots = u_N^{(2N-1)-}(0) = u_N^{(2N)-}(0) = 0$ and it is shown that $u(0) = 0$ and $\dot{u}(0) = 0$.

Let $v = \dot{u}$, then outside the closing contacts, $v$ is the solution of the linear system

$$M\ddot{v} + Kv = 0 \qquad (B.1)$$

Moreover, $v_N = \dot{u}_N$, from the assumption, $v_N$ satisfies $v_N^{(\ell)}(0) = 0$ for all $1 \le \ell \le 2N-1$. From the last equation of (B.1): $m_N \dot{v}_N + k_{N,N-1} v_{N-1} + k_{NN} v_N = 0$, and $k_{N,N-1} \neq 0$, it follows that $v_{N-1}(0) = \dots = v_N^{(2N-3)}(0) = 0$. Similarly, from the $(N-1)^{\text{th}}$ equation and by the assumption that $k_{N-1,N-2} \neq 0$, this yields $v_{N-2}(0) = \dots = v_{N-2}^{(2N-5)}(0) = 0$. In the end, the second equation and $k_{21} \neq 0$ give $v_1(0) = \dot{v}_1(0) = 0$. As a consequence, we have $v(0) = 0, \dot{v}(0) = 0$ and the corresponding solution is $v = 0, \dot{v} = 0$. Thus, $\dot{u} = 0$ and $\ddot{u} = 0$. Substitution of these identities into (1.1) yields $Ku = 0$. This induces $u = 0$ since $\det(K) \neq 0$. However, this contradicts the fact that $u_N(0) = d > 0$, i.e. the solution cannot rest at its equilibrium $0$.

The next part shows that $\text{mult}(\dot{u}_{N-1})(\tau) \le 2N-3$ when a sticking phase arises after $t = 0$. Based on the last equation of (1.1), $m_N \ddot{u}_N + k_{N,N-1} u_{N-1} + k_{NN} u_N = 0$, the end of the sticking phase at $t = \tau$ satisfies

$$k_{N,N-1} u_{N-1}(\tau) = -k_{NN} d \qquad (B.2)$$

$$k_{N,N-1} u_{N-1}(\tau + \delta) > -k_{NN} d, \quad \forall 0 < \delta < 1 \qquad (B.3)$$
The sticking system is then a ($N-1$) degree-of-freedom non-homogeneous system $\underline{Mü} + \underline{Ku} = \underline{c}$ where $\underline{c} = [0, \dots, 0, k_{N-1,N}d]^T$. Now, assume that $\text{mult}(\dot{u}_{N-1})(\tau) \ge 2N-2$, then $\dot{u}_{N-1}(\tau) = \dots = u^{(2N-2)}(\tau) = 0$. By a similar argument for the reduced ($N-1$) $\times$ ($N-1$) system with the last entry of the solution being $u_{N-1}$ instead of $u_N$, it follows that $\text{mult}(\dot{u}_{N-1})(\tau) \le 2N-3$. $\square$

REFERENCES

[1] Patrick Ballard. "The Dynamics of Discrete Mechanical Systems with Perfect Unilateral Constraints". *Archive for Rational Mechanics and Analysis*, 2000, 154(3):199-274. doi: 10.1007/s002050000105. [hal-00111308].

[2] Mario di Bernardo, Chris Budd, Alan Champneys, and Piotr Kowalczyk. *Piecewise-smooth dynamical systems: theory and applications*. Applied Mathematical Sciences. Springer Science & Business Media, 2008. ISBN: 9781846280399.

[3] Chris Budd and Felix Dux. "Intermittency in impact oscillators close to resonance." *Nonlinearity*, 1994, 7(4):1191-1224. doi: 10.1088/0951-7715/7/4/007. [hal-02520814].

[4] David Chillingworth. "Dynamics of an impact oscillator near a degenerate graze". *Nonlinearity*, 2010, 23(11):2723-2748. doi: 10.1088/0951-7715/23/11/001. [hal-01390258].

[5] Constantin Corduneanu. *Almost periodic functions*. Interscience tracts in pure and applied mathematics. Interscience Publishers, 1968.

[6] Jean-Michel Coron. *Control and Nonlinearity*. Vol. 136. Mathematical Surveys and Monographs. American Mathematical Society, 2007. ISBN: 9780821849187.

[7] Mats Fredriksson and Arne Nordmark. "Bifurcations caused by grazing incidence in many degrees of freedom impact oscillators". *Proceedings of the Royal Society*, 1997, 453(1961):1261-1276. doi: 10.1098/rspa.1997.0069. [hal-01297285].

[8] Stéphane Junca, Huong Le Thi, Mathias Legrand, and Anders Thorin. "Impact dynamics near unilaterally constrained grazing orbits". *9th European Nonlinear Dynamics Conference (ENOC)*. Budapest, Hungary, 2017. [hal-01562154].

[9] Mathias Legrand, Stéphane Junca, and Sokly Heng. "Nonsmooth modal analysis of a $N$-degree-of-freedom system undergoing a purely elastic impact law". *Communications in Nonlinear Science and Numerical Simulation*, 2017, 45:190-219. doi: 10.1016/j.cnsns.2016.08.022. [hal-01185980].

[10] Jaap Molenaar, John G. de Weger, and Willem van de Water. "Mappings of grazing-impact oscillators". *Nonlinearity*, 2001, 14(2):301-321. doi: 10.1088/0951-7715/14/2/307. [hal-02520794].

[11] Arne Nordmark. "Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators". *Nonlinearity*, 2001, 14(6):1517-1542. doi: 10.1088/0951-7715/14/6/306. [hal-01297283].

[12] Arne Nordmark. "Non-periodic motion caused by grazing incidence in an impact oscillator". *Journal of Sound and Vibration*, 1991, 145(2):279-297. doi: 10.1016/0022-460X(91)90592-8. [hal-02546986].

[13] Michelle Schatzman. "A class of nonlinear differential equations of second order in time". *Nonlinear Analysis: Theory, Methods & Applications*, 1978, 2(3):355-373. doi: 10.1016/0362-546X(78)90022-6. [hal-01294058].

[14] Jorge Sotomayor and Marco Antonio Teixeira. "Vector fields near the boundary of a 3-manifold". *Dynamical systems Valparaíso 1986*. Vol. 1331. Lecture Notes in Mathematics. 1988:169-195. doi: 10.1007/BFb0083073. [hal-02558845].

[15] Huong Le Thi, Stéphane Junca, and Mathias Legrand. "Periodic solutions of a two-degree-of-freedom autonomous vibro-impact oscillator with sticking phases". *Nonlinear Analysis: Hybrid Systems*, 2018, 28:54-74. doi: 10.1016/j.nahs.2017.10.009. [hal-01305719v2].

[16] Anders Thorin, Pierre Delezoide, and Mathias Legrand. "Nonsmooth modal analysis of piecewise-linear impact oscillators". *SIAM Journal on Applied Dynamical Systems*, 2017, 16(3):1710-1747. doi: 10.1137/16M1081506. [hal-01298983].

[17] Anders Thorin, Pierre Delezoide, and Mathias Legrand. "Periodic solutions of $n$-dof autonomous vibro-impact oscillators with one lasting contact phase". *Nonlinear Dynamics*, 2017, 90(3):1771-1783. doi: 10.1007/s11071-017-3763-z. [hal-01505888].
[18] Anders Thorin, Mathias Legrand, and Stéphane Junca. "Nonsmooth modal analysis: investigation of a 2-dof spring-mass system subject to an elastic impact law". *Proceedings of the ASME IDETC & CIEC: 11th International Conference on Multibody Systems, Nonlinear Dynamics, and Control*. Boston, Massachusetts, 2015. doi: 10.1115/DETC2015-46796. [hal-01185973].

[19] Phanikrishna Thota. "Analytical and Computational Tools for the Study of Grazing Bifurcations of Periodic Orbits and Invariant Tori". PhD thesis. Engineering Sciences [physics], Virginia Polytechnic Institute and State University, 2007. [tel-01330429].