Acoustic vibrational resonance in a Rayleigh-Plesset bubble oscillator

Highlights • Vibrational resonance (VR) is reported in a Rayleigh-Plesset gas bubble model.• The complex dynamics is described by a classical particle in a time-dependent single or double-well potential.• Single and multiple VR can occur by the variation of the parameters of the acoustic waves and the properties of the liquid.


Introduction
Nonlinear dynamical systems are ubiquitous in nature and they exhibit a diversity of fascinating phenomena. One of these is resonance, which may occur in both linear and nonlinear dynamical systems [1][2][3][4]. The analysis of resonance, especially in nonlinear systems, has received close attention in the scientific community, both pure and applied. The fundamental understanding of resonance is highly relevant to numerous practical applications in many branches of science, engineering, and medicine [1,[5][6][7][8]. In general, resonance implies the matching of two or more frequencies within a system and, in most dynamical systems, it can appear in various forms when the frequency of an external driving force is matched to the natural frequency of the system, when it gives rise to an enhancement of the system's response [1,[8][9][10]. The type of resonance depends largely on the matching method as well as on the nature of the external driving force that produces it [11,12]. Vibrational resonance, for example, occurs in dual-frequency-driven non-linear systems with distinct frequencies [7][8][9][10][13][14][15][16][17][18][19]; stochastic resonance (SR) occurs when one of the two forces is replaced with noise of appropriate intensity [20][21][22][23][24][25][26][27][28]; and chaotic resonance takes place when the high-frequency component of the driving force is replaced with the output from a chaotic system [29][30][31]. Other forms of resonance include coherence resonance [32][33][34], ghost resonance [21,35,36] and autoresonance [37,38], parametric resonance [39][40][41], as well as conventional resonance in the form of harmonic, subharmonic and ultraharmonic resonances [42][43][44].
Our focus here is on the phenomenon of vibrational resonance (VR), originally proposed by Landa and McClintock [9]. In VR, a nonlinear system driven by a low-frequency signal is also subjected to a second, but high-frequency harmonic excitation, referred to as the fast-signal, such that the frequency () of the high-frequency component is significantly greater than the frequency ( ) of the low-frequency component. Under these conditions, the response amplitude of the system at the frequency of the low-frequency component versus the amplitude of the fast signal exhibits an "inverted-U-shaped" curve, similar to the signal-to-noise ratio curve observed in the SR phenomenon. By tuning the amplitude of the auxiliary signal or other parameter of the system, the weak/low-frequency excitation can be greatly amplified by VR. Basically, the fast signal acts as a control input to the system through the adjustments of its amplitude, frequency and phase (See Ref. [1] and references cited therein for further details). The phenomenon has gained enormous research attention in the last two decades and has been extensively investigated due its several potential industrial and biomedical applications in a wide range of contexts including bistable systems [45][46][47][48], multistable systems [19,49], systems with various forms of potential structures [13,15], linearly damped systems [1,50], https://doi.org/10.1016/j.ultsonch.2020.105346 Received 10 February 2020; Received in revised form 6 August 2020; Accepted 9 September 2020 electrical circuits and time-delayed systems [5,6,51,52], as well as plasma models [10,53]. Although earlier VR studies focused on the parameters of the high-frequency component of the bi-harmonic force, some recent VR research has explored the effects of parameters such as nonlinear damping/dissipation [7,10,53,54], the complementary roles played by system's innate bifurcation parameters [13,48,55], the effect of time-delay [52,56] and fractional derivatives [57][58][59], antiresonance [60], the depth of a potential and the location of its minimum [61], the effects of potential deformation [62] and potential roughness [53], and mass-dependent VR [63], amongst others. More importantly, VR has been realized experimentally in vertical cavity surface emitting lasers and optical systems [45,46,49,64,65]. The aforementioned investigations have enriched our understanding of the VR phenomenon, its associated mechanism, and the roles played by relevant system parameters in inducing or promoting VR, as well as in suggesting potential applications.
Some intriguing properties, novel phenomena and potential applications have been explored and reported. For instance, Gandhimathi et al. [66] examined the relationship between SR and VR, while Morfu and Bordet [67] recently established a connection between phaselocking modes and VR. Rajamani et al. [68] observed a novel ghostvibrational resonance phenomena in a system driven by multi-frequency signals and extended their analysis to a network of oscillators. A novel vibrational antiresonance in which a strong signal-suppression occurred due to a large negative phase-shift in coupled oscillators was very recently reported by Sarkar and Ray [60]. Furthermore, the existence of subharmonic and superharmonic resonances [69] and higher harmonics [70] have also been reported. In a different sequence of developments, new methods of analysis in which nonlinear systems are driven by aperiodic forces, in contrast to periodic driving, have been investigated and has given birth to what is termed aperiodic vibrational resonance. There are a number of its variants, including re-scaled aperiodic vibrational resonance and the double-sampling aperiodic vibrational resonance methods [71,72]. In a very recent paper, a new approach to VR in terms of a spectral amplification factor was presented [73]. All these results are directed towards exploring the diverse potential applications of VR in, for instance, improving energy harvesting from electromechanical vibration devices [74][75][76] as well as in enhancing generalized energy detector approaches to weak random signals detection [77]. Another important application is concerned with signal processing. In this direction, VR has been proposed for the transmission, filtering and amplification of signals; such applications have been demonstrated both theoretically and experimentally [78][79][80][81]. In rotational machinery quite generally, bearings play crucial roles e.g. in mining, aerospace, railways, and metallurgy. When they become faulty, the equipment typically fails, sometimes catastrophically. In this context, VR has been proposed as a promising technique for the early noninvasive diagnosis of nascent bearing faults [82][83][84][85][86].
In the light of the forgoing, we investigate and analyze the VR phenomenon in a model of an amplitude-modulated, acousticallydriven gas bubble in an incompressible liquid which, to the best of our knowledge, has not previously been discussed. Bubbles appear in both the life sciences and natural sciences, as well as in technology [43], and their dynamics has received considerable attention in connection with acoustic cavitation [42][43][44] where bubbles are excited by an acoustic sound field mostly consisting of single or dual frequencies. In addition to the existence of bifurcation structures and chaotic behaviour, bubble dynamics is rich in special classes of resonance, and these have constituted one important focus of recent research. For instance, the nonlinear dynamics of bubbles under dual-frequency acoustic excitation was numerically investigated recently within a wide range of parameters by Zhang et al. [44]. Two significant resonance features, namely, combination resonance and simultaneous resonance were observed. Moreover, parameters such as the bubble radius, the acoustic pressure amplitude, and the energy distribution between wave components were reported to be of crucial importance in that they impacted significantly on the nonlinear bubble oscillation. Under multi-frequency acoustic cavitation, the dynamics of bubble oscillations also shows evidence of resonances [44]. However, despite the huge number of investigations of bubble oscillations, and the wide range of potential applications of multi-frequency excitations of gas bubbles, the dynamics of gas bubbles driven by dual-frequency forces is yet to be fully investigated and understood. In sonochemistry, fluid engineering and medicine (in particular, medical diagnostics and therapy), where acoustic cavitation due to bubble formation is well-known, a multifrequency acoustic sound field promotes system effectiveness and can also be useful in controlling the chaotic dynamics of the bubble oscillations, especially where they are undesirable [2,87,88]. In addition, the dual-frequency approach is employed in designing sonochemical reactors [89], synthesis with nanoparticles [90] and sonoluminescence [91], and for investigating bubble size distributions with void fraction in the open ocean [92], in studying the dynamic variations of bubble radii during oscillations [93,94], and in the enhancement of the ultrasound accuracy of biomedical diagnosis [95,96], and the efficiency of tumor therapy, including tumor ablation [97]. The dual-frequency approach can also accelerate bubble growth through mass transfer, leading to the generation of large bubbles [44,96,98]. We draw attention to a very recent study of dual-frequency-driven bubbles in a highdimensional parameter space based on GPU accelerated computations [99]. It represents the state-of-the art in applications of dual-frequency irradiation and also provides an extensive review ranging from increased sonochemical yields to decreased sonochemical efficiency due to the synergetic effects of dual-frequency driving. Consequently, there is a compelling need for further investigation of the nonlinear responses of dual-frequency driven oscillators to the effects of high-frequency components in the context of the VR of gas bubbles driven by an amplitude-modulated acoustic field. We undertake such an investigation here. In addition to the parameters of the modulating force, and on account of the multi-parametric nature of the bubble oscillator, the possible contributions of the parameters of the bubble and its surrounding medium (in this case liquid) to the occurrence of VR are investigated as they can in principle provide information essential for addressing other problems arising in acoustic cavitation. We provide numerical evidence for the occurrence of VR as well as for its enhancement by amplitude modulation over a wide parameter range. Our results are validated theoretically based on the method of separation of time-scales. The rest of the paper is organized as follows. In the following section (Section 2), we present the Rayleigh-Plesset bubble oscillator model -expressing the complex equation as the dynamics of a classical particle in a potential well of the Liénard type [100] for the first time. We provide evidence that, using the Rayleigh-Plesset bubble model, an acoustically-driven bubble oscillates in a timedependent single or double-well potential whose properties are determined by the density of liquid and its surface tension. Section 3analyzes and discusses the stability of equilibrium points. In Section 4, detailed analytical treatments of the resonance behaviour are given. Numerical simulation results, complemented by theoretically computed results are presented in Section 5. The paper is concluded with the identification of some highlights, in Section 6.

The model
The model of interest is the Rayleigh-Plesset equation for smallamplitude oscillations of a spherical gas bubble in an infinite incompressible liquid [101]. In the model equation, the bubble is treated as an oscillatory system (bubble oscillator) that incorporates various parameters of the surrounding liquid. In general, the size of the bubble and its gas/vapour contents are specified. Thus, for a spherical gas bubble, the defining variables for its size and shape can be consolidated into the time-dependent bubble radius R t ( ). In addition to R t ( ), other essential parameters include the equilibrium radius (R 0 ), the external pressure on the liquid (P ext ), and the bubble's internal pressure P in , as well as the polytropic exponent of the gas in the bubble ( ), the density of the liquid ( ), the liquid and thermal viscosities (µ µ , and l th ), and the surface tension ( ) of the liquid.
When driven by acoustic fields, the bubble dynamics may be likened to that of a damped-driven oscillator and, during the bubble oscillatory motion, the contained interior gas behaves polytropically. The polytropic exponent describes the thermodynamics state of the gas in the bubble. It depends on the acoustic frequency of the driving field, the bubble radius, and the nature of the gas. In general, if the bubble is sufficiently large or the driving frequency is sufficiently high, the bubble will behave adiabatically because there will be a significant temperature gradient across its radius [102]. In this case, the ratio of specific heats. If, on the other hand, the bubble is sufficiently small or the driving frequency is sufficiently low the bubble will behave isothermally. Intermediate thermodynamic states between these two extremes have been identified [102]. In principle, the VR phenomenon requires a very high-frequency driving force. It is therefore appropriate to assume the adiabatic state, so that = 1.3.
Moreover, it has long been known that the natural oscillations of a bubble are dependent on the effective damping mechanisms arising from the liquid and thermal viscosities. There is also a viscosity due to acoustic radiation to be considered, as well as damping due to mass diffusion and to the evaporation-condensation of vapour [103]. In experimental studies, however, damping due to mass diffusion and to the evaporation-condensation of vapour are often neglected [102]. While the other sources of damping are clearly understood, the origin of thermal viscosity is not straightforward. Unlike purely isothermal processes in which pressure and volume changes are mutually in phase, a phase difference exists in real bubbles and the gas can behave adiabatically. The latter processes are associated with a net flow of heat energy into the liquid when work is done on the gas, reducing its volume by increasing the pressure during compression. This work done on the gas is far greater than the work done by the gas in moving the surrounding liquid during expansion, thereby increasing the internal energy of the gas as well as transferring energy by conduction through the gas into the surrounding liquid. The net flow of heat energy is characterised by a thermal viscosity that has been reported to influence the free oscillations of the bubble strongly, particularly during the transition between the adiabatic and the isothermal states [103][104][105]. Recently, the impact of non-linear damping/dissipation on the occurrence of VR in nonlinear systems has also gained considerable research attention [7,10,53,54]. In view of this, and to account for the roles of all the individual factors contributing to damping of the bubble dynamics, and their influence on the VR phenomenon, it is essential to treat the liquid and thermal viscosities, as well as the acoustic radiation, independently in the derivation of the equation of motion. This contrasts with the traditional approach of defining an effective viscosity whose effective contributions may appear negligible in the VR analysis.
Here, we consider the Rayleigh-Plesset equation for bubble oscillations which can be written [43,101,[106][107][108]; Here P 0 is the ambient pressure, and and are the amplitude and angular frequency of the external acoustic excitation. To facilitate theoretical and numerical investigation of VR, it is useful to express Eq.
(1) as the classical dynamics/motion of a particle in a potential well. This enables us to elucidate the nature of the potential within which a bubble oscillates -an important issue in the study of dynamical systems that was overlooked in previous research on bubble dynamics. Thus we will, for the first time, obtain an exact expression for the potential well in which the bubble oscillates in the context of the Rayleigh-Plesset bubble model. We will follow the steps in [109], expressing the bubble radius R, in terms of a dimensionless variable x, i.e. = + R R x (1 ) 0 , and write P in as Eq. (2) is the continuity condition for the surface tension across the bubble interface, governing its forced radial oscillations [103].
is the pressure at the bubble wall in the liquid satisfying the continuity condition, and P p R R t ( , , ) 0 is the dimensionless deviation from the equilibrium pressure. For small-amplitude oscillations one considers thermal effects in the bubble content leading to strong damping on the transitional behavior of the bubble. At Blake's critical threshold [110], the equilibrium radius of the bubble depends only on this tension and, assuming the equilibrium radius of the bubble is in the neighborhood of the critical radius, the derivative of the tension with respect to the equilibrium radius of the bubble can be approximated as in dimensionless form [103,111,112]. With where l 0 . Using a binomial expansion, we can express the multiplicative term + x (1 ) 1 as : retaining up-to the quadratic term. After some algebraic manipulation, it is easy to show that Eq. (7) can be written as: Substituting for P in eq , in Eq. (4) into Eq. (8) and using the definition of F acc yields such that the Rayleigh-Plesset equation in Eq. (1) can be written in dimensionless form as cos , The equilibrium radius of the bubble, R 0 may take on values ranging from a micron or less (i.e. nanobubbles and microbubbles) up to several millimetres [42,43,111,113,114]. For underwater sparks or with explosives, the radius of the bubble may reach even larger values [43]. The dimensionless bubble parameter are given in Table 1  atm. Finally, we express Eq. (10) in terms of the bubble's potential: where V x ( ) is the bubble's potential .
We find that the Rayleigh bubble oscillator is a strongly nonlinear, parametrically excited, system characterized by higher-order nonlinear space-dependent dissipation of the Liénard type III [100], typical of fluidic plasma oscillators [115]. Parametric excitation arises from the Table 1 Calculated values of the system's dimensionless parameters for an equilibrium radius = R 10 0 mm. The values of the constants used are given in the text. 145.38 modulation of the system's potential when irradiated by acoustic fields, t cos of amplitude , so that the potential V x t ( , ) is time-dependent. Note that previous works linearised the Rayleigh-Plesset equation (Eq. (1)))) to that of a harmonic oscillator by neglecting all powers of x and of its derivatives higher than the first (see e.g. Refs. [102][103][104]112]). The corresponding truncation of higher-order nonlinearities in the potential and the damping functions impacted significantly on the nature and number of equilibrium points.
In the absence of the external acoustic field, the bubble would oscillate freely in the time-independent potential, The shape of the potential V x ( ) given by Eq. (14) describes the equilibria of the system and depends on the signs and values of the three dimensionless parameters, , and as we will show in the following analysis. (See the Appendix I for mathematical details).

