Non-feedback technique to directly control multistability in nonlinear oscillators by dual-frequency driving

A novel method to control multistability of nonlinear oscillators by applying dual-frequency driving is presented. The test model is the Keller–Miksis equation describing the oscillation of a bubble in a liquid. It is solved by an in-house initial-value problem solver capable to exploit the high computational resources of professional graphics cards. During the simulations, the control parameters are the two amplitudes of the acoustic driving at fixed, commensurate frequency pairs. The high-resolution bi-parametric scans in the control parameter plane show that a period-2 attractor can be continuously transformed into a period-3 one (and vice versa) by proper selection of the frequency combination and by proper tuning of the driving amplitudes. This phenomenon has opened a new way to drive the system to a desired, pre-selected attractor directly via a non-feedback control technique without the need of the annihilation of other attractors. Moreover, the residence in transient chaotic regimes can also be avoided. The results are supplemented with simulations obtained by the boundary-value problem solver AUTO, which is capable to compute periodic orbits directly regardless of their stability, and trace them as a function of a control parameter with the pseudo-arclength continuation technique.


Introduction
One of the common features of nonlinear systems is multistability [1]; that is, two or more attractors can coexist at a fixed parameter combination. Multistability appears in almost any field of science; for instance in biological systems [2,3], hydrodynamics [4], mechan-ical engineering [5], chemical reactions [6,7], neuron dynamics [8], climate dynamics [9] or social systems [10] to name a few. Multistability is also found in bubble dynamics and acoustic cavitation [11], the topic of the present study. In applications, a major challenge in nonlinear dynamics is how to control multistability and switch between coexisting states representing different system behaviour.
The mechanism of the emergence of multiple attractors can be rather different. In conservative systems, for example, one can obtain an arbitrarily large number of coexisting attractors by introducing small dissipation [12] (see also the Kolmogorov-Arnold-Moser (KAM) theory [13]). The number of attractors is roughly inversely proportional to the damping factor [14]. The appearance of homoclinic tangencies or the coupling of systems with a large number of unstable invariant manifolds can also lead to a high level of multistability. In these cases, even an infinite number of attracting periodic orbits called Gavrilov-Shilnikov-Newhouse sinks can coexist [15,16]. Another important mechanism leading to multistability is the coupling of a large number of identical systems, where the number of the attractors scales with the number of oscillators [17]. This results in a so-called attractor crowding that makes the system extremely sensitive to external noise leading to random hopping between many coexisting states. Several other mechanisms such as delayed feedback [18], parametric forcing [19] or external noise [20] can also induce multistability.
The importance of the control of multistability can be explained with the fact that different stable states represent different system performances. A good example is the large subject of chaos control [21][22][23][24][25], where the main aim is to eliminate unpredictable behaviour. Furthermore, there are various fields where the control of multistability is of paramount importance: laser technology [26], optical communication [27], cardiology [28], genetics [29] or ecology [30]. If multistability is undesired, then basically two approaches are possible: make the system monostable or preselect and stabilize an attractor against external noise. In case of a desirable multistability, where the control between different states (operations) on demand is important, the task is to properly select an attractor the system should approach and/or exclude certain undesired attractors from the dynamics. In these cases, the stabilization against external noise can also be crucial.
There are three main control strategies known in the literature, namely non-feedback, feedback and stochastic control [1]. Each has its own advantages and disadvantages. The non-feedback techniques are simple and easy to use, where the main aim is to kick the system to another stable state [31] or annihilate attractors by periodic perturbation (modulation) of a parameter or a state variable [32]. Unfortunately, it cannot be guaranteed that the system settles down onto the desired attractor. The kick introduces randomness to the control, and extremely long transients are also possible in the presence of transient chaos. In case of periodic modulation, there is no direct control over which attractors are annihilated, since usually the attractor with the smallest basin is destructed first. With increasing perturbation magnitude, the annihilation process continues with attractors having successively the next smaller basin. Applying a feedback control strategy, the desired attractor can be targeted directly; moreover, it is a reliable technique to control attractors against external noise. However, details are necessary about the state space [33] or even the solution of the attractor itself [34] or its Jacobian. This makes the feedback control technique unusable in certain problems. For instance, in an acoustically trapped bubble, it is hard to obtain state space information about the oscillating bubble and feed this information back to the resonator. The third kind of control (stochastic) also performs attractor annihilation, this time with the addition of external noise [35]. The drawback is similar as in case of the non-feedback technique: the attractors are destructed in the order of the size of their basin of attraction [1]. That is, the direct control is again lost. Attractors with very small basins are usually buried by the noise immediately, without the need of increasing the magnitude of the noise to some higher level.
The main aim of the present study is to propose a possible non-feedback control technique capable of driving a sinusoidally excited system to a pre-selected periodic attractor directly. It is based on the observation that a period-2 and a period-3 attractor can be continuously transformed into each other by adding a second, sinusoidal component to the excitation. The proper tuning of the excitation amplitudes is necessary; moreover, the fixed commensurate frequency combination must be adjusted according to the periods of the attractors being transformed. The advantages and disadvantages over other, already existing control techniques, the limitations and possible extensions to transformation of a period-n to a period-m attractor in general, and the application to other nonlinear oscillators are also discussed.
The test model is the Keller-Miksis equation well known in the bubble dynamics community. It is a second-order nonlinear ordinary differential equation describing the oscillation of a bubble in a liquid [11]. The parameter space of the present study is relatively large due to the dual-frequency driving (amplitudes and frequencies). The approach to obtain a global picture about the bifurcation structure in the fourdimensional driving parameter space within reasonable time is to exploit the high arithmetic processing power of professional graphics cards (GPUs) using an in-house initial-value problem (IVP) solver written in C++/CUDA C. The parameter scans are supplemented with results obtained by the boundary-value problem solver AUTO [36], which is suitable to compute periodic orbits directly regardless of their stability, and follow their paths with the pseudo-arclength continuation technique.

