Free Response of a Rotational Nonlinear Energy Sink Coupled to a Linear Oscillator: Fractality, Riddling, and Initial-Condition Sensitivity at Intermediate Initial Displacements

Inertial coupling of rectilinear motion of a linearly sprung primary mass to a rotational nonlinear energy sink (NES), in which a second mass undergoes viscously damped rotation about a vertical axis, can lead to “targeted energy transfer” (TET) in which kinetic energy of the primary mass is transferred to the second mass, and which is then dissipated [1–3]. This approach to TET has been demonstrated in experiments [1], and in simulations of structural response [4] and of flow-induced vibration of a circular cylinder [5–8]. Beyond use in “passive” vibration suppression, a rotational NES can also be used as a “switch” to turn on or suppress flow-induced vibration of a cylinder [9]. A recent review [10] covers use of rotational NES approaches in a variety of applications.

Here, for two combinations of α and σ considered earlier [1–3,11], we investigate the part of the MIC space for which A₀(α) < η_i < A₁(α), with three main objectives.

First, we show that for (α, σ, λ) = (1, 0.01, 0), complete dissipation occurs only in a narrow range of MICs, and that most MICs lead to orbitally stable semi-trivial solutions having a single rectilinear amplitude lying just below A₀(α), but with different angular displacements. On the other hand, for (α, σ, λ) = (0.1, 0.1, 0), there is a riddled set of MICs for which complete dissipation of initial energy occurs, demonstrating that different combinations of α and σ can give rise to very different results.

Second, we show that for both combinations of α and σ, there are portions of the MIC space in which the asymptotic solution is sensitive to initial conditions, through the mechanism of transient chaos during the free response of the system. We characterize that sensitivity and show how it gives rise to both fractal and riddled basins of attraction.

Third, we consider the case in which there is direct rectilinear damping (λ > 0) and show how a measure of settling time depends on λ for both combinations of α and σ. The complexity of the basins of attraction in the λ = 0 case is shown to carry over, through transient chaos, to the case with direct rectilinear damping, and can lead to very different settling times for nearby initial conditions.

As in our previous work [11], the focus is on identifying and understanding new dynamical behavior, rather than on investigating the entire parameter space.

The remainder of the paper is organized as follows. In Sec. 2, we briefly discuss how solutions of ordinary differential equation (ODE) systems sensitive to initial conditions are also sensitive to numerical integration details. For (α, σ, λ) = (1, 0.01, 0), we present trajectories for representative MICs in Sec. 3, including one case showing transient chaos, followed in Sec. 4 by a division of the MIC space into regions in which the asymptotic behavior of trajectories originating therein is distinct. We discuss fractality and riddling in the MIC space in Sec. 5. In Sec. 6, we present results for (α, σ, λ) = (0.1, 0.1, 0) and show that there is a riddled set of MICs for which all initial energy is dissipated, whereas for (α, σ, λ) = (1, 0.01, 0), no MICs in the range A₀(α) < η_i < A₁(α) are found for which all initial energy is dissipated. We consider the effects of nonzero direct rectilinear damping in Sec. 7, followed by a probabilistic treatment of the effects of initial-condition sensitivity on settling time in Sec. 8, and a discussion and conclusions in Secs. 9 and 10.

It has long been understood that for ODE systems exhibiting sensitivity to initial conditions, there is also sensitivity to the details and accuracy of numerical integration [13,14]. This has led to development of a concept dubbed the “critical predictable time” [15], beyond which a computed solution diverges from the exact solution. This duration of fidelity depends on order of accuracy of the integrator, time-step size, and precision of arithmetic [15,16] and can also depend on algorithmic details (such as the order of sequential additions) affecting the results of finite-precision arithmetic.

