## 2 Materials and Methods
### Stability Analysis of the Self-Excited System  
Stability analysis allows us to determine whether the system's quasi-static operating point is stable or unstable and, in the latter case, what the unstable modes are. However, in order to perform a stability analysis based on the study of the complex eigenvalues, it is necessary to linearize the equations governing the nonlinear dynamic behavior of the system (cf. Equation (6)) around the quasi-static equilibrium point (noted with subindex \"\(st\)\"). To do so, the friction coefficient \(\mu(s,N)\) is linearized around \((\alpha,N)\) using a first-order [PERSON] approximation, as follows:  
\[\mu(s,N) \approx\mu(\alpha,N)+\frac{d\mu(s,N)}{ds}|_{st}(s-\alpha)\] \[=\mu_{st}+\frac{d\mu(s,N)}{ds}|_{st}\frac{\dot{u}(t)}{V_{x}} \tag{11}\]  
Replacing this in Equation (6) and applying relation (9) in order to express velocity \(\dot{u}\) in terms of \(\dot{\mathbf{q}}\), leads to:  
\[\mathbf{\mu}\ddot{\mathbf{q}}(t)+\mathbf{c}\dot{\mathbf{q}}(t)+\mathbf{\gamma}\mathbf{q}(t )=\frac{N}{V_{x}}\bigg{(}\frac{d\mu(s,N)}{ds}|_{st}\bigg{)}\mathbf{\Phi}^{T}\mathbf{ \Phi}\dot{\mathbf{q}}(t) \tag{12}\]  
The linearized damping matrix \(\mathbf{c}_{lin}\) is defined such that:  
\[\mathbf{c}_{lin}=\mathbf{c}-\frac{N}{V_{x}}\bigg{(}\frac{d\mu(s,N)}{ds}|_{st} \bigg{)}\mathbf{\Phi}^{T}\mathbf{\Phi} \tag{13}\]The equation governing the system's nonlinear behavior (6) is linearized around the quasi-static equilibrium such that:  
\[\mu\ddot{\mathbf{q}}(t)+\mathbf{c}_{lin}\dot{\mathbf{q}}(t)+\gamma\mathbf{q}(t)= \mathbf{0} \tag{14}\]  
The associated eigenvalues problem satisfies the following equation, where \(\hat{\mathbf{q}}\) represents the complex mode shapes:  
\[\left(\lambda^{2}\mu+\lambda\mathbf{c}_{lin}+\gamma\right)\hat{\mathbf{q}}= \mathbf{0} \tag{15}\]  
This quadratic eigenvalue problem is reduced into the generalized eigenvalue problem, as follows:  
\[(\Lambda\mathbf{A}+\mathbf{B})\hat{\mathbf{\theta}}=\mathbf{0} \tag{16}\]  
with  
\[\mathbf{A}= \begin{bmatrix}\mathbf{c}_{lin}&\mu\\ \mu&\mathbf{0}\end{bmatrix} \qquad\mathbf{B}= \begin{bmatrix}\gamma&\mathbf{0} \\ \mathbf{0}&-\mu\end{bmatrix} \qquad\hat{\mathbf{\theta}}= \begin{bmatrix}\hat{ \mathbf{q}}\\ \lambda\hat{\mathbf{q}}\end{bmatrix} \]  
Vector \(\hat{\mathbf{\theta}}\) corresponds to mode shapes in state space. Because of the original eigenvalue problem being quadratic, numerical solution of Equation (16) gives two eigenvalues \(\lambda_{i_{1,2}}\) per mode. Furthermore, the instabilities associated to a decaying friction coefficient with respect to creepage are translated into an equivalent negative damping (hence the linearized damping matrix \(\mathbf{c}_{lin}\) in Equation (13)). Therefore, a convenient way to assess the degree of instability of mode \(i\) is to analyze the associated equivalent damping \(\xi_{i,eq}\), which is proportional to its growth rate. The following formulation is proposed to express eigenvalues \(\lambda_{i_{1,2}}\) as a function of the equivalent damping \(\xi_{i,eq}\) and the natural frequency \(\omega_{i}\) of mode \(i\):  
\[\lambda_{i_{1,2}}= \begin{cases}-\xi_{i,eq}\omega_{i}\pm i\omega_{i}\sqrt{1- \xi_{i,eq}^{2}}&\text{if }\xi_{i,eq}\leq 1\\ -\xi_{i,eq}\omega_{i}\pm\omega_{i}\sqrt{\xi_{i,eq}^{2}-1}&\text{otherwise}\end{cases}  \tag{17}\]  
By analogy to a positive growth rate \(\Re(\lambda_{i})>0\), a negative equivalent damping indicates that mode \(i\) is unstable. From now on, when referring to the most unstable mode, understand the mode for which \(\Re(\lambda_{i})\) is the highest, or as an equivalent, \(\xi_{i,eq}\) is the most negative. However, it can be pointed out here that the mode with the strongest real part is not necessarily the one that will contribute most to the generation of self-sustained vibrations. Indeed, linear conditions used in the stability analysis are not valid during transient oscillations.