Mathematical model
The test model studied (Keller-Miksis equation [11]) describing the evolution of the bubble radius reads where R(t) is the time dependent bubble radius; c L = 1497.3 m/s and ρ L = 997.1 kg/m 3 are the sound speed and density of the liquid domain, respectively. According to the general, dual-frequency treatment, the pressure far away from the bubble consists of a static component, P ∞ , and periodic components with pressure amplitudes P A1 and P A2 , angular frequencies ω 1 and ω 2 , and with a phase shift θ . The connection between the pressures inside and outside the bubble at its interface can be written as where the total pressure inside the bubble is the sum of the partial pressures of the non-condensable gas, p G , and the vapour, p V = 3166.8 Pa. p L denotes the pressure in the liquid at the bubble-liquid interface. The surface tension is σ = 0.072 N/m, and the liquid kinematic viscosity is μ L = 8.902 −4 Pa s. The gas inside the bubble obeys a simple polytropic relationship where the polytropic exponent is γ = 1.4 (adiabatic behaviour), the equilibrium bubble radius is R E = 10 μm and the static pressure is P ∞ = 1 bar. For numerical purposes, system (1)-(4) is rewritten into a dimensionless form by introducing the dimensionless variables and the equations are rearranged to minimize the number of coefficients. The final form of the numerical model iṡ where the numerator, N KM , and the denominator, D KM , are (1 + C 9 y 2 ) − y 1 (C 7 cos(2πτ ) and The coefficients are summarized as follows: The angular frequencies ω 1 and ω 2 are normalized by the linear, undamped angular eigenfrequency [37] of the unexcited system that defines the relative frequencies as With R E = 10 µm and with the other constants in (25) given before in this section, the eigenfrequency is f 0 = ω 0 /2π = 340 kHz.

Numerical implementation and parameter choice
The simplest technique to solve the system (8)-(9) is to take an initial-value problem solver, integrate the equa-tions forward in time and, after the convergence of the transient trajectory, record and save the important properties of the solutions found. Thereby, our main interest is in periodic solutions and their organization in the reduced driving parameter space discussed below. The first part (2048 T ) of the integration is regarded as initial transient and dropped followed by integration of the system further by 8192 T for the proper convergence of average quantities like Lyapunov exponent or winding number, where T is the period of the dualfrequency excitation function. The scheme employed is the explicit and adaptive Runge-Kutta-Cash-Karp method with embedded error estimation of orders 4 and 5 (the algorithm is adapted from [38]). The quantities stored are the points of the (global) Poincaré section (see the choice in Sect. 4 below), which are standard ingredients for bifurcation diagrams; furthermore, the maximum bubble radius and the maximum absolute value of the bubble wall velocity, which are important for applications; finally, the period, the Lyapunov exponent and the winding number of the attractors found, quantities that are essential for a detailed analysis of bifurcation structures. A strategy to represent the results of parametric studies involving high-dimensional parameter spaces consists in creating high-resolution bi-parametric plots, a rapidly spreading technique in the investigation of nonlinear systems with a high-dimensional parameter space [39][40][41][42][43][44][45][46][47][48][49]. The system studied here, a bubble in water with dual-frequency acoustic excitation, has a four-dimensional driving parameter space (P A1 , P A2 , ω R1 , ω R2 ). For simplicity, the phase shift is set to θ = 0. It is restricted in the present study by the requirement that the relative-frequency pairs obtained from ω R1 and ω R2 are rational. The pressure amplitudes P A1 and P A2 are taken as the main control parameters at fixed relative-frequency pairs. The range of each pressure amplitude is 0-5 bar, investigated at first with a resolution of 501 steps for an overview. The relativefrequency values ω R1 and ω R2 are chosen from the following set: Exploiting symmetries in the driving parameter space, this gives 8 i=1 i = 36 combinations of two frequencies. Observe that two orders of magnitude difference in the frequency range are covered and that, for instance, ω R1 = 1/5 and ω R2 = 1/1 means that in an experiment, the bubble is driven by 0.2 ω 0 and by the resonance frequency ω 0 . Computations are performed at every possible frequency combination of the above set, which means altogether 36 two-dimensional plots for every quantity to be studied (period, Lyapunov exponent, etc). Moreover, owing to the fine details in the diagrams, high-resolution surveys are necessary. In order to reveal possible coexisting attractors, 10 randomized initial conditions are used at each parameter set. Indeed, the 10 different trajectories can converge to different attractors with distinct properties (e.g. period). Altogether, a single plot contains 501 × 501 × 10 ≈ 2.5 million initial-value problems.
Even a single bi-parametric plot requires large computational capacities; thus, the high processing power of professional videocards (GPGPUs) is exploited. During the parameter studies, millions of simple (only second order) and independent ordinary differential equations are to be solved. Moreover, the algorithm employed needs only function evaluations. Therefore, our problem is well suited for parallelization on GPUs. The program code is written in a C++ and CUDA C software environment. The available GPUs are a Titan Black card (Kepler architecture, 1707 GFLOPS double-precision processing power), two Tesla K20 cards (Kepler architecture, 1175 GFLOPS) and eight Tesla M2050 cards (Fermi architecture, 515 GFLOPS). The application of doubleprecision floating-point arithmetic is mandatory in bubble dynamics due to the possibly large-amplitude, collapse-like response of the system [11]. The final, optimized code is approximately about 50 times faster on the Titan Black, about 30 times faster on the Tesla M2050 and about 120 times faster on the Tesla K20 than on a four-core Intel Core i7-4790 CPU. Surprisingly, the professional Tesla K20 card was more than two times faster than the gamer Titan Black card even though the theoretical peak processing power is higher for the Titan Black GPU. A thorough performance analysis between GPUs and CPUs solving large numbers of ordinary differential equations is beyond the scope of the present study. The interested reader is referred to the publications [38,[50][51][52] for more details.

Global Poincaré section and the period-reducing phenomenon
The parameters investigated in the present study are the pressure amplitudes (P A1 , P A2 ) and the relative frequencies (ω R1 , ω R2 ) of the dual-frequency excitation with phase shift θ = 0, i.e. C 12 = 0. The main control parameters are the pressure amplitudes at several fixed frequency combinations. The ratio of each frequency pair is rational; therefore, periodicity of the dual-frequency driving is guaranteed (quasiperiodic forcing is excluded).
As an example, a dual-frequency excitation is presented in Fig. 1 as a function of time together with its two periodic components at ω R1 = 3 and ω R2 = 2. For definiteness, the pressure amplitudes are set to P A1 = P A2 = 1 bar. The black and blue curves are the highand low-frequency components, respectively. Since the dimensionless time τ is defined by means of the first frequency component, see Eq. (5), the normalized period corresponding to the relative frequency ω R1 is T 1 = 1, indicated by the vertical black line in Fig. 1. The normalized period of the second component (relative frequency ω R2 ) can be calculated with the frequency ratio; that is, T 2 = 1/C 11 = ω R1 /ω R2 = 3/2 = 1.5 (blue vertical line). Compare these periods with the arguments of the harmonic functions in equation (10) and keep in mind that the phase shift is zero (θ = C 12 = 0). As only normalized periods will be used throughout the article, this very often occurring quantity will be simply called "period".
The sum of the two harmonic components results in a still periodic, but non-harmonic forcing function presented by the red curve in Fig. 1. The period of this function is T = 3T 1 = 2T 2 = 3, which is used as a global Poincaré section for the system (8)-(9). Consequently, during the numerical integration of the system with an initial-value problem solver, the continuous trajectories are sampled at time instants τ n = n · T (n = 0, 1, 2, . . .). Accordingly, the period, T N , of a periodic orbit is defined as T N = N · T, N ∈ IN, and the solution is called period-N orbit. Observe, however, that in the special cases of P A1 = 0 bar or P A2 = 0 bar, the period T , which defines the global Poincaré section, does not coincide with the period of the actual driving T 2 or T 1 . For P A1 = 0 bar, the dual-frequency forcing becomes a purely harmonic function shown by the blue curve in Fig. 1, and the continuous trajectories are sampled only after every second period of excitation, 2T 2 . Therefore, an originally period-2k solution is observed only as a period-k orbit. Similarly, when P A2 = 0 bar, the trajectories are sampled only after every third period of the excitation, 3T 1 (compare the black and red vertical lines in Fig. 1), and the period of an originally period-3k orbit is reduced from 3k to k. This period-reducing phenomenon, occurring when one of the two components of the forcing function is vanishing and only a purely harmonic excitation is left, plays an important role in the bifurcation structure of the system, discussed in more detail in Sect. 5; furthermore, period reduction necessarily takes place for other frequency pairs as long as the ratio of a pair is rational.

Results from a global scan
Throughout the following subsections, the fundamental mechanism capable of controlling multistability in harmonically driven nonlinear oscillators is presented. The basis is the period-reducing phenomenon discussed above, by which the dual-frequency treatment can generate "bridges" between different periodic attractors. First, a global picture is shown via a series of highresolution bi-parametric plots of Lyapunov-exponent and periodicity diagrams at different frequency combinations (the control parameters are the excitation amplitudes). Next, the details are demonstrated through a specific example with frequencies ω R1 = 3 and ω R2 = 2 (see again Fig. 1 for the driving signals). It must be emphasized that the amplitudes for both frequencies of the driving can be arbitrary; that is, there is no restriction to small-amplitude perturbations corresponding to the second sinusoidal component.

The branching phenomenon
Out of the data sets for the 36 relative-frequency combinations, four diagram pairs are shown in Fig. 2 (3,1) and (3,2) organized in rows and for two different quantities (Lyapunov exponent and period) organized in columns. The left column contains the Lyapunov-exponent diagrams, where the colour-coded area means chaotic solutions (positive exponent) and where in the greyscale domains there are periodic attractors (negative exponent). In the right column, the periods of the (converged) solutions are presented up to period-6. Periodic solutions higher than six including chaos can be found in the black domains.
In the upper two rows, the plots correspond to relative frequencies lower than or equal to the linear resonance frequency (ω R1,2 ≤ 1), while the lower two rows have relative frequencies higher than or equal to the linear resonance frequency (ω R1,2 ≥ 1).
Since detailed investigations of dual-frequency driven systems for arbitrary driving amplitudes with rational frequency ratio are absent in the literature, only very few information exists about the bifurcation structure shown in Fig. 2. It can be clearly seen that the dual-frequency driving causes a very complex interplay between the resonances originating from the singlefrequency driven system (horizontal and vertical axes in Fig. 2), see especially the first two rows of the figure. Their comprehensive topological analysis is beyond the scope of the present study; instead, it intends to add some pieces to our knowledge by looking for special features in the diagrams. In particular, we focus on a very specific phenomenon, the branching mechanism, turning out to be a special feature of dual-frequency driven nonlinear oscillators with rational frequency combinations, which can be used for the control of multistability.
Observe that on the horizontal axis in the third row and second column of Fig. 2, three period-3 branches emerge from the segment marked by P3B3 (period: 3, number of the branches: 3). Similarly, three branches appear from the segments P4B3 and P3B3 on the horizontal axis, and two branches are merged together at segments P2B2 and P3B2 on the vertical axis in the last row and second column of Fig. 2, better to be seen Fig. 2 List of high-resolution bi-parametric plots as a function of the pressure amplitudes at different relative-frequency pairs organized in rows. Left column: Lyapunovexponent diagrams separating the chaotic and periodic solutions. Right column: Periodicity diagrams, where periodic solutions are plotted with different colours up to period-6. Period-7 or higher solutions (including chaos) can be found in the black regions. The pressure amplitudes P A1 and P A2 are given in bar in the enlargement and extension in Fig. 3. Observe also that the number of the branches is exactly the same as the value of the corresponding relative frequency along the respective axis. This phenomenon will be referred to as branching mechanism throughout the paper and will be investigated in more detail for the case ω R1 = 3, ω R2 = 2.
In order to get a deeper insight into the branching mechanism at the relative-frequency pair ω R1 = 3, ω R2 = 2, computations were repeated with much higher resolution of the pressure amplitudes (Fig. 3). The ranges of the control parameters are P A1 ∈ [0, 8] (bar) and P A2 ∈ [0, 5] (bar) with resolutions 4001 and 2501, respectively. The number of initial conditions at a control parameter pair (P A1 , P A2 ) is increased from 10 to 5 × 5 = 25 defined on an equidistant grid on the y 1 -y 2 phase plane. Thus, 2501×4001× 25 ≈ 250 million initial-value problems have to be solved. This huge amount of computations is divided into 5 × 8 = 40 (ΔP A1,2 = 1 bar) segments and distributed to different GPUs.
At first sight, Fig. 3 contains only slightly more information about the branching mechanism than the bottom right diagram in Fig. 2, since it is hard to represent the many coexisting attractors in a single plot obtained by the increased number of initial conditions (in case of coexistence, the solution with the highest period is plotted). Investigating the system periodby-period, however, this high-resolution bi-parametric scan becomes very helpful to understand the branching mechanism (explored rigorously in the next sections); in particular, to identify the connected branches with different periods confidently, such as the ones originating from P4B3, P3B3, P2B3, P3B2 or P2B2, where, for instance, P2B3 and P2B2 are connected.
Typical colour combinations can be seen at the border of the branch families. The three green P2 bands are bordered each by a yellow (period-4) band at the right border (a period-doubled zone) followed by a thin black line (indicating period-8 and higher periods including chaos that are encoded that way). Similarly, the three dark blue P3 bands are bordered each by a light blue (period-6) band, again indicating the start of a perioddoubling sequence. In addition, the three yellow P4 bands are bordered by a thin black line, presumably also the beginning of a period-doubling sequence to be seen only at higher resolution. Note that also other dark blue areas are bordered on one side by a period-doubled zone of light blue followed by black. Thus, the branches themselves are period-doubling objects. It should also be noted that in any case, two of the three branches in each Pn B3, n = 2, 3, 4, turn to the left (lower driving amplitudes P A1 ) and one to the right (higher driving amplitudes P A1 ).

Coexisting period-1 solutions
Although the red period-1 domain in Fig. 3 shows no sign of multiple branches, the examination of its inner structure-the organization of the coexisting period-1 orbits-reveals important features of the branching mechanism. They turn out to be general for other periods as well. By filtering the results for period-1 orbits and by counting the number of the different attractors at a given control parameter pair, the structure of the coexisting period-1 solutions can be made visible in the P A1 -P A2 plane, see Fig. 4 with a magnification of the zone near the origin in the bottom panel. Domains with different colours represent areas with a different number of coexisting attractors up to four. These domains are bounded by saddle-node (SN, solid) and period-doubling (PD, dashed) curves computed by the boundary-value problem solver AUTO introduced in more detail later. The boundaries computed by AUTO are continuous curves. Some fractal-like shapes in the colour-coded domains are an artefact of the lack of a sufficient number of initial conditions or the presence of extremely long transients during the initial-valueproblem computation Crossing one of the boundaries means the increment or decrement of the number of the attractors by one.
Multiple stable period-1 orbits in a bi-parametric space usually originate via cusp bifurcations forming pairs of SN points (hysteresis) and giving rise to hysteresis zones. The overlapping of such hysteresis zones may produce domains with a high number of coexisting period-1 attractors [53]. The sole cusp bifurcation in the bottom panel of Fig. 4, however, cannot explain the mosaic-type pattern of the period-1 structure. The investigation of the two special cases, where one of the amplitudes of the dual-frequency driving P A1 or P A2 is zero (vertical or horizontal axis in Fig. 4), can help to understand this pattern.

5.3
The limiting case P A1 = 0 and the period-reducing phenomenon The bifurcation structure with P A1 = 0 bar is shown in Fig. 5, where the first component of the Poincaré section Π(y 1 ) is plotted as a function of P A2 . The colour code is the same as in the cases of Fig. 3 and the right column of Fig. 2. For instance, red, green and blue curves are period-1, -2 and -3 orbits, respectively. Some of the periodic orbits are also marked by N ×PX, where N is the number of the coexisting attractors of period X . The pitchfork, saddle-node and period-doubling bifur- The period-reducing phenomenon discussed in detail in Sect. 4 plays an important role in the period-1 bifurcation structure. Observe that the applied frequency combination is the same as in Fig. 1; furthermore, the first component of the dual-frequency driving is zero. Consequently, the driving is a purely harmonic function proportional to the blue curve in Fig. 1, and an originally period-2k solution is observed only as period-k orbit. Thus, the first bifurcation of the stable period-1 curve originating from the equilibrium solution y 1E (both amplitudes are zero) is a pitchfork point PF, instead of a period-doubling point PD, which gives birth of two period-1 branches 2 × P1. Only the subsequent bifurcations exhibit a Feigenbaum cascade (2 × P2 → 2 × P4 → · · · ). Similarly, the first bifurcation point in a cascade must be a symmetry-breaking bifurcation in case of a symmetric orbit in a symmetric nonlinear oscillator [54]. Two examples of bubble radius versus time curves y 1 (τ ) are plotted in the top panel of Fig. 6 so as to help understanding the period-reducing mechanism. The blue and black period-1 solutions are initiated at P A2 = 1 bar from the black dots marked by the bold numbers 1 and 2 in Fig. 5, respectively. The blue harmonic function in the bottom panel in Fig. 6 is the purely harmonic driving of the system, see also Fig. 1. It is clear that if the Poincaré sections are taken at time instances τ n = n · T (n = 0, 1, 2, . . .), which is used during the present simulations, the solutions are regarded as period-1 orbits. In contrast, if the Poincaré sections are taken at time instances τ n = n ·T 2 (n = 0, 1, 2, . . .), which is a common choice in case of a purely harmonic driving, the solutions are period-2 orbits. To prevent ambiguity in the definition of the periodicities, solutions will be referred to as "originally" period-N orbits if the latter definition of the Poincaré section is applied. Although the blue and black attractors in Fig. 6 are different (each has its own basin of attraction), the only difference between them is a phase shift in time by one period of the driving (Δτ = 1.5).
Altogether two types of bifurcation scenarios of periodic orbits are observable in Fig. 5. For odd basic periods (period before a bifurcation cascade takes place), the first point is a pitchfork bifurcation followed by a Feigenbaum period-doubling cascade (PF →   PD → PD → · · · ). Some examples are the already investigated period-1 curve P1 starting from y 1E , the period-3 orbits P3 originating from a SN point at P A2 = 1.5 bar or the period-5 branches P5 appearing again via a SN point at P A2 = 2 bar. Observe that the PF points have no influence to the colour, i.e. period. In case of even basic periods, the period is immediately halved (period-reducing phenomenon), and the bifurcation scenario is a classic Feigenbaum cascade without an initial PF point (PD → PD → · · · ). An example here is the originally period-6 solution 2 × P3 generated at P A2 = 2.5 bar. Through the discussion of Figs. 5 and 6, it has been shown that the definition of the Poincaré section has an important influence on the apparent periodicity of the solutions. However, problems only appear in the two limiting cases of either P A1 = 0 bar or P A2 = 0 bar. In the first special case of P A1 = 0 bar, a suitable choice can be either T 2 or T , see the bottom panel of Fig. 6. Since the addition of any small value to the first pressure amplitude P A1 will destroy the special structure introduced in Fig. 5, T 2 will be no longer a period of the dual-frequency driving but only T . Therefore, in order to keep the definition unified, the usage of T as a global Poincaré section is advantageous. This, at first sight only technical, extension of the definition to the limiting, purely harmonic-driving case is of great help in ordering and explaining the bifurcation scenarios of dual-frequency driven systems.
In order to illustrate this effect, it is useful to solve and compute complete bifurcation curves of periodic orbits with a boundary-value-problem solver. As an example, in Fig. 7A, the second component of the Poincaré section Π(y 2 ) of the period-1 curves is plotted as a function of the control parameter P A2 computed by the bifurcation analysis and continuation software AUTO [36]. Here, the first component of the pressure amplitude P A1 is still zero. AUTO is capable of tracking down whole bifurcation curves including the unstable part even if they contain multiple turning points, and it can detect the bifurcations and their types. This is the reason why AUTO is commonly used to study the bifurcation structure of nonlinear systems [5,53,[55][56][57][58][59][60][61]. In Fig. 7, the stable and unstable parts of the period-1 branches are indicated by the black and red curves, respectively. The curve originating from y 2E becomes unstable via a pitchfork bifurcation PF at P A2 = 0.19 bar producing two stable coexisting period-1 branches. These stable curves lose their stability approximately at P A2 = 3 bar by a period-doubling point PD.   Fig. 7B. Likewise, the curves P1 E,u 1 and P1 2 2 presented in Fig. 7A form again a single curve Fig. 7B composed by a stable and an unstable part separated by a saddle-node bifurcation (black dot) at the turning point. Elevating P A1 from 0.3 bar to 0.6 bar, a pair of SN points appears along the branch P1 E 1 + P1 1 2 forming a hysteresis. Observe that both curves P1 E 1 + P1 1 2 and P1 E,u 1 + P1 2 2 have their own separate evolution as the amplitude P A1 changes; consequently, the black and blue curves in the top of Fig. 6 must become different (not only by a phase shift in time) for P A1 > 0 bar as they have their own separate evolution as well. All the pitchfork structures in the bifurcation scenarios initiated from odd basic periods are broken up in the same way as presented in Fig. 7B, except for the formation of the hysteresis which not necessarily happens.

5.5
The limiting case P A2 = 0 Shortly in this section, the complex development of the period-1 structure in a wide range of the amplitudes P A1 and P A2 will be discussed in detail. Up to now, it is important only to understand how the period-reducing mechanism works for single-frequency driving, how it decomposes periodic solutions into multiple orbits of lower periodicities, and how these multiple orbits are separated under the perturbation of the amplitude of the other harmonic component. In order to understand the global picture, however, the investigation of the other special case, where the control parameter is P A1 and the amplitude of the second component P A2 is zero, is necessary. The corresponding bifurcation diagram is shown in Fig. 8, where the first component of the Poincaré section Π(y 1 ) is plotted as a function of the pressure amplitude P A1 . The colour code and the notation system are the same as in Fig. 5. At first sight-in view of (2) with θ = 0-both cases (P A1 = 0 and P A1 = 0) seem exchangeable. However, switching to the other frequency component, different bifurcation scenarios develop (if the frequencies are distinct).
There are two types of bifurcation scenarios of periodic orbits. If the basic (original) period is divisible by 3, then the period reducing takes place immediately decreasing the period of the orbits from 3k to k followed by a Feigenbaum cascade as usual (PD → PD → · · · ). Notable examples are the branches 3 × P1 and 3 × P3 both appearing via SN bifurcations at P A1 = 1.2 bar and at P A1 = 4.1 bar, respectively. For other basic periods, there is no period-reducing phenomenon, but only a period-doubling cascade. Examples of this scenario are the period-1 curve departing from the equilibrium solution y 1E , the period-4 orbits generated at   Fig. 4 by the yellow domain. The period-1 orbit P1 emerges from the equilibrium point y 1E and remains stable up to the pressure amplitude P A1 ≈ 7 bar. It loses its stability there via a period-doubling bifurcation. The originally period-3 family of solutions 3 × P1 appearing through a SN point is decoupled into three period-1 branches existing between the amplitudes P A1 ≈ 1.19 bar and P A1 ≈ 4.67 bar. Therefore, there is a wide range of amplitudes P A1 where four period-1 attractors exist. Figure 9 shows these coexisting orbits at P A1 = 2 bar starting from the black dots in Fig. 8. The dashed curve is the period-1 orbit related to the curve originating from y 1E . The black, red and blue solid curves are the orbits of the period-reduced 3 × P1 branches, where the only difference between them is again a phase shift in time, in this case of one or two periods T 1 (Δτ = 1 or 2).

Perturbation of the limiting case P A2 = 0
Again, the boundary-value-problem solver AUTO can help to understand the behaviour of the system under a small perturbation of the amplitude of the second, harmonic component of the driving (P A2 > 0). The bifurcation curves related to the four period-1 orbits can be seen in Fig. 10A for P A2 = 0 bar. The black and red curves mean stable and unstable orbits, respectively. The bifurcation of the saddle-node (SN) and period-doubling (PD) points are also detected. They are marked by dots (SN) and crosses (PD). Observe that the P A1 -ranges of the stable (black) segments in Fig. 10A coincide with the domain of existence of the corresponding parts of the red curves in Fig. 8. On the right hand side of Fig. 10, only the evolution of the curves 3 × P1 is presented by increasing the amplitude P A2 from zero to 0.1 bar. It is clear that while the saddle-node bifurcations (black dots) take place at the same value of P A1 when P A2 = 0 bar (see the vertical thin line in Fig. 10B), the locations of these bifurcation points differ when P A2 > 0. This means that each of the three bifurcation curves have their own evolution as the amplitude P A2 changes. Consequently, the black, red and blue solid curves in Fig. 9 become different (not only by a phase shift in time) for P A2 > 0.

Period 1: the complete picture
Due to the exhaustive discussion in the preceding sections, it is clear now how the system behaves under single-frequency driving when a second, harmonic component is added as a perturbation. The huge amount of results presented in Fig. 3 can help to establish a global overview of the evolution of all the period-1 solutions by means of filtering the results by period, here period-1. Figure 11 shows these period-1 solutions represented in a three-dimensional plot, where the second component of the Poincaré section, Π(y 2 ), the bubble wall velocity, is plotted as a function of the two main control parameters P A1 and P A2 . The different period-1 orbits can be decomposed into several layers indicated by the blue, black, red and green surfaces. In the limit cases when one of the pressure amplitudes is zero (planes Π(y 2 )-P A1 or Π(y 2 )-P A2 ), the points are plotted with bigger markers to be able to distinguish such special solutions more easily. These solutions, which form complete bifurcations curves presented in Figs. 7A and 10A, are marked by P1 m n , where n is the original period and m = 1, . . . , n is a serial number, P1 means period-1 solution. When both P A1 and P A2 are zero then the system is in equilibrium and the dimensionless bubble wall velocity and its Poincaré section y 2E are zero (it is the origin of the diagram). The black period-1 surface originating from this point is marked by P1 E 1 (instead of P1 1 1 , emphasizing the equilibrium origin E). There is only one curve of this type, which is a surface of small-amplitude bubble oscillations. Fig. 11 3D representation of the stable period-1 orbits (Poincaré section points of the bubble wall velocity) above the parameter plane P A1 -P A2 by filtering the results presented in Fig. 3 by period, here period-1. The period-1 orbits can be separated into four different layers marked by the blue, black, red and green surfaces. The indication of the curves in the Π(y 2 )-P A1 (Fig. 10) or Π(y 2 )-P A2 (Fig. 7) planes is P1 m n , where n is the original period and m = 1, . . . , n is a serial number, P1 means period-1 solution. The black surface originating from the equilibrium solution of the system is marked by P1 E 1 (instead of P1 1 1 to show the connection with the equilibrium point). (Color figure online) Figure 11 indicates a very interesting phenomenon, namely the originally period-2 solutions on P1 1 2 can be transformed into the originally period-3 solutions on P1 1 3 through the blue surface. Similarly, orbits on P1 2 2 can be transformed into orbits on P1 2 3 via the green surface. These are smooth transformations by means of the continuous change of the amplitudes P A1 and P A2 of the dual-frequency driving. Naturally, all these orbits are considered as period-1 attractors due to the specific choice of the global Poincaré section. Consider, however, the following scenario of parameter variations visualized by the supplementary movie file (Online Resource 1) where the top panel shows the dimensionless bubble radius curves y 1 as a function of the dimensionless time τ , and the bottom panel represents the dual-frequency driving signal P A (τ ) (excluding the static ambient pressure).
Firstly, let us start the investigation with a solution lying on the bifurcation curve 2 × P1 in Fig. 5 (P1 1 2 or P1 2 2 in Fig. 11) so that the amplitude of the single-frequency excitation is somewhere between P A2 = 1.4 bar and P A2 = 2.9 bar (P A1 = 0 bar) with relative frequency ω R2 = 2. Specifically, in Online Resource 1, the pressure amplitude is chosen to be P A2 = 2.5 bar. In this video, besides the originally period-2 (red) attractor, which is the initial state of the transformation, an originally period-3 (blue) attractor coexists regarded as a final target of the control problem. The colour code of these attractors is the same both in Fig. 5 and in the top panel of Online Resource 1; in addition, by comparing the top and bottom panels, the periodicities of the attractors are clearly visible.
Secondly, guide the period-2 orbit smoothly to the curve 3 × P1 shown in Fig. 8 (either P1 1 3 or P1 2 3 in Fig. 11) by varying the pressure amplitudes continuously. By the end of the second operation, the system is again driven by a single frequency (ω R1 = 3) with pressure amplitude between P A1 = 1.2 bar and P A1 = 4.7 bar (P A2 = 0 bar). In Online Resource 1, this transformation takes place along a circle in the P A1 -P A2 plane with radius 2.5 bar; that is, the final value of the first pressure amplitude is P A1 = 2.5 bar. The continuously changing trajectory during the transformation is marked by the black curve in the top panel.
Thirdly, keep the single-frequency driving but change both the amplitude P A1 and the relative frequency ω R1 of the driving (in general) to transform the solution from the red curve 3 × P1 in Fig. 8 to the blue curve P3 in Fig. 5. This third step is possible, since both the 3 × P1 and the P3 curves are related to the same subharmonic resonance of order 1/3, which forms a con-tinuous domain in the pressure amplitude-frequency parameter plane [53,61]. Therefore, a smooth transformation of a period-2 solution into a period-3 orbit of a single-frequency driven system is possible by a temporary addition of a second harmonic component to the driving. This is an efficient way to control multistability in harmonically driven nonlinear oscillators. According to the best of the authors' knowledge, this phenomenon has not yet been published in the literature. During this third step in Online Resource 1, the change of P A1 is not necessary and kept constant. Observe how the initial and the final parameter set become identical; that is, single-frequency driving with amplitude P A1,2 = 2.5 bar at ω R1,2 = 2. The only difference is the state of the system: operation on the period-3 (blue) instead on the period-2 (red) attractor.
In order for a better visualization of the surfaces presented in Fig. 11, hundreds of bifurcation curves are computed by AUTO departing from the results shown in Fig. 10A. These surfaces, composed by lines, are summarized in Figs. 12, 13 and 14. The stable and unstable parts are the black and red curves, respectively. As a skeleton, all the bifurcation curves of the limit cases (P A1 = 0 bar or P A2 = 0 bar) are presented by the blue (stable) and the green (unstable) curves, which can help to identify the surfaces more easily (compare teresis, see also the bottom panel of Fig. 4. An example for this hysteresis, which exists between the regions of P A1 > 0.33 bar and P A1 < 1.19 bar, is already introduced in Fig. 7B for P A1 = 0.6 bar. As expected, the red, standalone surface in Fig. 11 has no connections with other surfaces, see Fig. 13. The green surface in Fig. 11 together with its unstable counterpart forms a simple banded, dual-layer surface presented in Fig. 14.

Higher periods
If the smooth transformation works between certain period-2 and period-3 orbits, there must be other combinations of periodic solutions where the transformation is possible. In order to investigate these other possibilities, the number of the coexisting period-3 and period-4 attractors is plotted in the top and the bottom panels of Fig. 15, respectively. The colour codes are the same as in case of Fig. 4 (period-1 orbits). Let us focus only on those parts of the complex bifurcation structures where the transformations take place. In the top panel of Fig. 15, there are two branches along which the following conversion occurs: 2 × P3 = P6 ↔ 3 × P3 = P9; that is, an originally period-6 solution is transformed into an originally period-9 orbit. Similarly, the bottom panel of Fig. 15 shows four bands converting originally period-8 orbits into period-12 orbits (2 × P4 = P8 ↔ 3 × P4 = P12). As a generalization, the present study found that applying the relativefrequency combination ω R1 = 3 and ω R2 = 2, a transformation between period-2k and period-3k attractors is possible, here k = 1, 2, . . ..

Winding numbers and compatibility of periodic orbits
Naturally, not every kind of period-2k and period-3k orbits can be transformed into each other; in other words, there is a "compatibility" condition. Observe that the first PD points of the period-1 branches in Fig. 5 (2 × P1) and in Fig. 8 (3 × P1) are connected via the period-doubling curve PD shown in the top panel of Fig. 4. It is well known in the literature that the quantity called winding number w = nk/mk is invariant and rational near a period-doubling or a saddle-node point [62,63]. Here, mk is the period of a given attractor and nk is the number of the rotations of a nearby trajectory around the attractor during one period. Therefore, to each bi-parametric PD or SN bifurcation curve a single, rational winding number can be associated. When a transformation between period-2k and period-3k orbits is possible, then it must be possible between their corresponding bifurcation points as well, see, for instance, the PD curve in Fig. 4 which has the winding number w = 1/2. Consequently, the transformation is possible only if the winding numbers of such bifurcation points are compatible; that is, when they are equal.

Discussion
According to the literature, the only option to directly target an attractor in case of multistability is to apply a feedback control technique. The main drawback of the feedback control methods is that explicit knowledge of the state space (modified targeting method [33], bushlike paths to a pre-selected attractor [64], reinforcement learning [65]) or the properties of the targeted attractor (trajectory selection by a periodic feedback [34]) is necessary. Control of multistability is possible by delay feedback as well, which can be combined with chaos control easily [66]. In that case, the feedback stabilizes one of the unstable orbits embedded in a chaotic attractor. Unfortunately, it just stabilizes an unstable orbit and does not target directly an already existing attractor. Parenthetically, the feedback control can stabilize an attractor against external noise by employing the Jacobian of the attractor as a feedback [67]. However, the state of the system must be already near to the desired orbit; that is, the application of one of the above targeting techniques is mandatory. If feedback control is possible then the aforementioned controls work very well; otherwise, the application of a non-feedback method (including a stochastic one) is inevitable. The main advantage of the nonfeedback control technique presented in this study is the ability to target the desired attractor directly, which is not possible with traditional non-feedback controls. For instance, the attractor selection by short pulses [31] (kick the system) introduces randomness to the control; that is, it cannot be guaranteed that the final state of the system after the kick is the desired attractor. Moreover, extremely long transients in the presence of transient chaos [68] are also possible. With the continuous transformation of attractors such transients can be avoided. The attractor annihilation techniques by pseudo-periodic forcing [69] or by harmonic perturbation of a state variably or a parameter [32] cannot target an attractor directly either. Roughly speaking, both techniques destruct attractors in the order of the size of their basin of attraction with increasing forcing or perturbation amplitudes.
The applicability of the method presented here has also requirements/limitations. Firstly, in many cases, there is the requirement that the system should not be modified appreciably during the control. Therefore, all non-feedback techniques apply only small perturbations. It is absolutely not fulfilled in the present method as large deviations from the original system parameters is usually necessary. However, by the end of the control, the system goes back exactly to the same parameter set and additional control is no longer required to sustain the desired state (attractor). It must be also noted that if the basin of attraction of the targeted attractor is small then slow transformation may be necessary.
Secondly, due to the large parameter variation during the control, detailed knowledge on the topology of the periodic orbits in a large parameter space must be available. This is usually a minor problem since the bifur-cation structure of many nonlinear oscillators have the same basic building blocks [70]. If there is no information about a state variable, the orientation in the topology can be supported; for instance, by measuring and Fourier transforming an indirect quantity connected to the dynamics (in case of bubbles their emitted pressure signal). The peaks presented in the spectra can provide valuable information about the nature of the actual attractor [71], see also the studies of Lauterborn and Cramer [72] and Lauterborn and Suchla [73] where complete period-doubling scenarios could be measured in this way of harmonically forced bubbles represented via spectral bifurcation diagrams. Another requirement related to the topology is that the transformed trajectories must be compatible in terms of winding numbers as already discussed in Sect. 8.
It must be emphasized that it is possible to combine the present non-feedback control technique with other methods. For instance, with stochastic control, attractors with small basins can be annihilated to "simplify" the state spaces during the transformation. Moreover, even if feedback control is possible in an application, the method introduced here is still a good candidate for targeting a desired attractor, and after the continuous transformation, the attractor can be stabilized against external noise with the feedback control. Similarly, the delay feedback control can stabilize an unstable orbit that replaces the originally stable chaotic attractor. After that, the continuous transformation technique can be applied between the resulting coexisting stable states if necessary and if it is possible. In this way, chaos control and targeting an attractor directly is accomplished simultaneously. There are case studies where the application of two different kinds of controls may be successful when the usage of only one of the control fails [74].
Since the control technique is demonstrated numerically only on a specific system (Keller-Miksis equation describing bubble dynamics) and on a very specific example, many questions arise: Is it possible to "bypass" somehow the winding number compatibility condition? Is it possible to transform a chaotic attractor into a stable periodic orbit and vice versa? Is this method applicable to other systems? In order to answer these question correctly, additional detailed analyses are mandatory; however, preliminary conjectures can be made.
Even if a continuous transformation between two coexisting attractors does not exist applying a single-frequency pair, it may be possible to find a transformation by using subsequent frequency pairs one after another, and reaching the desired attractor step-by-step hopping onto "intermediate" attractors. The transformation of chaotic attractors into periodic orbits fits into the large topic of chaos control [21]. Since chaos control is possible with the addition of a second sinusoidal component to the driving [24,75], the main issue is that how chaos is annihilated (a smooth inverse Feigenbaum cascade or a sudden crises) during the parameter tuning, and what is the resulting periodic orbit which is interesting from the winding number compatibility point of view. It is well known that chaos control is possible merely by adjusting the phase shift between the sinusoidal components of the driving [76]; therefore, the effect of the phase shift must also be thoroughly investigated. It can also play an important role between the transformation of the period-2 and period-3 orbits itself shown in Fig. 11. The period-2 solutions can always be transformed into one of the period-3 orbits through the blue or the green surfaces. In contrast, if the period-3 attractor is in a phase in time described by the red surface then transformation is not possible in the opposite direction (from period-3 to period-2), see again Fig. 11.
It is very likely that the control technique of continuous transformation is possible for other systems as well. As it was already mentioned, the topology of many nonlinear oscillators is composed by the same building blocks. Moreover, in recent publications, a phenomenon called replication of periodic domains is reported by applying a second sinusoidal component to the driving [40,77], or a single periodic perturbation in case of maps [78,79]. We believe that the replication of periodic domains and the branching mechanism introduced in Sect. 5.1 have the same dynamical origin.
Finally, it is worth mentioning that any application using dual-frequency driving can benefit (at least indirectly) from the results presented in this study due to the very detailed investigation. Some possible applications are acoustic cavitation [80][81][82][83], Faraday waves [84], stability of traveling beams [85,86] or laser-driven dissociation of molecules [87].

Conclusion
In the present study, a new non-feedback technique to control multistability in harmonically driven nonlinear oscillators is introduced. It is based on the observation that two attractors which cannot be exchanged by tuning the system parameters may be transformed continuously into each other by temporarily adding a second sinusoidal component to the driving. The process to target a desired attractor directly was tested at a relevant example, the Keller-Miksis equation-a secondorder ordinary differential equation describing bubble dynamics-and visualized via an animation (Online Resource 1). If this technique is proven to be applicable to other systems as well, it may become an excellent candidate to accomplish non-feedback control without long transients and/or randomness.