**4. Bistability Analysis**

In all the above analysis, we did not consider the multistability in each issue to simplify the demonstration. In fact, for the special structure of symmetry, coexisting attractors exist in their own basins of attraction in phase space. Specifically, for symmetrical systems, when the symmetry is broken, a pair of symmetrical coexisting attractors usually show up.

System (2) is a rotational symmetric system, which can be proved by the invariance of transformation *x* → −*x*, *y* → −*y*, *z* → *z*, *u* → −*u*. Symmetric systems are prone to show multistability due to the effect of broken symmetry. In general, predicting multistability seems not easy in theory. A common method to identify multistability is using the basins of attraction based on the ergodic initial conditions. Alternative methods can resort to non-bifurcation manipulation, in which a linear transformation is performed on the basin of attraction to generate a dynamical dispersion for a fixed initial condition, which can reveal different coexisting symmetrical pairs by generating different average values [40].

When offset boosting is introduced from the variables *x* and *u*,

$$\begin{cases} \dot{\mathbf{x}} = \ -ay - (\mathbf{x} + d)\mathbf{z} - (\mathbf{u} - d), \\ \dot{y} = (z - c)(\mathbf{x} + d), \\ \dot{z} = -b - m(\mathbf{x} + d)y, \\ \dot{u} = k(\mathbf{x} + d) - y. \end{cases} \tag{7}$$

and the average values of variables *x* and *u* will change according to the offset control parameter *d*. The coexisting attractors are drawn into different areas of the basin since the basins of attraction of the coexisting symmetric pair of attractors also change according to the offset parameter, as shown in Figure 16. In Figure 16a, the averages of variables *x* and *u* are revised by the offset parameters, while the average values of variable *y* remains the same. From Figure 16b, the invariance of Lyapunov

exponents indicates the same structure of the symmetric pair of coexisting attractors. The typical phase trajectories of the symmetrical attractors of the system (2) under a pair of opposite initial conditions are shown in Figure 17.

**Figure 16.** Dynamical behaviors of system (7) with *a* = 5, *b* = 4, *k* = 0.5, *m* = 1 and initial condition [1, −1, −1, 1], when parameter *d* varies in [−5, 5]: (**a**) average values of the signals, (**b**) Lyapunov exponents.

**Figure 17.** Coexisting symmetrical chaotic attractors of system (2) with *a* = 5, *b* = 4, *c* = 1.3, *k* = 0.5, *m* = 1 with initial conditions IC1 = (1, −1, −1, 1) (green); IC2 = (−1, 1, −1, −1) (red).

To further verify the multistability in system (2), the basin of attraction is shown in Figure 18, which has the predicted symmetry and complex fractal structure. To show the types of attractors more clearly and comprehensively, the similar chaotic attractors are presented using an identical color in the basin of attraction. It can be clearly seen that there are two areas in different colors in the basin. The corresponding Lyapunov exponents of the two attractors are (0.2137, 0.0623, 0, -1.5761), and the Kaplan-Yorke dimension is 3.1751.

**Figure 18.** Basins of attraction of system (2) with *a* = 5, *b* = 4, *c* = 1.3, *k* = 0.5, *m* = 1 in plane of *z*(0) = −1 and *u*(0) = 0.

Both chaotic and hyperchaotic attractors show sensitivity to the initial condition, and furthermore, multistability and hyperchaos make the sensitivity more complicated. From two initial conditions in the same basin of attraction, even a slight difference results in great divergence in system (2), which is shown in Figure 19a. While from the two initial conditions in different basins of attraction, the slight difference leads to two separate phase trajectories, as shown in Figure 19b.

**Figure 19.** Dynamical behavior of system (2) with *a* = 5, *b* = 4, *k* = 0.5, *m* = 1 under different initial condition (**a**) *c* = 1; (**b**) *c* = 1.3.