Here, rather than using a high-order accurate integrator, a small time-step size, and high-precision arithmetic to precisely characterize exact solutions of Eqs. (1a) and (1b) over an extended duration, we instead fix the order of integration using an explicit (4,5) Runge–Kutta formula (implemented in matlab's ode45), control the time-step size through the relative and absolute error tolerances in the integrator, and use 64-bit arithmetic. We equate the integration tolerances, denote their common value by δ, and except as otherwise stated, set them at δ = 10⁻⁸. We do this for ensembles of initial conditions in small regions of the MIC space, so that when the solution is sensitive to initial conditions or integration details, we obtain information about the distribution of responses. This approach provides information analogous to what would be obtained by repetition of a physical experiment, where the initial condition in any realization cannot be exactly specified. When, as discussed in Secs. 3, 4, and 6, qualitatively similar results are obtained for different error tolerances, that strongly suggests that the computed behavior is ascribable to the mathematical model of the underlying physical system.

The asymptotic B_∞ (sometimes zero) is extracted by closely following the method previously described [4] using the Levenberg-Marquardt algorithm [17] with 50 iterations. Wavelet transforms were performed as previously described [11].

The dimensionless rotational damping parameter α = 1 lies to the left of all turning points of the stability boundary (see Fig. 1 of [11]), so for (α, σ, λ) = (1, 0.01, 0) orbitally stable solutions can exist for 0 ≤ B_∞ < A₀(1) = 1.18341… and A_2k−1(1) < B_∞ < A_2k(1) (k ≥ 1) (see Table 1 of [11]). Having considered 0 < η_i ≤ A₀(1) in Ref. [11], we restrict attention to 1.18341… = A₀(1) < η_i < A₁(1) = 4.2029…, which excludes harmonic solutions lying within the tongues of orbital stability, and thus restricts B_∞ to the range 0 ≤ B_∞ < A₀(1).

For η_i = 1.2 (slightly greater than A₀(1)) and θ_i = 0.1π, Figs. 1(a), 1(b), and 1(e) show that the amplitude of the rectilinear oscillation decreases asymptotically to about 1.1743, just below A₀(1). The angular displacement initially decreases to θ = 0.019π at τ = 2.78 (Fig. 1(c)), after which its amplitude grows transiently until τ ≈ 170, followed by asymptotic decay to zero. As discussed in Sec. 4.2, transient growth of the amplitude of dθ/dτ (Fig. 1(d)) is due to energy being transferred from rectilinear to rotational motion faster than it can be dissipated by the latter. Although the time series for the rectilinear and angular displacements and velocities are quite symmetric about zero, nonlinearity leads to marked asymmetry in D (defined as the RHS of Eq. (1a), and shown in Fig. 1(f)). To check sensitivity with respect to initial conditions, we compute a baseline trajectory with the nominal initial condition and integration tolerance (η_i, θ_i, δ) = (1.2, 0.1π, 10⁻⁸), and trajectories with reduced integration tolerance (1.2, 0.1π, 10⁻¹¹) and perturbed initial condition (1.2 + 10⁻⁸, 0.1π, 10⁻⁸), which we refer to as the BL, RIT, and PIC trajectories, respectively. They are graphically indistinguishable (Fig. S1 available in the Supplemental Materials).

For (η_i, θ_i) = (1.4, 0.1π), the response is shown in Figs. 2(a)–2(f). After θ rapidly decays to nearly zero, it undergoes small irregular oscillations, followed by rapid growth over 20 ≤ τ ≤ 68 (shown in more detail in Figs. S2(a)–S2(f) available in the Supplemental Materials on the ASME Digital Collection), and then decay to zero. In contrast to the (η_i, θ_i) = (1.2, 0.1π) case (Figs. 1(a)–1(f)), B varies little during the initial transient (0 ≤ τ ≤ 20). During transient growth of the oscillations in θ, both θ and η vary nearly symmetrically about zero, and D is again markedly asymmetric. The BL, RIT, and PIC trajectories are again graphically indistinguishable.

For (η_i, θ_i) = (2.5, 0.5π), Fig. 3(a) shows that θ initially decreases monotonically to −1.428π at τ ≈ 8.5, beyond which dθ/dτ > 0 until τ = 98.8 (Figs. 3(a)–3(c)), with a dominant frequency of 0.33, shown in the wavelet transform of Fig. 3(d). At τ = 98.8, the trajectory enters an interval of irregularity, during which Fig. 3(c) shows that dθ/dτ assumes both positive and negative values. This oscillatory behavior, which persists through about τ = 310, is broadband, with no dominant frequency (Fig. 3(d)), and has all of the indications of transient chaos. Beyond about τ = 310, θ undergoes regular oscillations with significant anharmonicity (Fig. 3(b); see also Figs. S3(a) and S3(b) available in the Supplemental Materials) for representative portions of the time series, dominated by frequencies 0.0800 and 0.2400 (Figs. S3(c) and S3(d), Table S1 available in the Supplemental Materials), with two higher frequencies (0.4000 and 0.5599) being less energetic and dying out more quickly as τ increases (Figs. S3(e) and S3(f), Table S1 available in the Supplemental Materials), and θ approaches 33π. During the entire evolution, the rectilinear motion (Fig. 3(e), Figs. S3(g) and S3(h) available in the Supplemental Materials) is quite regular, and similar to that shown in Figs. 2(a) and 2(b) for (η_i, θ_i) = (1.4, 0.1π).

The BL, RIT, and PIC trajectories, integrated with (η_i, θ_i, δ) = (2.5, 0.5π, 10⁻⁸), (2.5, 0.5π, 10⁻¹¹), and (2.5 + 10⁻⁸, 0.5π, 10⁻⁸), respectively, show that θ and dθ/dτ diverge near τ = 143 (Figs. S4(b) and S4(c) available in the Supplemental Materials), during the interval of transient chaos, leading to different θ_∞. On the other hand, time series of η(τ) are graphically indistinguishable (Fig. S4(a) available in the Supplemental Materials), with B_∞ (1.1712706…, 1.1712734…, and 1.1712707… for the BL, RIT, and PIC trajectories, respectively) being nearly independent of these changes in η_i and δ. The strong similarity of the angular velocity's wavelet transform (Fig. 3(d), Figs. S5(a)–S5(c) available in the Supplemental Materials) for all three trajectories during transient chaos is consistent with the qualitative similarity of the time series shown in Figs. S4(a)–S4(c) available in the Supplemental Materials. Before and after the interval of transient chaos, the wavelet transforms are nearly identical for the BL, RIT, and PIC trajectories.

To characterize irregular behavior of the angular velocity of the BL trajectory during 100 ≤ τ ≤ 300 (shown in Figs. 3(c) and S4(c) available in the Supplemental Materials), we compute the correlation dimension d_corr, following the approach of Lai and Tél [18], with details provided in Sec. A of the Supplemental Material. The slope of the logarithmic plot of the correlation integral C_N(m) versus correlation distance ɛ_corr, denoted by D₂(m) (Fig. S6(a) available in the Supplemental Materials), has a limiting value d_corr of about 3.3 as the embedding dimension m = 8, 16, 24, …, 56 increases (Fig. S6(b) available in the Supplemental Materials). This value indicates that the irregular oscillation of dθ/dτ in the BL time series is transiently chaotic. For the same time interval, the correlation dimensions of the RIT and PIC trajectories (not displayed) are 3.57 and 3.45, respectively, showing that improving the integration and perturbing the MIC do not change d_corr significantly.

For (α, σ, λ) = (1, 0.01, 0) with A₀(1) < η_i < A₁(1) and 0 < θ_i < π, Figs. 4(a) and 4(b) show B_∞ and θ_∞ for 684 MICs at nonuniformly incremented η_i and θ_i, selected from 5160 MICs (Tables S2(a) and S2(b) available in the Supplemental Materials). We see that B_∞ depends continuously and nonmonotonically on η_i and θ_i over this range. Every trajectory with 1.28 < η_i ≤ A₁(1) leads to an orbitally stable time-harmonic with amplitude near B_∞ = 1.171, as discussed in Sec. 4.3.

Figure 4(b) shows that the separation between Regions IIA (where θ_∞ = 0) and IIB (where θ_∞ = π) [11] for 0 < η_i < A₀(1) = 1.18341… extends to η_i = 1.64. Region I (where θ_∞ varies from 0 at the Region IIA boundary to π at the Region IIB boundary, with B_∞ = 0 throughout) continues to separate Regions IIA and IIB over the entire range A₀(1) < η_i ≤ 1.64, becoming exceedingly narrow as η_i approaches A₀(1) from below.

For 1.64 < η_i < A₁(1), Fig. 4(b) shows that dependence of θ_∞ on η_i and θ_i is more complex than at smaller η_i, with |θ_∞| increasing rapidly as η_i increases. We refer to this part of the MIC space as Region III (Fig. 5(a)). We next discuss behavior in the parts of Regions IIA and IIB above η_i = A₀(1) (Secs. 4.2 and 4.3), and in Region III (Secs. 4.3 and 4.4).

For several θ_i, Figs. 6(a) and 6(b) show the η_i-dependence of the asymptotic amplitude B_∞ over 1.18 ≤ η_i ≤ 1.26, in which the tendency to an η_i- and θ_i-independent asymptote near B_∞ = 1.171 is evident as η_i increases. We note the behavior for θ_i = 0.55π, where the sharp drop in B_∞ (Fig. 6(a)) is due to the “depression” in Region I and adjacent parts of Regions IIA and IIB [11], required by the fact that B_∞ varies continuously in Regions IIA and IIB to meet the value B_∞ = 0 at their boundaries with Region I, wherein every MIC leads to complete dissipation of initial energy. That the depression in B_∞ shown in Fig. 6(a) extends only a fraction of the way to B_∞ = 0 reflects the extreme narrowness of Region I for these larger η_i, and the difficulty of finding MICs in or sufficiently close to Region I [11]. Except for θ_i = 0.55 (near the depression), relatively little energy is dissipated for 1.18 ≤ η_i ≤ 1.26, and B_∞ differs from η_i by no more than 7% (Fig. 6(b)).

For θ_i = 0.1π, 0.475π, and 0.7π, Figs. 6(b) and 6(c) show how B_∞ and E_D,rot depend on η_i, with B_∞ attaining a maximum near η_i = 1.2004. (Similar results are found for other θ_i.) To understand this maximum, we consider time series of dθ/dτ for θ_i = 0.475π with η_i = 1.18, 1.19, 1.2004, and 1.21, shown in Figs. 7(a)–7(d). For η_i < 1.2004, three prominent early extrema (whose magnitudes exceed the plotting range in Figs. 7(a)–7(d); see Figs. S7(a)–S7(d) available in the Supplemental Materials) during 0 < τ < 10 are followed by a rapid drop of the amplitude of dθ/dτ, leading to a more gradual oscillatory decay to zero. Here, η_i,max = 1.2004 is the initial displacement at which the maximum B_∞ occurs. For η_i = η_i,max, the three early large-amplitude excursions of dθ/dτ are followed by small-amplitude oscillations quickly settling down to |dθ/dτ| < 0.009. For η_i > η_i,max = 1.2004 (Fig. 7(d)), the three early extrema are followed by a significant transient increase in the amplitude of dθ/dτ for τ > 10, with the instantaneous rate of energy transfer to the rotational mode exceeding the rotational dissipation rate. For such trajectories, the rotational mode essentially serves as a “buffer” or reservoir, temporarily storing kinetic energy transferred from the rectilinear motion, pending later dissipation by rotation.

As mentioned in Sec. 4.1, every trajectory with 1.28 < η_i ≤ A₁(1) = 4.2029… leads to a semi-trivial solution with B_∞ near 1.171. For θ_i = 0.1π, 0.475π, and 0.7π and 51 values of η_i in the smaller range 1.4 ≤ η_i ≤ 1.6, Tables 1 and S2(a), available in the Supplemental Materials, show the mean value, variance, and maximum variation of B_∞ for five different integration tolerances. At each of these 153 MICs, we used a second-order Shanks transformation [19] to extrapolate the B_∞ values for the five integration tolerances as δ → 0, and then calculated the mean, variance, and maximum variation of those extrapolates. Over this range of η_i, Table 1 clearly shows that (a) the mean value of B_∞ rapidly converges to approximately 1.1712734 and (b) the variance and the gap between the maximum and minimum values of B_∞ tend to zero as δ → 0. From this, we see that there is no variation in B_∞ over 1.4 ≤ η_i ≤ 1.6. We refer to solutions with B_∞,s,1 = 1.1712734… as being time-harmonic “special” solutions, with the understanding that there is a countable set of them, each having θ_∞ equal to a different integer multiple of π. (The subscript “1” refers to the fact that this is the lowest-lying, “ground-state” value, lying below A₀(1).)

For θ_i = 0.1π, 0.475π, and 0.7π, there remains a dependence of B_∞ on η_i for η_i less than about 1.35 (Table S3 available in the Supplemental Materials). For θ_i = 0.1π, ∂B_∞/∂η_i < 0 (Figs. S8(a)–S8(c) available in the Supplemental Materials), with a magnitude that rapidly decreases as η_i increases over the range 1.2 ≤ η_i ≤ 1.30, for all δ. For 1.30 ≤ η_i ≤ 1.34, |∂B_∞/∂η_i| is orders of magnitude smaller than at η_i = 1.20, decreases rapidly with decreasing δ, and is insensitive to δ if δ is sufficiently small. Beyond η_i ≈ 1.34, |∂B_∞/∂η_i| decreases by roughly one order of magnitude for each decadal decrement in δ. For θ_i = 0.1π, the fact that B_∞ assumes a value independent of η_i as δ is reduced for 1.34 ≤ η_i ≤ 1.64 clearly shows that B_∞ = B_∞,s,1 in this range.

The question naturally arises as to how a large range of MICs lead to solutions with a single rectilinear amplitude (but different values of θ_∞), given that such solutions are orbitally, but not asymptotically stable. What these MICs have in common is that each leads to one of a countable number of attracting manifolds in the part of the phase space where E > [A₀(1)]²/2, with each in turn leading to an orbitally stable solution with precisely one amplitude, namely, B_∞,s,1, where E is the energy defined in Eq. (5) of Ref. [11]. The details of these manifolds are beyond the scope of this work.

Immediately above η_i = 1.64 and extending to the bottom of the first tongue at η_i = A₁(1) = 4.2029… lies Region III, in which θ_∞ can assume values beyond the range 0 ≤ θ ≤ π. Each MIC in Region III leads to an asymptotic solution with (B_∞, θ_∞) = (B_∞,s,1, nπ), where n is an integer. For η_i < 2.1 and 43 values of θ_i (Table S2(b) available in the Supplemental Materials), Fig. 8 shows that θ_∞ is clustered near zero, with a range whose extent (maximum θ_∞ less minimum θ_∞) increases with increasing η_i. As η_i increases beyond about 2.29, the range of θ_∞ splits into positive and negative ranges, which for 2.4 ≤ η_i ≤ A₁(1) = 4.2029… maintain (a) approximately equal extents as they shift to progressively larger magnitudes and (b) approximate symmetry about θ_∞ = 0.5π.

The departure of θ_∞ from the much simpler behavior found below Region III reflects the fact that most MICs in Region III give rise to transient chaos. The increase in |θ_∞| with increasing η_i (Fig. 8) follows directly from the fact that increasing η_i leads to longer persistence of chaotic behavior, because at larger η_i, the greater difference between the initial energy corresponding to η_i and the energies corresponding to asymptotic values of B_∞ results in longer persistence of transient (rotating) chaos.

For the increments of η_i and θ_i shown in Table S2(b) available in the Supplemental Materials, the distribution of θ_∞ values in Region III suggests that we characterize the scales of MIC regions in which the asymptotic solution is invariant to initial-condition perturbations. We use the approach described in Sec. B of the Supplemental Material. For a nominal MIC (η_i,nom, θ_i,nom) and each of nine perturbation magnitudes ɛ_p = 10^−K (4 ≤ K ≤ 12), we consider (η_i,nom, θ_i,nom) and four nearby MICs (η_i,nom, θ_i,nom − ɛ_p), (η_i,nom, θ_i,nom + ɛ_p), (η_i,nom − ɛ_p, θ_i,nom), and (η_i,nom + ɛ_p, θ_i,nom)). We denote by K_min the smallest K (corresponding to the largest ɛ_p) for which all five trajectories give the same asymptotic solution (B_∞, θ_∞), and denote the local scale of this basin of attraction by $s(B∞,θ∞)=10−Kmin$⁠. For each ɛ_p such that the five trajectories do not give a single (B_∞, θ_∞), we say that the asymptotic solution is sensitive to initial conditions on the scale of ɛ_p. We have not iterated between this value and the next-largest value of ɛ_p, which would more precisely characterize the scales of MIC regions in which the asymptotic solution is invariant.

