**1. Introduction**

The dynamics of systems affected by friction are most often studied in the context of friction-excited vibrations (FIV). Prominent examples for FIV in mechanical structures and machines range from brake systems [1–4], clutches [5], drill strings [6] to artificial hip joints [7] and others. FIV often arise through positive energy feedback from a friction interface with the structure, i.e., through self-excitation [8–10]. Sub-critical Hopf bifurcations [11,12] and isolated solution branches [13–15] are a common observation in those systems, such that bi- and multi-stable systems have been reported numerously [14,16,17]. The computation of those nonlinear responses (periodic, quasi-period orbits, chaotic trajectories) is a well-established field of research [18–21], mostly resulting in the identification of complicated bifurcation diagrams [11,13,22–24]. The stability of the solutions is usually assessed by local Lyapunov-type stability metrics [25,26]. Hence, the stability statement is often a binary one that measures the state's robustness against *small* perturbations. However, the actual size of *permissible perturbations*, i.e., those for which the trajectory would still return back to the state, is not given. In a multi-stable system configuration, the long term steady-state behavior thus depends on the choice of initial conditions or the size of instantaneous perturbations. Once the system enters another basin of attraction, severe jumping phenomena may occur. Typically, such jumps are related to a change from a stable steady sliding state to high-intensity periodic vibrations or stick-slip cycles [27–29], or from one periodic solution to another periodic solution [11].

This work investigates a rather novel technique denoted as *basin stability* to estimate the size of the system's basins of attraction in a subset of the state space. The basins' size estimation can then be considered a global stability metric, i.e., indicating how likely the system is to end up on one of the co-existing stable steady-states. Therefore, those probabilities add new insights to the rather binary stability statements derived from local perturbation-based approaches. We study a friction oscillator excited by a falling friction slope and a second oscillator excited through binary flutter instability. Our results indicate that the basin stability analysis is a robust and easily applicable model-agnostic technique. It can reveal the *actual* picture of the long-term behavior for a given set of perturbations, thus augmenting classical bifurcation and stability studies. Using the basin stability analysis, some solutions can even be ruled out if one can guarantee strict control over the instantaneous perturbations to system trajectories or operating conditions.

#### **2. The Concept of Basin Stability**

We study nonlinear dynamic systems

$$
\dot{\mathbf{x}} = \mathbf{f}\left(\mathbf{x}, t\right), \quad \mathbf{x} \in \mathbb{R}^N \tag{1}
$$

with the states **x** (*t*) in the *N*-dimensional state space. The long-term asymptotic behavior is denoted as *attractor* A [30] throughout this work. Typically, the Lyapunov spectrum <sup>Λ</sup> = [*λ*1, . . . , *<sup>λ</sup>N*] is assessed to characterize the linear stability of a state **x** against small perturbations. For fixed points, the Lyapunov exponents are equivalent to the system's eigenvalues derived from the complex eigenvalue analysis (CEA). The real parts of the eigenvalues indicate linear stability to a *small* perturbation about the fixed point. The sizes of the real parts indicate the strength of attraction (*λ* < 0) or rejection (*λ* > 0) for stable or unstable directions in state space, respectively. However, the eigenvalues do not encode a piece of information about the permissible size of perturbations that are still attracted by the fixed point. While this is not an issue for systems that feature only a single stable solution, the situation is different for systems featuring multiple stable solutions. For these systems, local stability concepts may have only a limited validity: a non-small perturbation of a state can result in a jump to another attractor. Hence, global stability concepts are required to assess the size of permissible perturbations, i.e., to characterize the basins of attraction for all solutions. The basin of attraction

$$\mathcal{B}\left(\mathcal{A}\right) = \left\{ \mathbf{x}\_0 \in \mathbb{R}^N \mid \lim\_{t \to \infty} \mathbf{x}\left(t\right) = \mathcal{A}, \quad \mathbf{x}\_0 = \mathbf{x}\left(t = 0\right) \right\} \tag{2}$$

denotes the subset of states that converge to the same attracting set A. The basin boundaries are related to unstable solutions of the system which represent separatrices of the basins in state space. Depending on the size and shape of its basin, an attractor can be more or less robust against non-small perturbations. There are multiple ways to compute the basins of attraction, e.g., through Lyapunov functions [31]. These methods are known for some canonical, low-dimensional, and well-studied systems. However, they are not readily available, or straight-forward to compute, for any generic and high-dimensional nonlinear dynamical system, such as frictional oscillators which are studied in this work.

