## 2 Methods
The spectrum of the Lyapunov exponent for the linearized Equation (9) is estimated through [PERSON]'s algorithm using the QR-factorization orthonormalization [6; 24]. In this model, Equation (6) represents a nonsmooth dynamical system with discontinuous impact effects. The dynamic system, including the discontinuity at the moment of impact, can be rewritten from [PERSON]'s method. The discontinuity occurs at the point of impact (\(t=t_{i}\)). \[t_{i-1}<t<t_{i}\text{: }\dot{\mathbf{z}}=\mathbf{f}_{i}(\mathbf{z}),\ \mathbf{z}(t_{i-1})=\mathbf{z}(t_{i-1}^{+}) \tag{10}\]  
\[t=t_{i}\text{: }0=\mathbf{h}(\mathbf{z}(t_{i}^{-  
\[\widetilde{t}_{i}<t<\widetilde{t}_{i+1}\text{: }\dot{\widetilde{\mathbf{z}}}= \mathbf{f}_{i+1}(\widetilde{\mathbf{z}})\text{, }\widetilde{\mathbf{z}}(\widetilde{t}_{i})= \widetilde{\mathbf{z}}(\widetilde{t}_{i}^{+}) \tag{19}\]  
where each interval of discontinuities is smooth, and \(\mathbf{h}(\mathbf{z})\) and \(\mathbf{g}(\mathbf{z})\) are the indicator function and the transition condition, respectively. The plus and minus denote the right- and left-sided limits, and  
\[\delta t_{i}=\widetilde{t}_{i}-t_{i}=-\frac{D\mathbf{h}(\mathbf{z}_{i}^{-}) \delta\mathbf{z}_{i}^{-}}{D\mathbf{h}(\mathbf{z}_{i}^{-})\mathbf{f}_{i}( \mathbf{z}_{i}^{-})} \tag{20}\]  
\[\delta\mathbf{z}_{i}^{+}=D\mathbf{g}(\mathbf{z}_{i}^{-})\delta\mathbf{z}_{i} ^{-}+[D\mathbf{g}(\mathbf{z}_{i}^{-})\mathbf{f}_{i}(\mathbf{z}_{i}^{-})- \mathbf{f}_{i+1}(\mathbf{z}_{i}^{+})]\delta t_{i} \tag{21}\]  
in which  
\[D\mathbf{h}(\mathbf{z}_{i}^{-})=\left.\frac{\partial\mathbf{h}(\mathbf{z})}{ \partial\mathbf{z}}\right|_{\mathbf{z}=\mathbf{z}_{i}^{-}},\ D\mathbf{g}( \mathbf{z}_{i}^{-})=\left.\frac{\partial\mathbf{g}(\mathbf{z})}{\partial \mathbf{z}}\right|_{\mathbf{z}=\mathbf{z}_{i}^{-}} \tag{22}\]  
are the Jacobian matrix of the indicator function and transition condition at point \(\mathbf{z}_{i}^{-}\), respectively, and \(\mathbf{z}_{i}^{-}=\mathbf{z}(t_{i}^{-})\) and \(\mathbf{z}_{i}^{+}=\mathbf{z}(t_{i}^{+})\). For the impact oscillator with [PERSON]'s model of contact, the Jacobian matrix of transition condition and indicator function becomes the following matrix:  
\[D\mathbf{h}(\mathbf{z}_{i}^{-})=\left[ \begin{array}{cccc}1&0&0&\end{array}  \right]^{\mathbf{T}},\ D\mathbf{g}(\mathbf{z}_{i}^{-})=\mathbf{I}. \tag{23}\]  
The recurrence plot (RP) is used as a tool to visualize the recurrence pattern of a dynamical system. A recurrence matrix is expressed by a series of solution vectors. The corresponding RP is based on the following recurrence matrix as follows:  
\[\mathbf{R}_{i,j}=\|\mathbf{z}_{i}-\mathbf{z}_{j}\|,i,j=1,\ldots,N \tag{24}\]  
where \(\|\cdot\|\) is an L-2 norm, \(N\) is the measured points, and \(\{\mathbf{z}_{i}\}_{i=1}^{N}\) is a system's trajectory in its phase space. Here, the element of phase space indicates the possible state of the system for the time evolution law. In such cases, phase space needs to be reconstructed. The most popular method of reconstruction is the time delay method  
\[\mathbf{z}_{i}^{\text{A}}=\sum_{j=1}^{m}u_{i+(j-1)\
u}\mathbf{e}_{j} \tag{25}\]  
where \(u_{i}=u(i\Delta t)\) is the discrete-time series, \(\Delta t\) is the sampling rate, \(m\) is the embedding dimension, \(\
u\) is the time delay, and \(\mathbf{e}_{j}\) are unit vectors. The reconstruction does not change the dynamical properties, and the reconstructed phase space can be expressed through an appropriately selected embedding dimension and time delay. In general, the time delay can be determined appropriately using the Mutual information method. During time-delay reconstruction, all self-crossing of the trajectory for the attractor's dimension \(d_{A}\) can disappear if it is set to the embedding dimension \(d>2d_{A}\). From a physical point of view, it is very important to determine the minimum embedding dimension to minimize the numerical calculation of the Lyapunov exponent.