For 1.64 ≤ η_i ≤ 2.14, B_∞ = B_∞,s,1 for all MICs, so that $s(B∞,θ∞)$ actually measures the scale over which θ_∞ is the same at the central point and its four perturbed neighbors. Figure 9 shows that $s(B∞,θ∞)$ generally decreases with increasing η_i, so that θ_∞ becomes more sensitive to initial conditions as η_i increases. Again, this is because higher initial energy leads to longer persistence of transient chaos, during which the NES mass rotates, thus increasing the mean value of |θ_∞|, and decreasing the likelihood that nearby initial conditions will lead to trajectories with θ_∞ equal to the same integer multiple of π. Moreover, the fraction of the initial conditions for which $s(B∞,θ∞)<10−12$ increases with increasing η_i, suggesting that basins of attraction for each asymptotic solution (B_∞, θ_∞) = (B_∞,s,1, nπ) become riddled, as discussed quantitatively in Sec. 5. We also note that for fixed η_i between 1.7 and 1.94, $s(B∞,θ∞)$ decreases as θ_i approaches 0 or π. This can be explained by examining dθ/dτ for trajectories with MICs (η_i, θ_i) = (1.84, 0.8π) and (1.84, 0.99π) (Figs. S9(a) and S9(b) available in the Supplemental Materials, respectively), where temporally chaotic behavior persists considerably longer for θ_i = 0.99π (until roughly τ = 170) than for θ_i = 0.8π (until roughly τ = 95). This suggests that in Region III, temporally irregular behavior will persist longer as θ_i approaches 0 or π, consistent with greater sensitivity of the asymptotic state to perturbations, resulting in smaller values of $s(B∞,θ∞)$⁠.