Equilibrium and stability
To analyze the potential structure, we first revisit briefly the wellknown theory of bubbles dynamics in relation to our model as shown in Fig. 1(a) [116][117][118][119][120][121]. Accordingly, the oscillation of gas bubbles can be classified into three distinct regions relative to the special threshold which we summarise below.
(I) If the ambient pressure is higher than the vapour pressure of the liquid, the freely oscillating bubble posses only one stable equilibrium state, denoted as case I in Ref. [120]. (II) When the ambient pressure falls below the so-called Blake's critical threshold or is in its neighbourhood, there exist two equilibria, one of which is stable while the other is unstable. This condition was classified as case II in Ref. [120]. (III) The trajectory of the system could tend to infinity whenever the ambient pressure is lower than the Blake's critical threshold so that the dynamics is neither stable nor unstable, and without any local/ global maxima/minima. This condition was classified as case III in Ref. [120].
Remarkably, case II and case III are separated by the Blake's critical threshold which is a special limit that determines the stability boundary of the bubble system and also specifies the limit where bubbles tend to grow infinitely. Given the approximations employed in our model, our results can be considered to be in agreement with the existing literature [116][117][118][119][120][121].
Using our proposed bubble oscillator model, we plot the dependence of the dimensionless bubble radius x on the dimensionless potential parameter which controls the potential structure as shown in Fig. 1(b). It is obvious that the three distinct regions enumerated above are predicted by our model in a similar manner to the way in which the ambient pressure was used in previous works. Thus, one of the novelty of our model is its ability to employ the critical ambient pressure ( = ) estimated at the Blake's threshold by equating = from Eq. (11) in predicting the potential well conditions, thereby providing an alternative approach to the exploration of bubble dynamics. Notably, the value of P cr 0 depends on the liquid surface tension which can be appropriately controlled [120], the equilibrium bubble radius (R 0 ) and the thermodynamic process defined by the polytropic parameter .
If the ambient pressure P 0 is greater than the estimated critical pressure ( > P P cr 0 0 ), the dynamics of the system satisfies case I corresponding to a single-well potential condition for < 0. Furthermore, depending on the values of P 0 , there is a narrow regime where case I (for < 0) overlaps with case II resulting in the coexistence of single-well (case I) and double-well (case II) oscillations which suggest the existence of multistability of attractors recently reported in coupled bubbles [122]. In addition, the dynamics when = 0 or slightly greater than zero does not translates to the dynamics for = P 0 0 . This is obvious in Fig. 1(a), and considering that is also dependent on the ambient pressure P 0 . The double-well condition ( > 0) which we termed case II arises if < P P cr 0 0 , in agreement with the existing literatures [116,120]. Finally, when becomes extremely large such that , every trajectory of the system tends to infinity, and the case III, in which there is no equilibrium state is satisfied [116,120].
The completeness of the proposed prediction is clearly summarised in Figs. 2 and 3 where we have plotted the potential V x ( ) for different values of , and . A single-well ( < = = 0, 0) or a single-hump ( > = = 0, 0) as shown in Fig. 2(a) illustrates the special case of V x ( ) in which = = 0, corresponding to a harmonic oscillator, d and 0 being respectively the damping parameter and square of the linear resonance frequency investigated earlier [102][103][104]112]. It admits only two types of potential structure satisfying the general case I -a single-well for < 0, and a single-hump for > 0.
When higher nonlinear terms are considered two equilibrium states appear. The potential V x ( ) can be a single-well-single-hump as depicted in Fig. 2 0, 0, corresponding to one stable equilibrium located at = x 0 and one unstable equilibrium state located at Fig. 3(a). Previous reports had shown that the system potential may consist of a smooth background potential plus a rapidly and randomly varying potential that may alter the dynamics of the system when the particles are insufficiently energetic to overcome the net potential barriers. When a dynamical system with energetic particles encounters such potential, the system is capable of annulling the effects on its trajectories [53]. This case of > 0 (double-well-double-hump potential) is therefore consistent with the stable and unstable equilibria depicted as case II provided that the potential barrier at the hump (H 1 ) located at < x 0 is significantly higher than the barrier (H 2 ) near = x 0. Under this condition, energetic particles can subdue the potential barrier created in the neighbourhood of this unstable equilibrium (H 2 ). The inset, Fig. 3(b) depicts the zoom of the asymmetric structure of the potential well, showing two wells located near = x 0.10 and = x 0.35 and sandwiched by relatively small barrier at = x 0. Fig. 3(c) shows that our model satisfies case III for sufficiently low ambient pressure below Blake's threshold, showing no equilibrium point. In this case, the particles of the bubble system gain sufficient energy to overcome all the potential barriers due to the infinitely increasing bubble radius as . The discussions above can be analytically justified based on theoretical verification of the stability properties using the Jacobian matrix. The unperturbed nonlinear Rayleigh bubble oscillator derived from Eq. (10) is written as: The equilibrium points, P x y ( , ) i of the autonomous system (15), where x and y are the fixed points are obtained by equating the L.H.S to zero and solving the resultant algebraic equation. Thus, This yields the following four equilibrium points: Remarks 3.1. We make the following remarks on the equilibrium points Depending on the sign of , two cases exist for the equilibrium points P P , It is conjectured that the system can also admit additional equilibrium solutions depending on the choice of truncation of the binomial expansion of + x (1 ) 1 in connection with higher-order nonlinear terms of the potential functions. Research shows that up-to-triple well solutions are possible in higher-order nonlinear systems in contrast with the familiar single and double equilibrium solutions reported earlier for lower-order potential functions [53,115,[123][124][125]. Similarly, a flip-flop between hard-spring and soft-spring bistabilities due to higher-order truncation of the Toda oscillator was observed and analyzed by Goswami [126] while, earlier, a third-order approximation reduced the Toda oscillator model to the Heńon-Heiles (HH) Hamiltonian system which is non-integrable in contrast to the original integrable Toda Hamiltonian [127][128][129]. Although the nonlinear dynamics, including the bifurcations and the evolution of diverse orbits as well as details of the equilibrium points of our derived bubble oscillator, is beyond the scope of the present paper, we comment that the features of our derived bubble potential function are consistent with the previously reported results [116][117][118][119]130,131,121].