The basin stability proposed by Menck et al. [32,33] is a global stability concept for complex systems that aims at measuring stability against non-small perturbations by a volume-based probabilistic approach. Conceptually, the basin stability measures the volumetric share of all basins of attraction in a hypervolume of the state space. For a computationally feasible solution, a distribution *<sup>ρ</sup>* (**x**) of perturbations is drawn from a reference subset Q ⊂ <sup>R</sup>*<sup>N</sup>* of the state space, representing a set of states to which the system may be pushed to through non-small perturbations with R Q *ρ* (**x**) *d***x** = 1. For each perturbation, the steady-state behavior of the dynamical system is obtained through time-marching integrations. Then, the fraction of perturbed states that converged to the specific attractor A denotes an estimate for the basin stability S<sup>B</sup> (A), i.e., [32,34]

$$\mathcal{S}\_{\mathcal{B}}\left(\mathcal{A}\right) = \int \kappa\_{\mathcal{B}}\left(\mathcal{A}\right)\left(\mathbf{x}\right)\rho\left(\mathbf{x}\right)d\mathbf{x}, \quad \kappa\_{\mathcal{B}}\left(\mathcal{A}\right)\left(\mathbf{x}\right) = \begin{cases} 1, & \text{if } \mathbf{x} \in \mathcal{B}\left(\mathcal{A}\right) \\ 0, & \text{otherwise} \end{cases}.\tag{3}$$

Here, *κ*<sup>B</sup> denotes an indicator function that classifies a steady state solution **x** (*t*) to belong to the attractor A. Therefore, S<sup>B</sup> (A) is an estimate for the volume share of the basin of attraction B<sup>A</sup> given the reference subset Q ⊂ <sup>R</sup>*<sup>N</sup>* sampled by *<sup>ρ</sup>* (**x**) [32,35]. Naturally, for a *<sup>k</sup>*-multi-stable system, the basin stability values of all *k* attractors add up to unity ∑ *k <sup>i</sup>* S<sup>B</sup> (A*i*) = 1. The size, i.e., the number of states, of the dynamical system to be studied by the basin stability is only limited by computational power for the Monte Carlo simulations. As the basin stability computation can be considered a repeated Bernoulli experiment [32], the standard error of the basin stability estimate is

$$e = \frac{\sqrt{\mathcal{S}\_{\mathcal{B}} \left(1 - \mathcal{S}\_{\mathcal{B}}\right)}}{\sqrt{n}} \tag{4}$$

which can be used to find a subset Q that ensures a low standard error. Recently, systems with fractal, riddled, and intermingled basin boundaries were studied [33] indicating the robustness of the basin stability concept. All basin stability computations in this work were obtained from the open-source package bSTAB [36] available at https://github.com/TUHH-DYN/bSTAB/tree/v1.0.

Figure 1 displays a schematic for illustrating the practical computation of basin stability values. A nonlinear dynamical system with two states **x** = [*x*1, *x*2] is studied (In fact, the system is the single-degree-of-freedom frictional oscillator to be discussed in Section 3). The system exhibits three solutions: A stable equilibrium position (**x**EP), an unstable periodic orbit (**x**UPO), and a stable limit cycle (**x**LC). The distribution of perturbations *ρ* (**x**) is chosen such that all solutions are contained in Q and *n* = 100 samples are drawn uniformly at random. The trajectories starting from *n*EP = 37 states in the basin BEP converge towards the equilibrium position, while *n*LC = 63 states are located in the basin BLC and thus converge to the stable limit cycle. As a result, the basin stability estimates are S<sup>B</sup> (EP) = 0.37 and S<sup>B</sup> (LC) = 0.63, respectively. Because the separatrix, which is the unstable periodic orbit, is explicitly known for the system, the basin volumes can be determined analytically. The exact volumetric fractions of BEP and BLC in Q are 0.3275 and 0.6725, respectively. Therefore, the basin stability computed from *n* = 100 samples is a good approximation for the system at hand (Appendix A.2 indicates that *n* ≈ 300 samples are required for a very close approximation of the analytical results).

**Figure 1.** Schematic of the basin stability calculation. In the two-dimensional state space, two stable attractors EP (equilibrium position) and LC (limit cycle) co-exist. The respective basins of attraction BEP and BLC are separated by an unstable periodic orbit (indicated by the dashed line). The steady-state behaviors of *n* = 100 randomly sampled states are used to estimate the volume shares of the basins of attraction in the subset Q. The resulting basin stability estimates are S<sup>B</sup> (EP) = 0.37, S<sup>B</sup> (LC) = 0.63 for this example.