We checked the estimate $s(B∞,θ∞)=10−12$ at (η_i,nom, θ_i,nom) = (2.09, 0.95π) by computing trajectories emanating from 49 points on a uniform grid within the square having edge length 3 × 10⁻¹² in the (η_i, θ_i) plane and centered about (η_i,nom, θ_i,nom), at which (B_∞, θ_∞) = (B_∞,s,1, π). For this nominal MIC, all 49 trajectories lead to (B_∞,s,1, π), strongly suggesting that the calculated $s(B∞,θ∞)$ results from the original four perturbed MICs ((η_i,nom − 10⁻¹², θ_i,nom), (η_i,nom + 10⁻¹², θ_i,nom), (η_i,nom, θ_i,nom − 10⁻¹²), and (η_i,nom, θ_i,nom + 10⁻¹²)) lying in a simply connected part of the region of attraction.

Comparison of Figs. 9 and S10 available in the Supplemental Materials shows that the variation of $s(B∞,θ∞)$ with η_i and θ_i is qualitatively very similar for the integration tolerances δ = 10⁻⁸ and 10⁻¹⁰. As discussed in Sec. 2, this suggests that the observed behavior is a property of the underlying ODE system and is not an artifact of limited-accuracy numerics.

For (α, σ, λ) = (1, 0.01, 0), every Region III MIC leads to an asymptotic solution (B_∞, θ_∞) = (B_∞,s,1, nπ), with a single rectilinear amplitude (Sec. 4.3). For each integer n, we refer to the resulting set of MICs as a “true” basin of attraction (Sec. C of the Supplemental Material). We also refer to “amplitude-only” basins of attraction, each being the set of MICs giving rise to a single value of B_∞, regardless of θ_∞. For (α, σ, λ) = (1, 0.01, 0), one such amplitude-only basin is the set of MICs for which B_∞ = B_∞,s,1, with θ_∞ = nπ regardless of n. We perturb each member of a random ensemble of MICs using a series of discrete magnitudes ɛ_p, as in Sec. 4.4, and define g_true(ɛ_p) as the fraction of MICs for which the asymptotic solution is sensitive to MIC perturbations on the scale ɛ_p. The slope of a logarithmic plot of g_true(ɛ_p) versus ɛ_p provides the uncertainty exponent a_true (Sec. C of the Supplemental Material). We can then use the values of a_true to identify fractal behavior (a_true not close to zero) and riddling behavior (a_true close to zero) [20].

Case 1: Portions of the MIC space where basins of attraction are neither fractal nor riddled. In Region III, we consider a rectangular neighborhood centered about (η_i,cen, θ_i,cen) = (1.68, 0.4π), a nominal MIC leading to (B_∞, θ_∞) = (B_∞,s,1, 0). In the neighborhood (1.67 ≤ η_i ≤ 1.69, 0.397π ≤ θ_i ≤ 0.403π), Fig. 10(a) shows the distribution of θ_∞ for 3000 randomly chosen points, with bands corresponding to θ_∞ = −π and π. Values of θ_∞ for a separate set of 700 randomly selected MICs in the smaller rectangle (1.684 ≤ η_i ≤ 1.687, 0.400π ≤ θ_i ≤ 0.402π) shown in Fig. 10(b) suggest that there are no fractal features between the two bands, or that the scale of such features is very small.

Case 2: Portions of the MIC space with fractal basins of attraction. We next compute the uncertainty exponent a_true in the square neighborhood of edge length 2 × 10⁻⁶ centered about (η_i,cen, θ_i,cen) = (1.8, 0.5π) using 400 randomly selected MICs in the square and four perturbed MICs at each of these for ten values of the perturbation magnitude ɛ_p in the range 10^−9.5 ≤ ɛ_p ≤ 10⁻⁷. The distribution of θ_∞ for all 16,400 MICs (Fig. 11(a)) reveals a series of bands corresponding to five values of θ_∞, qualitatively similar to fractal basins of attraction found previously in a different context [21]. Because B_∞ = B_∞,s,1 is constant over this rectangle, these represent true basins of attraction for (B_∞,s,1, nπ), each with a different n. Almost all points lie within bands corresponding to θ_∞ = 2π, 3π, 4π, and 5π. Between the θ_∞ = 3π (red) band intercepting the horizontal axis near θ_i = 0.5π − 8 × 10⁻⁷, and the θ_∞ = 2π (blue) band intercepting the vertical axis near η_i = 1.8 − 7.5 × 10⁻⁷, the fractal nature (established by the exponent a_true = 0.563 corresponding to the slope of the log [g_true(ɛ_p)] versus logɛ_p line shown in Fig. 11(b)) is more complicated, indicating incompletely resolved fractal features. For perturbations to the MICs greater than typical scales of fractality in the MIC space, Fig. 11(b) shows that the computed value of log g_true is nearly independent of log ɛ_p.

Case 3: Portions of the MIC space with riddled basins of attraction. For the square neighborhood with edge length 2 × 10⁻⁹ centered about (η_i,cen, θ_i,cen) = (2.2, 0.5π), no banding is apparent in the distribution of θ_∞ (Fig. S11 available in the Supplemental Materials) for 10,000 MICs (400 randomly selected MICs, with each serving as the center for four perturbed MICs at each of six values of ɛ_p in the range 5 × 10⁻¹³ ≤ ɛ_p ≤ 10⁻¹⁰). For no value of ɛ_p are all five values of n (one at the central MIC and four at the perturbed MICs) identical for even one of the 400 randomly selected MICs. This gives g_true = 1 for all ɛ_p, so that a_true = 0, clearly demonstrating that the basins of attraction are riddled.

For α = 0.1, nonrotating time-harmonic solutions are orbitally stable for asymptotic rectilinear amplitudes 0 ≤ B_∞ < A₀(0.1) = 0.4638… and A_2k−1(0.1) < B_∞ < A_2k(0.1) k ≥ 1 [11]. For initial rectilinear displacements 0.4638… = A₀(0.1) < η_i < A₁(0.1) = 3.7616…, the only range of orbitally stable oscillation amplitude is 0 ≤ B_∞ < 0.4638…. Here, we discuss sensitivity associated with transient chaos for (α, σ, λ) = (0.1, 0.1, 0), giving rise to fractal and riddled behavior in the MIC space. We also find a wide range of MICs leading to asymptotic solutions with B_∞ < 0.1, including zero-energy solutions of the type found in Region I for both combinations of α and σ [11], but not found for η_i > A₀(1) if (α, σ, λ) = (1, 0.01, 0).