Acoustic vibrational resonance
Now, we return to the main focus of this paper. The bubble oscillations can be highly complex and chaotic. We note that the media and the type of application determine parameters such as the density, viscosity, surface tension, and diameter of the bubble [132]. To investigate the occurrence of VR in Eq. (10), a higher amplitude and fast periodic driving force is coupled to the acoustic sound field, such that the low-amplitude acoustic field is modulated by a high-frequency cosine signal, f t cos . Under the exposure of the high-frequency modulating signal, a single bubble emit light -a phenomenon known as single bubble sonoluminescence (SBSL). Thus, our choice of high-frequency modulation signal stems from their roles in SBSL and other realworld bubble cavitation experimental applications in which only certain parameters are available for perturbation. Hence in practice we can only expect to perturb the external forcing terms. Parameters such as the surface tension and liquid viscosity are strictly determined by the host medium and therefore not available for perturbation since it would be impracticable e.g. to open up the skin during high-intensity focused ultrasound (HIFU) tumor treatment to change the properties of the tumor medium [132]. Here, we define , and its amplitude is greater than that of the weak driving acoustic force. In the following analysis, we denote by the coefficient of the quadratic damping force in Eq. (10) to obtain,

Analytical description of VR
We now employ the method of direct separation of the dynamics, to separate the motion of the bubble system into fast and slow motions. Thus, we obtain a set of integro-differential equations consisting of the systems' equation of slow motion whose response can be varied by simply adjusting the parameters of the high-frequency driving input. The response amplitude, Q, defined as the ratio of the amplitude A L to the frequency f is obtained by solving the equation of slow motion. Thus the solution x t ( ) of the bubble oscillator given by (18) is assumed to be a superposition of the solutions t ( ) for the slow dynamics whose frequency is and of t ( ) for the fast motion with frequency with , in the form: Here t ( ) is a periodic function of period = T 2 and is a periodic function in the fast time = t with period 2 . The mean value of with respect to fast time is given by The next step is to derive a system of two coupled integro-differential equations for the variables and from the main equation of the system (18) Eq. (25) is the sought-after analytic expression for the slow oscillation of the bubble, with the parameters of the fast signal. Eq. (25) will be employed in computing the theoretical response amplitude Q of the bubble system at the lower frequency . To obtain the equation of fast motion, we subtract Eq. (25) , and .
Thus, we can write Eq. (28) as where V ( ) eff , the effective potential of the system, is given as Hence, by neglecting the nonlinear term in Eq. (34), we obtain a linear system written as which can be re-written as where, = C t ( cos ) such that,  Fig. 3(a) for zero and small amplitudes, f. With increasing f the shape of the potential changes and, for = f 335, it becomes a single well with no resemblance to the system's potential presented in Fig. 3(a). It is evident that the effective potential of the slow motion depends on the amplitude f and frequency of the fast motion. In the absence of a high frequency driving force, i.e. by setting the amplitude = f 0 in Eqs. (26), then Eq. (31) reduces to Eq. (13), the model's potential.
One of the most striking results of this paper is that there is also a contribution to the effective resonance frequency of the system due to the parameters of the fast-oscillating force as given by Eq. 4 . This contribution was either not found, or more usually overlooked, in most earlier analyses of VR. Thus, our results here differ significantly from previous ones, where VR was found in systems with linear and nonlinear dissipation, such as the asymmetric Duffing and plasma oscillators, in which VR has been extensively studied. For small amplitude oscillations, the dynamics of the oscillating gas bubble explicitly depends on the acoustic damping constant and the natural resonance frequency of the bubble motion [44,111,133,134]. We remark at this point, that our interest in the parameter of the natural resonance frequency of the system is that it principally defines the various types of bubble resonance (e.g. main resonance, harmonics, subharmonics and combination resonances) observed under multi-frequency acoustic excitation [44]. In addition, it dictates the shape of the potential wells of the system.
We will consider both the single-well and double-well configurations of the system's potential, noting that in the latter case there is the possibility of observing multiple resonance peaks, and bearing in mind that VR occurs in systems with bistable and multistable potentials where for certain parameter regimes multiple resonance can be obtained [1,13,14,18,55,61,62,135]. Note that the effective resonance frequency can play a direct and significant role analogous to the effective nonlinear dissipation reported in Refs [7,10,54,53] in the enhancement of signals by modulating the parameters of the fast signal. The effective resonance frequency parameters, which can be computed from the second term of Eq. (31), therefore dictate the equilibrium points of the slow motion. When = 0 in Eq. (31), the effective potential function reduces to,