Figures 12(a)–12(f) show η, dη/dτ, θ, dθ/dτ, B, and D for trajectories emanating from (η_i, θ_i) = (3.4, 0.3π) and (3.4, 0.5π). During 0 ≤ τ ≤ 187 (Figs. 12(c) and 12(d)), dθ/dτ oscillates about zero, apparently chaotically, for both MICs, with relatively large oscillations giving rise to nonmonotonic decay of B. For both MICs, this is followed by unidirectional rotation, in opposite directions, over 187 ≤ τ ≤ 322, during which there is no significant oscillation in the decay of B (Fig. 12(e)). After that, the trajectory emanating from (η_i, θ_i) = (3.4, 0.5π) undergoes slow decay until the initial energy is completely dissipated (with (B_∞, θ_∞) = (0, −36.58π)), while the (η_i, θ_i) = (3.4, 0.3π) trajectory leads to a semi-trivial solution with (B_∞, θ_∞) = (0.0331, 45π).

To characterize the irregularity of rotational motion shown in Fig. 12(d), we compare the BL trajectory integrated with (η_i, θ_i, δ) = (3.4, 0.5π, 10⁻⁸) to the RIT and PIC trajectories, integrated with (η_i, θ_i, δ) = (3.4, 0.5π, 10⁻¹⁰) and (3.4 + 10⁻⁸, 0.5π, 10⁻⁸), respectively. Time series of η, θ, and dθ/dτ (Figs. S12(a–c) available in the Supplemental Materials) show that the first discrepancy in dθ/dτ is apparent at τ = 66 (Fig. S12(c) available in the Supplemental Materials), beyond which the three trajectories undergo distinct irregular transients and approach different asymptotic states. The BL and RIT trajectories lead to B_∞ = 0 and the PIC trajectory to B_∞ = 0.0258, in each case with a different θ_∞. For (η_i, θ_i) = (3.4, 0.3π), the BL, RIT, and PIC trajectories lead to (B_∞, θ_∞) = (0.0331, 45π), (0, −28.38π), and (0.0327, 27π), respectively.

Figures 13(a) and 13(b) show wavelet transforms of dη/dτ for the BL trajectories emanating from (η_i, θ_i) = (3.4, 0.3π) and (3.4, 0.5π), with the most energetic frequency being about 0.16 (near 1/(2π)), in both cases consistent with the power spectra of dη/dτ shown in Figs. 13(c) and 13(d). The wavelet transforms of dθ/dτ ≡ Ω for the two MICs are shown in Figs. 13(e) and 13(f). Their much greater complexity, compared with the wavelet transforms of dη/dτ, and lack of a dominant frequency in the power spectra (Figs. 13(g) and 13(h)), is consistent with temporally irregular rotational motion (Fig. 12(d)).

The wavelet transforms and sensitivity to initial conditions strongly suggest that this irregularity, with no dominant frequency, is transient chaos. To lend support to this, we calculate correlation dimensions for the BL, RIT, and PIC trajectories, integrated with (η_i, θ_i, δ) = (3.4, 0.3π, 10⁻⁸), (3.4, 0.3π, 10⁻¹⁰), and (3.4 + 10⁻⁸, 0.3π, 10⁻⁸), respectively, over 0 ≤ τ ≤ 150, and get d_corr = 2.311, 2.235, and 2.187, respectively. For (η_i, θ_i) = (3.4, 0.5π), over 0 ≤ τ ≤ 150, we get d_corr = 2.209, 2.283, and 2.187, for the BL, RIT, and PIC trajectories, respectively. These results (d_corr > 2) and the excellent agreement among them, clearly identify the observed irregular behavior as transient chaos, and establish such behavior as a feature of the underlying ODE system.

For (α, σ, λ) = (0.1, 0.1, 0) and (1, 0.01, 0), the 0 < η_i ≤ A₀(α) portion of the MIC space can be divided into three parts [11], with B_∞ a positive and continuous function of η_i and θ_i in Regions IIA (where θ_∞ = 0) and IIB (where θ_∞ = π). These regions are separated by Region I, wherein B_∞ = 0 and θ_∞ varies continuously over 0 ≤ θ_i ≤ π.

For (α, σ, λ) = (1, 0.01, 0), we showed in Sec. 4.1 that Regions IIA and IIB (and by necessity, the extremely narrow Region I separating them) extend above η_i = A₀(1), with θ_∞ = 0 or π for η_i extending to 1.64. Beyond that value, transient chaos led to riddling and fractality in the MIC space with respect to θ_∞. For (α, σ, λ) = (0.1, 0.1, 0), however, Tables S4(a) and S4(b) available in the Supplemental Materials show that Regions IIA and IIB extend only slightly beyond A₀(0.1) = 0.4638, with values of θ_∞ and B_∞ showing sensitivity to initial conditions for η_i as low as 0.66. For (α, σ, λ) = (0.1, 0.1, 0), Figs. 14(a) and 14(b) show the distributions of B_∞ and θ_∞ for 3403 MICs (83 values of η_i over A₀(0.1) < η_i ≤ A₁(0.1), for each of 41 values of θ_i over 0 < θ_i < π) in Region III (see Fig. 5(b)), extending upward from η_i = A₀(0.1), with Tables S4(a) and S4(b) available in the Supplemental Materials showing all computed B_∞ and θ_∞ at those MICs and an additional 4067 MICs (at 166 values of η_i and 45 values of θ_i). (Extension of Regions IIA and IIB slightly beyond A₀(0.1) = 0.4638 is more readily apparent in Tables S4(a) and S4(b) available in the Supplemental Materials than in Figs. 14(a) and 14(b).) For 0.46 ≤ η_i ≤ 0.74, B_∞ is uniformly positive and depends continuously on the MICs, with only a small fraction of the initial energy being dissipated, qualitatively similar to the results for (α, σ, λ) = (1, 0.01, 0) (Fig. 4(a)). For η_i greater than about 0.66, there are MICs, e.g., (η_i, θ_i) = (0.66, 0.575π) and (0.68, 0.05π), for which over 97% of the initial energy is dissipated (Table S4(a) available in the Supplemental Materials). For η_i greater than about 0.74, the initial energy is completely dissipated for many MICs (Fig. 14(c), Table S4(a) available in the Supplemental Materials). These results have no analogue for (α, σ, λ) = (1, 0.01, 0).

For (α, σ, λ) = (0.1, 0.1, 0), Fig. 14(d) shows that beyond η_i = 0.68 the range of θ_∞ largely splits into two parts, whose extents increase as η_i increases up to about 2.16, similar to results for (α, σ, λ) = (1, 0.01, 0) (Fig. 8). This is because at larger η_i, the greater differences between the initial energy corresponding to η_i and energies corresponding to asymptotic values of B_∞ again result in longer persistence of transient chaos. For η_i > 2.16, the maximum |θ_∞| does not increase appreciably with η_i, unlike in the (α, σ, λ) = (1, 0.01, 0) case, where |θ_∞| increases nearly linearly as η_i increases. The distribution of θ_∞ at a given η_i is also quite different, with the gap between negative and positive values of θ_∞ being relatively small for η_i beyond about 2.16 for (α, σ, λ) = (0.1, 0.1, 0) and less well defined, compared with the large and clear gap for (α, σ, λ) = (1, 0.01, 0), which widens as η_i increases.

For (α, σ, λ) = (0.1, 0.1, 0) with 0 < η_i ≤ A₁(0.1), there are no true basins of attraction, as defined in Sec. 5. This is because the B_∞ = 0 solutions found for this range of η_i can have any value of θ_∞ (i.e., θ_∞ is a continuous variable and need not be an integer multiple of π), while for the semi-trivial solutions (with θ_∞ = nπ) B_∞ is a continuous variable (i.e., there are no ground-state time-harmonic special solutions to which a large number of MICs are led, as there are for (α, σ, λ) = (1, 0.01, 0)).

The fact that B_∞ = 0 solutions exist and are sensitive to initial conditions is likely to be of considerably greater interest than the value of θ_∞. We thus define the amplitude-only “zero-energy basin” as the set of all MICs for which B_∞ = 0, regardless of θ_∞.

On a rectangular grid of 34 uniformly incremented η_i in the range 0.46 ≤ η_i ≤ 3.76, and 25 θ_i (19 uniformly incremented in the range 0.05 ≤ θ_i ≤ 0.95π and three each in the ranges 0.0001 ≤ θ_i ≤ 0.01π and 0.99 ≤ θ_i ≤ 0.9999π), Fig. 15 shows the scales of the zero-energy basin, denoted by $sB∞=0$⁠, computed following the method introduced in Sec. B of the Supplemental Material, with five values of ɛ_p = 10^−2J (2 ≤ J ≤ 6). In general, $sB∞=0$ decreases with increasing η_i, consistent with evolution of the zero-energy basin from a fractal structure to one that is riddled. For the same 850 MICs integrated with δ = 10⁻¹⁰, the distribution of $sB∞=0$ (Fig. S13 available in the Supplemental Materials) is generally the same as found using the default value δ = 10⁻⁸. This lack of dependence on integration tolerance again indicates that this behavior is a property of the ODE system and is not caused by numerical artifacts.

To further characterize the zero-energy basin, we extract the uncertainty exponent, denoted by $aB∞=0$⁠, from randomly selected MICs lying in the zero-energy basin. Computation of $gB∞=0(εp)$⁠, defined as the fraction of these MICs, which when perturbed with magnitude ɛ_p, give rise to one or more trajectories with B_∞ ≠ 0, follows the approach described in Sec. C of the Supplemental Material.

Figures 16(a)–16(c) show distributions of B_∞ in three rectangular neighborhoods, each with 10,201 points (101 uniformly incremented values of both η_i and θ_i). For the neighborhoods (1.30 ≤ η_i ≤ 1.32, 0.65π − 0.01 ≤ θ_i ≤ 0.65π + 0.01) and (1.61 − 10⁻⁶ ≤ η_i ≤ 1.61 + 10⁻⁶, 0.45π − 10⁻⁶ ≤ θ_i ≤ 0.45π + 10⁻⁶), Figs. 16(a) and 16(b) reveal structure qualitatively similar to fractal basins of attraction found in other systems [21]. On the other hand, for (2.96 − 10⁻⁶ ≤ η_i ≤ 2.96 + 10⁻⁶, 0.7π − 10⁻⁶ ≤ θ_i ≤ 0.7π + 10⁻⁶), Fig. 16(c) reveals a random, fine-scaled zero-energy basin. In each neighborhood, we independently generate random MICs (i.e., not from the 10,201 points shown) until 400 lead to B_∞ = 0, with distributions in each neighborhood shown in Figs. S14(a–c) available in the Supplemental Materials. For each of at least seven ɛ_p values spanning a range differing by at least a factor of 100, we perturb each MIC to generate four additional MICs ((η_i,b − ɛ_p, θ_i,b), (η_i,b + ɛ_p, θ_i,b), (η_i,b, θ_i,b − ɛ_p) and (η_i,b, θ_i,b + ɛ_p)). The fraction of the 400 zero-energy trajectories having initial-condition sensitivity at a given scale ɛ_p is denoted by $gB∞=0$⁠. The results (Figs. 16(d)–16(f)) show that $log(gB∞=0)$ is a linear function of log ɛ_p, giving rise to $aB∞=0=0.851$ for (1.30 ≤ η_i ≤ 1.32, 0.65π − 0.01 ≤ θ_i ≤ 0.65π + 0.01), and to $aB∞=0=0.220$ for (1.61 − 10⁻⁶ ≤ η_i ≤ 1.61 + 10⁻⁶, 0.45π − 10⁻⁶ ≤ θ_i ≤ 0.45π + 10⁻⁶), corresponding to a fractal zero-energy basin in those two neighborhoods. For (2.96 − 10⁻⁶ ≤ η_i ≤ 2.96 + 10⁻⁶, 0.7π − 10⁻⁶ ≤ θ_i ≤ 0.7π + 10⁻⁶), $aB∞=0=1.10×10−3$ is nearly zero, showing that the zero-energy basin is riddled in this neighborhood.

For λ > 0, initial energy is completely dissipated for every initial condition, so we can extract a settling time (τ_0.01, the time required to dissipate 99% of the initial energy) and the fraction of initial energy dissipated by rotational motion (β) for every trajectory.

For (α, σ, λ) = (1, 0.01, λ) and θ_i = 0.5π, Figs. 17(a)–17(d) show the η_i-dependence of B_∞ for λ = 0, and of λτ_0.01, β, and θ_∞/π for λ = 10⁻⁵, 10⁻⁴, 10⁻³, 3 × 10⁻³, and 10⁻², in each case over the range 1.1 ≤ η_i ≤ 4.3.

For the well-spaced values of λ considered, Figs. 17(b) and 17(c) show that as a general trend at each η_i, λτ_0.01 increases monotonically and β decreases monotonically as λ increases. (From Fig. 17(b), it can also be seen that as a general trend at each η_i, τ_0.01 decreases monotonically as λ increases.) But the irregularity of the η_i-dependence of λτ_0.01 (and hence τ_0.01) and β leads to situations for which λ₁ > λ₂ implies neither λ₁τ_0.01(η_i, λ₁) > λ₂τ_0.01(η_i, λ₂), τ_0.01(η_i, λ₁) < τ_0.01(η_i, λ₂), nor β(η_i, λ₁) < β(η_i, λ₂). An example is the case in which τ_0.01 is greater for (α, σ, λ) = (1, 0.01, 0.01004) than for (α, σ, λ) = (1, 0.01, 0.01), even though the direct rectilinear damping parameter is larger in the former case. As a general trend, β increases as η_i increases, suggesting that more initial energy leads to rotation which is more vigorous, more sustained, or both, compared with the increase in the strength of the rectilinear motion. But this general increase of β with increasing η_i is not monotonic, and the local “dips” in β become more prominent in absolute and relative terms as λ increases. Correspondingly, τ_0.01 generally decreases as η_i increases. The irregular η_i-dependence of λτ_0.01, τ_0.01, and β is associated with the irregular η_i-dependence of θ_∞ for η_i > 1.64 (Fig. 17(d)) in the (α, σ, λ) = (1, 0.01, 0) case (with no direct rectilinear dissipation, discussed in Secs. 4.4 and 5), which is in turn rooted in transient chaos associated with sensitivity to initial conditions.