Numerical description of VR
For full numerical integration of the nonlinear Rayleigh bubble oscillator, it is expedient to write Eq. (19) as a set of two coupled, firstorder, autonomous, ordinary differential equations (ODEs) of the form; where the response Q of the system can be computed from the expression such that the response to the frequency is defined as the amplitude of the sine and cosine components of the output signal, and the numerical values of A S and A C are related to the Fourier spectrum of the time series of the variable x computed at the frequency as

Analyses of resonance
We now analyze vibrational resonance in the bubble oscillator under different sets of conditions. The amplitude of the system's response is = Q S 1 ; where = + S C ( ) r 2 2 2 0 2 2 from Eq. (38). In the linearized equation of the system, the linear dissipation C 0 appearing in Eq. (36) is the linear damping coefficient which, for a given liquid of density , is dependent on both the frequency , and the amplitude f, of the fast signal, as well as on the fundamental parameters of the medium (i.e. the thermal (µ ) th and liquid (µ l ) viscosities, respectively). Moreover, the contribution of the damping force is negligible in determining the dynamics of the system compared to the effect of the liquid's surface tension, the related natural frequency of the system ( ); and variations in the parameters of the fast signal (g or ) produce the same effect as . It follows that a variation in the parameters on which the natural frequency of the oscillating gas bubble depends would produce additive and multiplicative effects on the appearance of VR in the system in accordance with C 0 . We set = W r 2 2 , so that = + S W C 2 0 2 suggests that Q will attain its maximum when S is minimum, that is, when = W 0 (such that, . Consequently, resonance will occur. For ,the value of becomes too large (i.e. since ) so that C 0 0 , and the system reduces to a linear differential equation with no dissipation, implying that VR does not occur.
Since the acoustic wave is weak ( 1), and the amplitude of the bubble oscillation is considered to be small [111] Eqs. (31) and (39) enable us to analyse the conditions for the occurrence of VR in each case. This shows that single, double, or multiple VR peaks can occur in the system under both conditions. Moreover, the appearance of the resonance peaks could be achieved by adjusting either the parameters of the fast signal (f or ), or the parameters on which the natural frequency of the bubble oscillator depends .

Resonance with a single-well potential
We fix parameter values as follows: and = 1.5. We choose < < 0, 0 and > 0, for which V x ( ) is a single-well potential as shown in Fig. 3(a). In order to examine the response amplitude of the system, Eqs. (41) was computed from the numerically integrated Eq. There is a marked increase in response amplitude as the value of decreases. In Fig. 5, the resonance peaks appear when > f 500. This is because the resonance frequency, 0 r when the value of is appreciably large since and < 0, showing that the enhancement of the systems response can also be achieved by modifying either the surface tension of the liquid or its density , both of which are dependent on , as shown by Eq. (10). To further illustrate the occurrence of VR for the conditions of single-well potential, the response Q is plotted against f in Fig. 6 for different values of = 0.5, 1.0, 2.0, and 3.0 with 245, 0.01, 145, 2 and = 32.9. In this case ( Fig. 6(a)), the system also responds significantly to the variation of , even at the lower values, e.g. = 0.5. Here, the resonance peaks occur at lower values of f but with smaller magnitudes compared to Fig. 5. Noticeable also is that increase in the values of the frequency of the acoustic signal shifts the peaks rightwards to higher values of f, rather than creating more peaks as depicted in Fig. 6(b).

Resonance with double-well potential
For the double-well case, subject to the condition: > < 0, 0 and > 0, Fig. 7 shows the variation of Q as a function of f for different values of , (145, 200, 270 and 300) with other parameters fixed. Here, Q exhibits a number of resonance peaks, with two distinct peaks which approaches the limiting value of 0.4. The response amplitudes Q in this case are smaller in magnitude than those in Fig. 5. However, the resonance occurs at lower values of f, including the second peaks. Notably, an increase in the value of produces three remarkable effects: (i) a shift in the resonance peaks towards lower values of f; (ii) a reduction in the system's response; and (iii) elimination of one of the resonance peaks. This is further illustrated in Fig. 8 where we have plotted Q against f for other values of and for a different set of the bubble's parameters. As becomes larger, typically when 550, the small resonance peaks that manifest themselves within the main resonance tend to disappear -a feature that is fully captured in Fig. 9 where Q is plotted against in the weak f regimes. For all f values shown, two distinct peaks are clearly visible when < 400. However, for > 400, the bubbles cease to vibrate. This can again be inferred from the expression for r at resonance. At resonance, induces new resonance peaks. For = 10, only one pronounced resonance peak occurs as shown in Fig. 10(a). In the double-resonance cases shown in Fig. 10(b and c)  ). In Fig. 10(d), when = 18, three resonance peaks appears with increased response. It has earlier been reported that in pure mechanical systems vibrational resonance occurs in the overdamped case due to minimization of the resonant frequency r , while in the underdamped cases, the mechanism is related to local minimization in r 2 2 [13,15]. This can arise either when the resonant frequency is tuned to match the lowfrequency component of the input signal, or when the input lowfrequency component matches the resonant frequency of the system. In contrast to these earlier reports, for the acoustically excited bubble system considered here, the resonant frequency r gets locally maximized via as a result of the double-well potential condition. This follows from the fact that the resonant frequency is dependent on both f and .
To examine the effect on the response amplitude Q of various properties of the liquid in which the bubble oscillates, we first illustrate the dependence of Q on for four different values of the liquid density related parameter ( 6.5, 7.0, 7.5, and 8.0), as shown in Fig. 11. Here, it is clearly evident that increasing the value of increases the system's response. The import is that, for constant ambient pressure P 0 and fixed equilibrium radius of the gas bubble, the bubble's response to high-frequency modulation is strongly dependent on the liquid density. The response is large when the density of the surrounding liquid is small, and vice versa. Furthermore, in Fig. 12, the response amplitude is shown as a function of for different values of f. Evidently, at the lower values of f ( = f 1.5), the system's response is greatly enhanced, implying that the nature and properties of the liquid can be utilized effectively in optimizing the response of the system. Although some properties of the liquid are constrained, nevertheless, in real world applications involving cavitation phenomenon, where the initial bubble radius R 0 , the pressure P 0 and the fundamental frequency are all determined by the type of application, consideration of appropriate liquids with specific properties (i.e. viscosity, surface tension and density) would provide a way of optimizing the system response.
To further analyse the occurrence of VR in the system subject to the double-well contitions, the dependence of Q on is plotted in Fig. 13 (42), are indicated by marker points. much less than 3. For = 5 ( Fig. 13(b)), the system resonates and further increase in the value of the low acoustic frequency component increases the response amplitude as well as inducing more resonance peaks. Note that this feature differs significantly from the behaviour found in the case of the single-well potential. This signifies a strong dependence of the bubble oscillator's VR on the surface tension of the confining fluid medium as well as its density, both of which are related to the normalized resonance frequency parameter, . It implies that the occurrence of VR in a bubble system depends directly on parameters defining the natural frequency .
The curved surface of the bubble oscillator's response amplitude Q as a function of the natural resonance frequency , and the amplitude f of the high-frequency acoustic field, are shown in Fig. 14. In addition to the effects of the parameters of the high-frequency forcing, one can immediately observe that the effect of high-frequency acoustic field is optimized by the parameter to enhance bubble response. This enhancement is significantly pronounced when the surface tension-dependent parameter ( ) is such that > 100 and < f 200. The occurrence of resonance induced by the parameter implies that system's response can be controlled by altering . In practical terms, the liquid's surface tension may be altered by adding impurities to the liquid. Fig. 15, shows a three-dimensional (3D) plot of the dependence of the system's response Q on the frequency of the low frequency acoustic field, and the natural frequency of oscillations . Evidently, multiple resonance peaks can occur with appropriate choice of the parameter values for the high frequency acoustic driving force and, in the low frequency range: , with 10 20. The 3D plot is characterized by two distinct and well-pronounced mountains separated by a valley wherein Q takes on near-zero values. While the bigger mountain occurs in the parameter regime < < < < 10 20, 100 300, the smaller mountain occurs in the neighbourhood of < < 15 20 and < < 7 100. As varies, the resonant response Q also varies. For 20, the peaks widen by increasing the spaces between them, together with a drop in the system's maximum response. This is similar to the earlier reported vibrational resonance mechanism in a quintic Duffing oscillator with multiple potential wells [14,15]. To complete the discussion    we refer to Figs. (14) and Fig. (15) illustrating the occurrence of multiresonance for 145, which is in good agreement with our condition for the occurrence of resonance in the system obtained from the conditions imposed on the potential wells. This has earlier been established in mechanical systems where resonance behavior is dependent on the shapes of the system's potential [8,10,[13][14][15][16][17][18][19]. We remark that regimes of strong resonance mountains can also be observed even at lower values of the amplitude f when the response Q is plotted as functions of both the perturbation parameters in the range ( giving rise to some narrow response regime for > 10 , and consequently will lead to the further disappearance of the resonance peak.

Conclusion
We have examined the oscillations of an acoustically-forced gas bubble in an incompressible liquid, using the Rayleigh-Plesset bubble model. We showed that the bubble oscillation could be expressed as the dynamics of a particle in a time-dependent single or double-well potential whose properties are determined by the density of the liquid and its surface tension. The contribution of the acoustic field to the potential of the system was analysed. We considered the case of a dual-frequency acoustic driving force, consisting of a low frequency acoustic wave and a high frequency force as amplitude modulator. We investigated the occurrence of vibrational resonance (VR) numerically and validated our results theoretically, and found excellent agreement. In our analysis of VR phenomenon, both the the single-and double-well potential configurations were treated, and the conditions for the occurrence of VR obtained. For the single-well condition, the system showed a significant response on variation of the system parameter , associated with the    liquid density. In this case, single resonance peaks were observed at higher values of f ( f 500) -the driving amplitude. In contrast, using the double-well potential configuration it was observed that, unlike the well-known traditional VR theory reported in the literature [8,9,45,49,136], amplification of the system's response is achieved by varying , instead of the parameters (f and ) of the high frequency acoustic force. This stems from the fact that the effective resonance frequency of the bubble oscillator depends on both the natural resonance frequency of the system and the parameters of the high frequency acoustic field. In addition, it was found that the appearance and disappearance of VR in the bubble system can be controlled by appropriate choice of the low frequency, .
Remarkably, multiple/dual-frequency bubble systems have many advantages over their single-frequency equivalent, with an abundance of practical applications as enumerated in the introduction. In contrast to numerous studies where multiple or dual-frequencies were achieved by additive application of external acoustic waves to the bubble oscillator, with the attendant need for large ancillary equipment to provide the external fields, dual-frequencies can be achieved instead through amplitude modulation of single-frequency acoustic waves. The acoustic vibrations of the bubble oscillator demonstrated here can be exploited to facilitate the nucleation of bubbles in biological tissues, enhancing the bubble growth rate as well as generating extremely violent forces during bubble collapse. The growth processes occur during cavitation and can be controlled by the parameters defining the properties (e.g. viscosity and surface tension) of the medium in which the bubble oscillates and the characteristics of the applied ultrasound field (amplitude and frequency). Specifically, control of the applied ultrasound energy can be achieved through the dynamic multifrequency radiation force usually created using ultrasound equipped with time-dependent high-frequency periodic amplitude modulation [137]. Thus, our simulated results could have direct applications to the production of polymer foams and to the ultrasonic degassing of liquids, as well as to the fabrication of porous ceramics -all of which undergo bubble nucleation events. Indeed, an amplitude modulation mechanism has been experimentally demonstrated recently for efficient enhancement of acoustically-driven micromixing and microcentrifuging processes, which has the advantage of eliminating the requirement for large ancillary equipment necessary for the provision of multiple external fields [138]. Finally, in biological applications, the parameters of the ultrasonic fields can be adjusted selectively to ensure that damage to larger vessels is significantly reduced or to favour the retention of extracellular matrix;, which can be achieved by modulating the amplitude of the acoustic wave using high frequency modulating forces.
The fixed points x y ( , ) of the autonomous bubble system are obtained by setting with the discriminant = 4 ( ) 2 . The nature of the fixed points ± x is determined by the sign of and can be examined for < 0 and > 0.
A. When < 0, there are six possible cases. In summary, it follows from A and B that = ± ± x 2 are two non-trivial fixed points and the two solutions + x and x corresponds to two equilibrium points P 3,4 ( ± , 0 2 ). Thus, the unperturbed system has four equilibrium points: + P P P P 0, 0 , , 0 , 2 , 0 , 2 , 0 .  indicating the concurrences of our numerical and analytical results.