The importance of direct rectilinear damping is particularly evident for MICs just above η_i = A₀(1), where for λ = 0, B_∞ lies just below A₀(1) (Fig. 17(a)) and little rotational dissipation occurs. For these MICs and small nonzero λ, direct rectilinear damping is responsible for almost all dissipation, as is evident from the fact that λτ_0.01 varies relatively little for each small λ considered. At larger η_i, however, considerable rotational dissipation occurs during passage of the λ = 0 trajectories to final amplitudes B_∞ just below A₀(1), as reflected in the increased values of β as η_i increases (Fig. 17(c)).

For (α, σ, λ) = (0.1, 0.1, λ), with 0.46 ≤ η_i ≤ 3.76 and θ_i = 0.5π, Figs. 18(a)–18(d) show the η_i-dependence of B_∞ for λ = 0, and of λτ_0.01, β, and θ_∞ for λ = 10⁻⁵, 10⁻⁴, 10⁻³, 3 × 10⁻³, and 10⁻². For each measure, we see a division in behavior for values of η_i lying below and above about 0.91. In the initial range of η_i shown, λ = 0 values of B_∞ lie somewhat below A₀(0.1) = 0.4638…. The irregular η_i-dependence in each plot is due to transient chaos associated with sensitivity to initial conditions. For the smallest η_i shown, λτ_0.01 is nearly independent of λ, corresponding to inverse dependence of τ_0.01 on λ for small η_i. The long times τ_0.01 (ranging from about 400 for λ = 10⁻² to about 4 × 10⁵ for λ = 10⁻⁵; see Fig. S15 available in the Supplemental Materials) are associated with the fact that relatively little rotational dissipation occurs for λ = 0 before a semi-trivial solution with amplitude just below A₀(0.1) is approached. For λ > 0, dissipation in these small-η_i cases is in large part due to rectilinear damping, for which it is easily shown that energy decay is exponential in time, with exponent λ. For the smallest η_i shown, β is small, because B_∞ values for λ = 0 lie just below A₀(0.1) = 0.4638…, with little rotational dissipation occurring. As η_i increases to 0.68, β generally increases, because for λ = 0, rotational dissipation brings down the initial energy (which increases quadratically with η_i) to values below A₀(0.1). For the smaller λ values, τ_0.01 varies strongly with η_i between η_i = 0.83 and about 0.96, due to transient chaos. As η_i increases beyond about 0.91, the λ = 0 values of B_∞ decrease rapidly to below 0.1, while for λ > 0, τ_0.01 decreases to a nearly λ-independent value below 100 at about η_i = 1.01 (Fig. S15 available in the Supplemental Materials). Beyond η_i = 1.01, τ_0.01 generally increases as η_i increases (with irregularity due to transient chaos) for the λ considered, with a weak dependence on both η_i and λ. From maximum values between 0.8 and unity, β generally decreases with increasing η_i (again with irregularity due to transient chaos).

We next focus on qualitative differences between the rectilinearly damped results for the two combinations of α and σ. We begin with the nonmonotonic dependence of τ_0.01 and β on λ for α = σ = 0.1 shown in Figs. 18(b) and 18(c), compared with the monotonic dependence on λ for α = 1 and σ = 0.01 (Figs. 17(b) and 17(c)) found for the same values of λ. Figures 17(a) and 17(d) show that for (α, σ, λ) = (1, 0.01, 0) and θ_i = 0.5π, B_∞ varies smoothly (and only weakly) over the range of η_i, and that θ_∞ is constant for A₀(1) ≤ η_i ≤ 1.65, with |θ_∞| increasing monotonically over 1.65 ≤ η_i ≤ 4.2. On the other hand, for (α, σ, λ) = (0.1, 0.1, 0) and θ_i = 0.5π, Figs. 18(a)–18(d) show that transient chaos gives rise to irregular dependence of both B_∞ and θ_∞ on η_i over the entire range shown. This irregularity carries over to irregular dependence of θ_∞ on η_i for small nonzero λ, again associated with transient chaos, which in turn leads to irregular dependence of τ_0.01 and β on η_i. Comparing the distributions of θ_∞ values for λ = 0 and λ > 0 shows the importance of the λ = 0 results to cases where direct rectilinear damping occurs.

Here, we illustrate how the effects of initial-condition sensitivity on the settling time measure τ_0.01 can be treated probabilistically. The discussion is intended to highlight the statistical nature of the problem, rather than to provide a detailed discussion of an application-specific approach.

For α = σ = 0.1 with direct rectilinear damping (λ > 0), we showed in Sec. 7 that τ_0.01 displays significant initial-condition sensitivity, related to the transient chaos found in the λ = 0 case. Because initial conditions in experiments can never be known or specified precisely, the variation of settling time in the neighborhood of a nominal initial condition will be of interest in applications where a rotational NES is used to damp the free response of a linear oscillator,

For α = σ = 0.1, and 441 uniformly distributed MICs in the neighborhood (1.999 ≤ η_i ≤ 2.001, 0.5π − 0.001 ≤ θ_i ≤ 0.5π + 0.001), Figs. 19(a)–19(e) show cumulative distribution functions (CDFs) for τ_0.01 at the five nonzero values of λ considered in Sec. 7. (Figure 19(f) shows the CDF for B_∞ when λ = 0.) Table 2 shows some statistics of these distributions. The minimum and mean values τ_0.01 decrease as λ increases, as do τ_0.01(F = 0.95) and τ_0.01(F = 0.8), below which 95% and 80% of each distribution lies. The middle 90% of τ_0.01 values, τ_0.01(F = 0.95) − τ_0.01(F = 0.05), has width less than 18% of the mean for each λ. Results in Table 2 also show that, for each λ, 90% of the settling times lie within 11% of the mean. From the standpoint of probabilistic design, stronger rectilinear damping reduces mean and minimum settling times. Almost all settling times are tightly clustered, but a small number of “outlier” MICs have considerably larger τ_0.01 (e.g., for λ = 10⁻²).

For a given combination of α, σ, and λ and a given nominal initial condition, CDFs like these can be used to determine the probability that τ_0.01 will be less than a specified value. For example, for (α, σ, λ) = (0.1, 0.1, 10⁻²) and a nominal MIC of (η_i, θ_i) = (2, 0.5π), Table 2 shows that for 95% of the initial conditions, the dimensionless settling time will not exceed 116.35. For each value of the direct rectilinear damping parameter λ, τ_0.01 at the 95th percentile of the settling time distribution exceeds the value at the 80th percentile by less than four percent. This statistic, along with the CDF plots themselves, shows that, although there is a fairly broad distribution of τ_0.01, there are relatively few “outlier” MICs with very long settling times.

Beyond the results in Sec. 7 showing that τ_0.01 and its initial-condition sensitivity can depend on α, σ, λ, and the nominal MIC, it is also clear that the statistics can depend on the size of the neighborhood considered. For α = σ = 0.1 and the same nominal initial condition, (η_i, θ_i) = (2, 0.5π), Figs. 20(a)–20(e) and Table 3 show that when the neighborhood is smaller (with area 0.01% as large as the neighborhood discussed above), the CDFs are narrower. Figure 20(c) and Table 3 show that for (α, σ, λ) = (0.1, 0.1, 10⁻³), the CDF is nearly vertical (with a variance of 0.01 versus 18.49 in the larger neighborhood), corresponding to negligible variation in τ_0.01. This is a classic example of fractal behavior in the MIC space, where reducing the size of the neighborhood eventually eliminates initial-condition sensitivity as the smallest fractal scale is reached. On the other hand, for (α, σ, λ) = (0.1, 0.1, λ) with λ = 10⁻², 3 × 10⁻³, and 10⁻⁴, reducing the neighborhood area by 99.99% reduces the variances from 162.90 to 74.81, from 58.57 to 35.53 and from 62.51 to 54.50, respectively (Tables 2 and 3). These more modest reductions suggest that the smallest fractal scale (possibly with riddling between fractal bands) has not been reached. Finally, for (α, σ, λ) = (0.1, 0.1, 10⁻⁵), the distribution is much narrower for the smaller neighborhood (with variance of 1.24 compared with 67.08 for the larger neighborhood), but still considerably larger than for λ = 10⁻³, suggesting that the smaller neighborhood doesn't quite lie in a single fractal band.

Solutions exhibit significant sensitivity to initial conditions for both combinations of the dimensionless parameters considered. This is in contrast to earlier results [11] at smaller initial rectilinear displacements. (For (α, σ, λ) = (0.1, 0.1, 0), the two trajectories computed by Gendelman et al. [1] in what we call Region III were for MICs for which the response is not sensitive to initial conditions.) This sensitivity is manifested in fractality or riddling of the basin of attraction of zero-energy solutions (regardless of θ_∞) for (α, σ, λ) = (0.1, 0.1, 0). In the case where there is direct rectilinear damping, this sensitivity gives rise to sensitivity of the settling time, as seen in τ_0.01, the time required to dissipate 99% of initial energy.

Sensitivity also extends to initial conditions that are not motionless. For (α, σ, λ) = (1, 0.01, 0), Table S5 available in the Supplemental Materials shows sensitivity of θ_∞ to initial conditions in a small neighborhood of (η_i, v_i, θ_i, Ω_i) = (2.2, 0.1, 0.5π, 0), and a strong bias to positive integer multiples of π. Sensitivity persists when the initial rectilinear and angular velocities are both nonzero. For (α, σ, λ) = (1, 0.01, 0), and a small neighborhood of initial conditions near (η_i, v_i, θ_i, Ω_i) = (2.2, 0.1, 0.5π, 0.1), Table S6 available in the Supplemental Materials shows that the final angular displacements are again sensitive to initial conditions and are strongly biased to positive multiples of π. For (α, σ, λ) = (0.1, 0.1, 0), Tables S7(a) and S7(b) available in the Supplemental Materials show that the zero-energy solutions found for MICs in Region III also occur for nonmotionless initial conditions, in a neighborhood where the asymptotic solution is sensitive to initial conditions.

From the standpoint of designing a rotational NES for vibration suppression, this work has three implications. First, the fact that two combinations of the dimensionless parameters give results so qualitatively different emphasizes the importance of understanding how the response depends on those parameters. Second, one should either attempt to avoid conditions (combinations of α and σ) where behavior is fractal or riddled for expected values of the initial conditions, or use a design methodology that deals with the sensitivity, especially when the settling time is sensitive to initial conditions. One broad class of approaches to handling such sensitivity is probabilistic. This includes approaches that adapt methods of uncertainty quantification [22] to the case in which there is strong sensitivity to initial conditions [23], as well as those that sample the initial-condition space over some range and provide data to develop a probability distribution for the response (e.g., final amplitude in the case of no direct rectilinear damping, or more generally the time to dissipate a given fraction of initial energy), as discussed in Sec. 8. A second class of approaches would consider “best” and “worst” cases, where one performs computations for a sufficiently large ensemble of initial conditions, and then uses the upper and/or lower bounds of an appropriate response measure (e.g., the settling time) in design. Finally, we expect that the sensitivity reported here will carry over to systems with additional degrees of freedom, such as the dynamical models of two- and nine-story model buildings for which Saeed et al. [4] performed simulations of vibration suppression by a rotational NES.

For a rotational nonlinear energy sink coupled to a linear oscillator, we examined the part of the motionless initial-condition space (with initial rectilinear and angular velocities both zero) in which the initial rectilinear displacement ranges from the lowest amplitude at which purely rectilinear time-harmonic motion is unstable absent direct rectilinear dissipation, up to the next-largest amplitude at which such motion is orbitally stable. For one combination of the dimensionless rotational damping and coupling parameters, no motionless initial conditions lead to zero-energy solutions, while for a second combination (with the rotational damping parameter decreased by 90% and the coupling parameter increased tenfold), the basin of attraction for zero-energy solutions is riddled. For both combinations, the dissipable fraction of initial energy depends on initial conditions. For one combination, a wide range of initial conditions in a basin that can be fractal or riddled are led to time-harmonic “special” solutions with a single rectilinear amplitude lying just below the lowest branch of the stability boundary. When direct rectilinear damping occurs, sensitivity to initial conditions carries over to the settling time, which can sometimes show extreme sensitivity.

The authors gratefully acknowledge support from National Science Foundation (NSF) Grant CMMI-1363231.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper. No data, models, or code were generated or used for this paper.

m =

embedding dimension, used to calculate correlation dimension

v =

dimensionless rectilinear velocity of primary mass

B =

$η2+(dη/dτ)2$

d_corr =

correlation dimension for chaotic transients

A_k =

the kth critical value of the rectilinear amplitude of a nonrotating, time-harmonic solution, separating orbitally stable solutions from unstable solutions

C_N(m) =

correlation integral for each m

D₂(m) =

slope of the logarithmic plot of the correlation integral C_N(m) versus correlation distance ɛ_corr

α =

dimensionless rotational damping parameter

β =

fraction of energy dissipation attributable to rotational damping

δ =

absolute and relative tolerance in numerical integration

ɛ_corr =

correlation distance

η =

dimensionless rectilinear displacement

θ =

angular position of NES mass

λ =

dimensionless direct rectilinear damping parameter

σ =

dimensionless parameter coupling rotational to rectilinear motion

τ =

dimensionless time

Ω =

dimensionless angular velocity

i =

denotes an initial value

s =

denotes the amplitude of a “special” asymptotic semi-trivial solution

∞ =

denotes an asymptotic value

0.01 =

denotes the value of τ at which 99% of the initial energy has been dissipated.

BL =

baseline trajectory, computed with nominal MIC and default integration tolerances

PIC =

perturbed initial-condition trajectory, computed with perturbed value of the initial rectilinear displacement

RIT =

reduced integration tolerance trajectory, computed using a smaller integration tolerance