Delay-induced vibrational resonance in the Rayleigh–Plesset bubble oscillator

We examine the impacts of time-delay and phase shift between two acoustic driving forces on vibrational resonance (VR) phenomena in the oscillations of a spherical gas bubble. Using the approximate method of direct separation of the motions, we obtain the equation of slow motion and the response amplitude, and we validate the theoretical predictions with numerical simulations. We find that the response amplitude of the system at the lower frequency varies periodically with respect to the phase shift. When the phase shift consists of an even number of periods, it can be optimized to enhance the system’s response in the relevant parameter space of the high-frequency driving force. In addition to the enhancement of the VR peak by variation of the phase shift, our results show that the time-delay also plays a significant role in the bubble’s response to dual-frequency acoustic driving fields. It and can be exploited either to suppress drastically, or to modulate, the resonance peaks, thereby controlling the resonances. Our analysis shows further that cooperation between the time-delay and the amplitude of the high-frequency component of the acoustic waves can induce multiple resonances. These results could potentially be exploited to control and enhance ultrasonic cleaning processes by varying the time-delay parameter in the presence of phase shifted dual-frequency acoustic waves. Moreover, it could be employed to achieve improved accuracy in ultrasonic biomedical diagnosis and tumour therapy, as well as for targeted delivery of reagents transported within bubbles.


INTRODUCTION
Nonlinear systems driven by external forces are abundant in both the natural and technological contexts [84], and many of them are known to be time-delayed.They range from physical, electronic and communication engineering systems to ecological and epidemiological models [5,26,48,61,65].A nonlinear dynamical system is said to be time-delayed when the evolution of its state variable is dependent on part of its history.Time-delay may appear in different forms due to processing delays leading to memory effects, or as a result of finite velocities of signal propagation, and it has been a subject of intense research for a long time [48,65]; the existence of a phase shift between a system's periodic response and the periodic force(s) driving it is of course well known in physics, including in stochastic resonances [25].
In the analyses of this interesting class of systems, the effects of time-delay and phase shift have been investigated independently.On the one hand, time-delay can indeed impact significantly on the system's dynamics and general efficiency, leading to the occurrence of intriguing varieties of dynamical phenomena including, but not limited to, strange nonchaotic attractors, hyperchaos, novel bifurcations, amplitude and oscillation death, stochastic dynamics, synchronization enhancement and pattern formation [48,65].On the other hand, the transmission, detection and enhancement of signals and coded information, as well as of energy, may be associated with time or with phase or with both.Thus, investigating the combined effects of time-delay and phase shift on dual-frequency-driven bubble oscillators in the context of vibrational resonance (VR) is essential for understanding a wide range of industrial applications of the physics of bubbles which, in many instances, are time-delayed [35,36,60,67].Many acoustic cavitation effects require that the complex behaviours, often manifesting as chaotic dynamics, are suppressed to ensure optimization of the bubble processes.Such applications are found in, for example, material science where bubbles can exert strong forces able to lift extraneous dirt particles from surfaces.Research on bubble dynamics is also motivated in part by the need to improve the efficiency of sonochemical reactors as well as for increasing the maximum photocurrent by means of dual-harmonics [38].In addition, dual-frequency-driven bubbles can be used in synthesis with nanoparticles [57], sonoluminescence [45], measurement of bubble density in water [62,78,79], measurement of fluid pressure [76], measurement of bubble size distributions [81], measurement of dynamic variations of bubble radii during oscillations [27,70], and in ultrasound to improve the accuracy and efficiency of biomedical diagnosis, especially in tumour therapy [3,32,77,89,[97][98][99][100].
The action of periodic driving forces on nonlinear systems can, in general, lead to the emergence of special types of nonlinear phenomena such as nonlinear resonances, which has been a longstanding research problem due to their occurrence in nearly every field [5,65,84].Among the various forms of nonlinear resonance occurring in mechanical, electrical and biological systems, stochastic and vibrational resonance phenomena are of particular interest and importance.They are strongly connected because dual-frequency forcing is required in each case for their occurrence [9,47,85,86].Stochastic resonance is induced by noise driving [24,28]; while vibrational resonance is a deterministic phenomenon that occurs when there is cooperation between lowfrequency and high-frequency harmonic forces driving a nonlinear system [1,12,49,65].During the last two decades, following its first observation by Landa and McClintock [49], VR has gained continuously increasing research attention due to its potential industrial and biomedical applications in a wide range of contexts, including bistable systems [2,13,21,22], multistable systems [12,64], systems with various forms of potential structures [1,42], linearly damped systems [5,65], electrical circuits and time-delayed systems [20,40,90,91], as well as plasma models [50,71,72].Early studies of VR focused on the parameters of the high-frequency component of the dual-frequency force.Recently however, attention has also been given to other features and parameters of the system including e.g.nonlinear damping/dissipation [50,71,72], the effect of time-delay [17,88] and fractional derivatives [23,33], antiresonance [74], the depth of a potential and the location of its minimum [11], the effects of potential roughness [50] and potential deformation [87], as well as variation in mass [80].More importantly, VR has been realized experimentally in vertical cavity surface emitting lasers and optical systems [12][13][14][15].These investigations have enriched our understanding of the VR phenomenon, its mechanism, and the roles played by relevant system parameters in inducing or promoting VR, as well as in suggesting potential real-life applications in areas, such as ecological systems [41], liquid crystals [30], neural dynamics [8], energy detectors [69] and energy harvesting from mechanical vibrations [18,19], electronic circuits and logic gates [44,56,82,83], a nano-electromechanical resonator [16], and a thermo-optic optomechanical nanocavity [55].
Studies of vibrational resonance in bubble oscillators driven by amplitudemodulated acoustic fields [59] are relatively new, although a diversity of complicated resonance patterns due to variations of amplitude and frequency, such as combination and simultaneous resonances, were reported a little earlier, involving the acoustic cavitation of bubbles under dual-frequency acoustic fields [98].In neither of these works [59,98] was consideration given to the effect of time-delay on the bubble oscillation.Furthermore, the bubble oscillator model in [59] did not consider the impact of a phase shift between the acoustic waves.To the best of our knowledge, the combined influences of time-delay and phase shift on VR have not hitherto been examined or reported.Investigating such effects is important because it promises to illuminate the industrial applications of cavitation effects beyond the already well-known applications of delayed dynamical systems in other fields such as in lasers, telecommunication devices, and in uncovering certain mechanisms of brain dynamics [20,40,90,92].For instance, timedelayed VR in bubbles oscillating in liquids within appropriate parameter regimes could enable control of the time interval between the sub-micron bubble growth and collapse, which could be exploited to enhance the efficiency of ultrasonic cleaning.Our analysis of VR is based on a time-delayed nonlinear Rayleigh-Plesset bubble model driven by two acoustic waves, with a large difference in their frequencies, and with a phase shift between the two driving acoustic waves.The effect of the phase shift on the dynamics of bubble oscillations should not be overlooked.It is well known that phase shift can impact on the acoustic pressure distribution, the energy dissipation, and the production of radicals; and it can also be employed to eliminate chaotic oscillations of the bubbles in connection with chaos control techniques based on periodic perturbations [89,94,95].Note that the application of acoustic cavitation in cleaning processes involves the activation of sub-micron bubbles: dirt particles are lifted from surfaces through the action of strong forces exerted by the bubbles [7,10].
Motivated by the above considerations, as well as by the numerous potential medical and industrial applications of time-delay and phase shift on bubble dynamics, the present paper focuses on the effects of constant time-delay, as well as of phase shift, on the vibrational dynamics of bubble oscillations.Its highlights include observation of a periodic variation in the response amplitude of the bubble oscillator with respect to the phase shift; and enhancement of the VR peak with variation of the phase shift and time-delay, in a cooperative fashion.The latter can be exploited, either to drastically suppress, or to modulate the resonance peaks, thereby controlling the resonances.The rest of the paper is organized as follows: In Section 2, the bubble oscillator model under investigation is presented and described.Section 3 discusses both the theoretical and numerical results.The paper is concluded in Section 4.

Model Description
We now consider small amplitude oscillations of a spherical gas bubble in an incompressible liquid.The Rayleigh-Plesset equation [68] can be used to derive the equation of motion of the bubble, taking into consideration relevant parameters of the surrounding liquid.The properties of the liquid are appropriately specified in formulating the model equations; while those of the spherical gas bubble are variables such as its size and shape, including the instantaneous bubble radius R, whose variation with time can be determined.Other parameters of interest are the bubble's equilibrium radius R 0 , the external pressure in the liquid P ext , and the pressure inside the bubble P in .In addition, we denote the polytropic exponent of the gas in the bubble as κ; while ρ, µ l , µ th , and σ are the density, the liquid and "thermal" viscosity, and the surface tension of the liquid, respectively.
The Rayleigh-Plesset equation [29,52,58,63,68] for bubble oscillations in incompressible liquids can be expressed as; where and The periodic function ǫ cos ωt is considered to be a weak acoustic driving force.P in in Eqn.
(2) is given as, with Substituting Eqns.( 2)-( 5) into Eqn.(1), and redefining the variable R as R = R 0 (1 + x), x being a dimensionless quantity proportional to the bubble radius R, the Rayleigh-Plesset equation (1) in terms of the system's potential then becomes [59] ẍ with and the time-dependent potential V (x, t) given as The dimensionless parameters in Eqns.( 6) and (7) are defined as; From Eqn. (7), β, given by Eqn.(8) is the square of the natural resonant frequency of the system.This closely resembles the results reported by Ida [39] and Zhang [96].In the literature, √ β is also known as the partial or linear natural frequency of the bubble [39,59,96,99].
For acoustically driven bubble oscillators, it is reasonable to assume that the displacement of the bubble is very small (x ≪ 1) because the amplitude of the driving force is in practice usually weak (ǫ ≪ 0.1) [39,59,97].Setting η = 3  2 , and assuming that β and γ are much greater than δ, with ǫ ≪ 0.1, the bubble oscillation can be reduced to with the potential written as, In the absence of the external acoustic fields and time-delay, the system potential given by Eqn. ( 10) admits two types of potential structure, dependent on the values of β, γ and λ.For instance, with β = 145, λ = 14.5 and γ = 32.9, the potential is a single-wellsingle-hump potential.However, when β = 14.5, λ = 14.5 and γ = −32.9, the potential is an asymmetric double-well-double-hump potential.Figure 1a(i) depicts a single-wellsingle-hump potential structure with β = 145, λ = 14.5 and γ = 32.9, while Fig. 1b(i) shows an asymmetric double-well-double-hump potential for β = 14.5, λ = 14.5 and γ = −32.9.For further details, see our recent paper [59], where the potential structure and its connections to the stability of the system were comprehensively discussed.
Throughout the present paper, we employ parameters for which the potential is of single-well form in order to illustrate the impacts of time-delay and phase shift.
If the amplitude of the acoustic force, δǫ cos ωt in Eqn. ( 9) is such that f = |δǫ|, we can then write Eqn.(9) as Eqn. (11) describes the dimensionless nonlinear Rayleigh-Plesset bubble oscillator, with β being the natural oscillation frequency of the system.α 0 , α 1 , α 2 , α 3 and α 4 are the linear and nonlinear damping parameters arising from the thermal and liquid viscosity of the liquid in which the bubble oscillates [96].
To investigate the occurrence of VR in the system given by Eqn.(11), the system is perturbed with a high-frequency acoustic force g cos(Ωt + φ), with Ω being the high-frequency component of the perturbation and φ a phase shift between the two commensurate acoustic driving forces.In acoustic cavitation applications, or in fields where the cavitation effect is very pronounced, suppressing the chaotic oscillations in bubble dynamics is of primary importance in controlling the bubble response and increasing its predictability.In this regard, the multiple-frequency excitation approach has numerous advantages.For instance, in sonoluminescence and sonochemistry, it has been employed in eliminating standing waves, reducing the cavitation threshold, and improving the chemical reaction efficiency, as well as for the effective control of cavitation processes [98,99].
We now focus on the effects of time-delay and phase shift on the VR phenomenon.Accordingly, a constant time-delay parameter (ζ) is introduced into Eqn.(11) in addition to high-frequency acoustic signal g cos(Ωt + φ), so that the dynamics becomes Bubble oscillations in a liquid occur at a natural oscillation frequency that depends on the difference between the average bubble radius and the equilibrium radius corresponding to that pressure.Following a pressure change, the equilibrium radius is approached with an effective time lag.We take account of this lag by introduction of the constant time-delay parameter (ζ) associated with the linear natural frequency term of the oscillating bubble.Because pressure variations could lead to additional resonance, variations in the time-delay could also shed more light on the bubble behavior.This is why we have chosen to introduce ζ only in the linear term x in Eqn. ( 12) which is connected directly to the parameter defining the natural resonant frequency of the bubble.Indeed, it turns out that the delay parameter ζ enhances significantly the bubble's response to dual-frequency acoustic driving fields.On the one hand, the phase shift can be introduced and adjusted appropriately in accordance with phase control techniques for nonlinear systems employing bi-periodic signals [89,94,95].In bubble dynamics, it has previously been shown that the growth or collapse of bubbles is dependent on the phase difference between the sound fields [31,34].In the following analysis and discussion, we show that the response amplitude of the bubble system at the low frequency varies periodically with respect to the phase shift when the phase shift consists of an even number of periods, and can be optimized to enhance the system's response in the appropriate parameter space of the high-frequency driving force.

Results and Discussions
We begin by examining the structure of the effective potential V (χ) which has been derived in the Appendix and given by Eqn.(39).V (χ) is shown in Figs.1a(ii) and b(ii) for the corresponding parameter values of Fig. 1a(i) and b(i) respectively.For the five values of the amplitude of the fast acoustic wave, g = 20, 150, 250, 330 and 400 considered, the effective potential mimic the systems potential with g = 20 as depicted in Figs.1a(ii) and b(ii).As expected from Eqn. (39), C 4 and C 5 reduces to β and γ, respectively, for small values of g.The system's effective potential illustrated in Figs.1a(ii) and b(ii) for higher values of g indicate that decrease in the amplitude of high frequency acoustic wave, increases the degree of skewness of the system's potential.It is also obvious from the figures (1a(ii) and b(ii)), that further increment in the value of g would energize the particles of the bubble, thereby enabling them to overcome the potential barrier created in the neighborhood of x = 0 and x < 0.More details on the structure of the effective potential and the relationships between the potential barriers (humps) and the energy of the bubble particles can be found in Ref. [59].
Next, we examine both the analytically and numerically computed results for the response amplitude.In practice, the analytical response amplitude of the system can be estimated from the values of χ obtained by direct integration of Eqn.(38) (see Appendix ).This can be done using any standard numerical integration scheme, such as the fourth order Runge-Kutta scheme which we have employed here.This approach is particularly effective for strongly nonlinear systems.Alternatively, the analytical response amplitude Q can also be obtained by considering the system's vibrations due to a deviation Y = χ−χ * of the slow motion χ from the equilibrium point χ * .Due to the complexity of the bubble dynamics, and for convenience in examining the influences of the phase and time-delay on the occurrence of VR in the system described by Eqn.(12), the analytical response amplitude was computed from Eqn. (38).The numerical results, on the other hand, were achieved by first expressing the bubble oscillator Eqn.(12) as the following set of two coupled first-order delayed differential equations (DDEs) written as Eqn. (12) was then numerically integrated and the response amplitude Q and phase θ, respectively, were computed from the expressions: and The response of the bubble oscillator given by is 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 follows; In Eqn. ( 16), T = 2π ω is the oscillation period of the low frequency input signal with n = 1, 2, 3 . . .being the number of complete oscillations.A zero initial condition was used for the various computations, and the other fixed parameter values were α 0 = 0.08, α 1 = 0.12, α 2 = 0.16, α 3 = 0.08 and α 4 = 0.04.In addition, the following parameter values were fixed while computing the resonance diagrams; ω = 5, Ω = 20ω, β = 145, γ = 32.9,λ = 14.5, f = 0.01 and g = 5.In all the plots, the continuous curves represent the numerically-computed Q from Eqn. (13) using Eqn.(14), while the analytically calculated Q from Eqns. ( 38) and ( 14), are indicated by dashed lines with marker points.

Phase shift effect on vibrational resonance
We begin our numerical investigation of the VR phenomenon in the delayed bubble oscillator by considering the effect of the phase shift between the two acoustic fields on the response of the system, which is one of the main foci of this paper.It has been reported that bubble growth or collapse depends on the phase difference between acoustic waves [31,34], and that the phase shift can be adjusted to control the chaotic dynamics of acoustically driven bubbles, thereby giving rise to regimes where their chaotic behavior could be stabilized [4,37].In the numerical simulation results and discussions that follow, we show that the phase shift can impact significantly on the occurrence of VR.This result has not been reported previously.Due to the complexity of the bubble dynamics and for convenience, the influence of the phase and the timedelay on the occurrence of VR in the system described by Eqn.(12) was examined by integrating Eqn.(38) to obtain values of x(t) -from whence the theoretical response amplitude Q was computed.In both the numerical and theoretical calculations, the system response Q was calculated from Eqn. (14).Fig. 2(a) demonstrates excellent agreement between theoretical prediction and the numerical simulation results.However, the observed features of the bubble's response to phase shift show significant departures from earlier reports on VR mechanisms [43,72,90,93].In Fig. 2(a), it is clearly evident that the response amplitude exhibits a periodic pattern in multiples of 2π, indicating a significant response to variations in the phase shift.In addition, it is obvious from the analytical expression given by Eqn.(38), that the phase parameter φ only appears in C 7 = g sin(φ), implying that the parameter C 7 is also 2π-periodic.This is illustrated in Fig. 2(a) by broken-lines with green hexagram markers.By comparing the numerically and analytically computed system response amplitude Q φm with the variation of C 7 , we can validate the 2π-periodicity of Q as shown in Figure 2(a), assuming that φ 0 = 0 and φ m = φ 0 +2mπ, with m = 1, 2, 3, 4, . ... The continuous lines in Fig. 2(a) shows the repeated periodic response amplitudes and Q φ 5 corresponding to φ = 0, 2π, 4π, 6π, 8π and 10π, respectively.The repeated pattern is further confirmed in Figure 2(b) where we have plotted the dependence of Q on the driving amplitude g for different values of φ.
At this point, it is worth mentioning that in most of the oscillators investigated for VR, the impact of phase shift has been neglected or overlooked because they do not exert a meaningful influence on the occurrence of VR.In the bubble oscillator, however, the scenario is different.Maximum peaks occurred and were observed for phase shifts corresponding to nπ 2 , n being odd integers.Although the system can respond appreciably to the dual-frequency acoustic fields for all values of n, the periodic nature of the resonance curves shows that only odd periods of nπ 2 give rise to a maximal response.This implies that, in situations where optimization of the VR phenomenon is essential in the bubble oscillator, odd values of n should be chosen.This complements the well established fact that, in bubble systems, the driving frequency of the external acoustic wave matches the bubble's natural frequency when the phase shift between the bubble's pulsation and an external sound wave is π 2 [39,54].Correspondingly, Q = 1 √ S in the linearized equation of the system given by Eqn.(42)(see Appendix ), where S = (ω 2 r −ω 2 ) 2 +C 2 0 ω 2 , with C 0 being the linear dissipation parameter.Therefore, Q attains its maximum provided S is minimum, implying that (ω 2 r − ω 2 ) 2 + C 2 0 ω 2 is minimum.Thus, ω r = ω at resonance, and from Eqn. (40), C 0 can only approach zero (C 0 → 0).This implies that the contributions of the thermal and liquid viscosity are negligible compared to the parameters of the fast acoustic fields, but not exactly zero, hence, validating the condition for the occurrence of VR.Otherwise, VR would not be observed because the condition would reduce Eqn.(42) to a linear differential equation.
To complete this section, we reveal in Fig. 3(a) the dependence of the response Q on both the phase shift φ and the amplitude g of the high-frequency component of the acoustic waves.We computed this in the range (φ, g) ∈ [(0, 10π), (0.0, 600)], with Ω = 20ω; while other parameters were fixed as before.The red areas in Fig. 3(a) indicate regimes of strong enhancement corresponding to phase differences of nπ 2 , with g = 400.The system resonates for g > 100, with a peak appearing at g = 400, and approaches a nonzero limiting value for g > 500.This is obvious in Fig. 3(b) where we have zoomed a small region of Fig. 3(a) to capture the hidden features within a short range of the phase difference.The plot expands a suitable parameter space with a good choice of φ and g values, yielding optimal enhancement of the system's response.This could have applications in sonochemistry, where selection of the parameter values of the high-frequency component of the acoustic field is desirable to optimise performance.14), while the analytically calculated Q were obtained using Eqn.(38) and Eqn. ( 14), are indicated by marker points with broken lines.

Time-delay Vibrational resonance
Figure 4(a) plots the response amplitude Q versus the amplitude of the high frequency driving field g for different values of time-delay (ζ = 0.10, 0.20, 0.60, 1.50 and 3.0), with φ = 1.17π keeping other parameters fixed.Here, the response amplitude Q of the delayed bubble system is greatly enhanced as indicated by the first peak of the curves corresponding to ζ = 0.10.In addition, the system's response is greatly improved, attaining a peak value of Q ≥ 70, spanning through a wide range of amplitudes of the fast acoustic field [g ∈ (0, 500)], indicated by the almost saturation window, and quickly dropping to zero when g > 500.However, increasing the value of the time-delay ζ drastically suppresses Q with corresponding compression in the range of g values for which enhancement occurs.At ζ = 1.5, the shape of the resonance peak changes and shifts towards lower g values.In addition, Fig. 4(b) shows the dependence of the system's response on the phase shift for six different values of the time-delay (ζ = 0.50, 0.60, 0.70, 0.80, 0.904, and 1.00) in the short range [φ ∈ (0, 3π)].In this case, the system exhibits a different resonance pattern for smaller values of g -attaining a plateau of semi-circular shape, implying that the resonance peaks as well as their shapes can be controlled effectively by adjustment of the time-delay, ζ.As clearly illustrated in Fig. 4(b), the periodically varying resonance shown in Figs.2(a) and 3(a), is optimally enhanced for even periods corresponding to nπ 2 .Moreover, increasing the value of timedelay suppresses the system's response, whereas the response is enhanced at the odd periods of nπ 2 , with fixed Q = 9.5 even when ζ varies.Figure 5(a) displays the response amplitude Q versus the amplitude of the high-frequency acoustic field for different values of the time-delay.Here, the response amplitude is significantly enhanced at lower values of ζ, and experiences gradual suppression with increase in ζ.This effect is similar to that shown in Fig. 4(b) for fixed value of g, showing that the system's response at the lower frequency ω can be controlled by a gradual variation of the time-delay ζ.
Note that Heckman et al. [36] showed earlier that time-delay in the oscillations of one bubble during the process of signal transmission to another bubble through their surrounding liquid medium (i.e.delayed coupling) shifts the stability of the bubble dynamics from stable equilibrium to an unstable state, provided that the delay time is sufficient.However, in most multistable mechanical systems the dynamics of the system is independent of the time delay.For instance, Rajasekar and Sanjuán [66], showed theoretically that both the dynamics and the response amplitude of the Duffing oscillator are independent of the time-delay parameter, even when coupled.This significant departure from previous results may be attributable to the complex and nonlinear structure of the interaction between bubble oscillators, which we will report in a forthcoming paper.It thus suggests that the bubble oscillator, although exhibiting a plethora of properties that are common to all nonlinear systems, has certain features that are peculiar to itself and therefore requires individual research attention.The occurrence of time-delayed vibrational resonance phenomena in bubbles implies that certain cavitation effects could be controlled by utilizing time-delay.Thus, oscillation  14), while the analytically calculated Q from Eqn. (38) and Eqn.(14), are indicated by marker points with broken lines.control of the bubble dynamics could be employed to adjust the resultant response of the system.Notably, for fixed value of the delay, the response amplitude Q plotted as a function of g is enhanced by increasing φ, as depicted in Fig. 5(b).
Next, we discuss the effect of a variation in amplitude of the high-frequency component of the acoustic field on the dependence of the response amplitude Q on the time-delay parameter.This is shown in Fig. 6(a), for φ = 1.17π keeping other parameters fixed.It is evident that, maximizing the amplitude of the high-frequency component of the acoustic waves suppresses the system's response amplitude, with the occurrence of resonance peaks at lower values of ζ (ζ ≤ 0.5).Remarkably, when g becomes large (typically g ≥ 5.00), multiple VR peaks appear.This is further illustrated in Fig. 6(b) which presents a three-dimensional plot showing the dependence of Q on the amplitude of the fast acoustic field g and the time-delay parameter ζ.The parameter regime in which the multiple resonance peaks occur could be identified.For instance, for g ≥ 500 and ζ ≤ 2.0, two resonance peaks arise.Moreover, for g → 1000 with ζ > 2.0, a third resonance peak gradually emerges.Thus, additional resonance peaks, with considerable enhancement of response amplitude of the output signal of the system at the low frequency ω occur only for large values of g and low values of the delay time ζ.

Concluding Remarks
In this paper, we have examined the impacts on VR phenomena of having a constant time-delay in the oscillations of a dual-frequency-driven spherical gas bubble, under conditions where there exists a phase shift between the low and high frequency acoustic driving forces.We presented theoretical results based on the method of direct separation of the fast and slow motions, complemented with numerical simulations.We analyzed VR, focusing on the influence of time-delay and phase shift on its occurrence.In addition to enhancement of the VR peak by variation of the phase shift, we find that the response amplitude of the system at low frequency varies periodically with respect to the phase shift when the phase shift consists of even-periods, implying an optimization of the system's response.Furthermore, time-delay also plays a significant role in the response enhancement and can be exploited either to suppress drastically or to modulate the resonance peaks, thereby controlling the resonances.Further analysis reveals that the cooperation between the time-delay and the amplitude of the high-frequency component of the acoustic waves can induce multiple resonances.
These results could be exploited to control ultrasonic cleaning by varying the time-delay parameter in the presence of phase shifted dual-frequency acoustic waves.Moreover, the efficiency of the cleaning process might be significantly enhanced by an appropriate choice of constant phase shift.Moreover, improved accuracy in ultrasound biomedical diagnosis and tumour therapies could be achieved through time-delayed vibrational resonance.With the cooperation of the amplitude of the high-frequency acoustic field, the delivery of reagents transported within the bubble that is targeted at combating the tumors can be optimized.
Looking to the future, we note that the results reported herein have opened up several interesting questions for further investigation.These include: how time-delay and phase-shift in VR impacts on the Blake's threshold; whether there is a maximum delay that the bubble can accommodate in the course of radial expansion or contraction.If the time delay on x becomes relatively large, is there threshold value beyond which the bubble will collapse?Investigations of these questions and related matters will be developed in detail elsewhere.and m ψ = F ( χ, χ, t) − F ( χ, χ, ψ, ψ, t) + Φ( χ + ψ, χ + ψ, t, τ ) ( − Φ( χ + ψ, χ + ψ, t, τ ) , where is the function of the system's variables from which one can compute the response amplitudes of the system.The first of these equations, namely Eqn.20 arises from the substitution of Eqn.(4.1) into Eqn.(17), averaging both sides with respect to the fast time τ according to Hypothesis 4.2.The second equation is then obtained by subtracting Eqn.(20) from Eqn. (17).
This technique, known as the method of direct separation of motion has been explicitly described in Refs.[6,65], and employed successfully to a broad range of theoretical analysis of VR [71-73, 75, 87].In general, one will often have to seek the approximate solution for the fast component ψ * , in the form of a small number of harmonics.If ψ is considered to be small compared to χ, then F and Φ may be linearised with respect to ψ (and possibly χ) to find a solution.
Throughout our analysis of VR in this paper, it is assumed that the fast motion ψ * is asymptotically stable so that the potential V (x) is well-defined over a certain range of initial conditions of the fast variables ψ.If there are several stable fast motions, then, the potential V (x) will depend on which motion is considered.In the analysis, Eqn.23 is the main equation that gives the vibrational mechanics of the system [6].It is an equation for the slow dynamics with an effective potential due to the fast dynamics.Now, we seek the solution of the equation by splitting the motion of the bubble system into fast and slow dynamics according to Hypothesis 4.1 and 4.2, thereby enabling us to obtain two integro-differential equations for the time-delayed bubble oscillator.One of them describes the slow motion of the system whose response can be modulated by the variation of the parameters of the high-frequency acoustic field.The response of the system, denoted by Q and defined as the ratio of the amplitude A L of the bubble oscillator to the frequency f is obtained by solving the equation for the slow motion.
Using the time-delayed bubble oscillator given by Eqn.(12), we can express Eqn.(18) as where ζ is a constant time-delay with Here, χ(t) is periodic with period T = 2π ω while ψ is periodic in the fast time τ = Ωt with period 2π.Its average value with respect to fast time τ is given by Eqn (19).
The goal is to derive a system of two-coupled integro-differential equations for the time-delayed variables χ and ψ from the main equation of the system (12), focusing on the slow components of the motion which is of principal interest here.Substituting Eqns.( 18) and ( 25) into Eqn.( 12), we have Averaging both sides of Eqn. ( 27) over the period of fast time [0, 2π Ω ] and substituting the mean value of ψ as expressed in the Eqn.(19), considering that ψ is a rapidly oscillating periodic function of the fast time, we can then re-express Eqn.(27) as: We can further simplify Eqn.(28) using the expression for the mean values in Eqn.(19).To proceed, we make the following remark and definition: Remarks 4.1 Since the phase shift is time-independent so that the biharmonic acoustic waves are out-of-phase, the average value of the fast signal over the period is sin φ sin Ωt, and approaches zero as φ → 0 when the biharmonic forces are in-phase.This allows us to make the following definition: If however the biharmonic forces are in-phase (φ = 0), then, the hypothesis (Eqn.(4.2)) does not hold.The approximate expression in Eqn.(30) clearly shows that the contributions of φ to the system's response is non-trivial.Using Eqns.( 29) and (30) above, Eqn.(28) reduces to Eqn. ( 31) is the first of the two coupled equations for the variable χ, for the slow oscillation of the bubble.It can be used in computing the theoretical response, Q of the system at the lower frequency ω.The equation of fast oscillatory motion can be obtained by subtracting Eqn. ( 31 Eqns. ( 31) and ( 32) are the two integro-differential equations of motion, namely, for the slow motion in χ and for the fast motion in ψ, respectively.Here, our interest lies in the equation of motion of the slow dynamics of the system Eqn.(31) which could be modulated appropriately by varying the parameters of the fast acoustic field to achieve VR.Thus, we assumed that the component ψ is much more faster than the slow component χ, while the components of the slow dynamics (i.e χ and χ) are considered frozen in Eqn.(32).Application of the inertial approximation approach, and assuming that ψ ≫ ψ ≫ ψ, enables Eqn.(32) to be approximated as ψ = g cos(Ωt + φ), (33) which has a solution such that, Eqn. (37) can be written in an elegant form as; so that the effective potential of the time-delayed system becomes where, and C 4 is the effective natural frequency defined in Eqn (40).
where ω r = √ C 4 is the natural resonance frequency.Eqn. ( 42) has a steady-state solution, Y ζ = A L cos(ωt + θ), where the response amplitude A L can be expressed as and the phase The response amplitude Q can be calculated from

Figure 1
Figure 1: a(i) The single-well-single-hump potential structure of the bubble oscillator with β = 145, γ = 32.9,λ = 14.5; a(ii) effective potential corresponding to slow motion of the system for five different values of the amplitude of the fast acoustic field (g = 20, 150, 250, 330 and 400), with Ω = 20ω and other parameters fixed as in a(i); b(i) the asymmetric double well potential of the dimensionless bubble oscillator with β = 14.5, γ = −32.9,λ = 14.5 ; b(ii) The effective potential structure of the bubble oscillator for five different values of the amplitude of the fast acoustic field (g = 20, 150, 250, 330 and 400), with Ω = 10ω and other parameters fixed as in b(i).

Figure 2 :
Figure 2: (a) Response Q against φ the phase shift of the fast acoustic force with ω = 5, Ω = 20ω, β = 145, γ = 32.9,λ = 14.5, f = 0.01, g = 5 and ζ = 0.5 with other parameters fixed.Continuous curves represent the numerically-computed Q from Eqn.(13) using Eqn.(14), while the analytically calculated Q from Eqn. (38) and Eqn.(14), are indicated by marker points with broken lines; (b) Repeated pattern of the response Q as a function of the amplitude g of the fast acoustic field with other parameters fixed as in (a).

Figure 3 :
Figure 3: (a) Surface plots of Q against φ and g with ω = 5, Ω = 20ω, β = 145, γ = 32.9,λ = 14.5, f = 0.01 and ζ = 0.5 showing that the system resonates for g > 100, with a peak located at g = 400, and approaches a nonzero limiting value for g > 500.The red areas indicate regimes of strong enhancement corresponding to phase differences of nπ 2 .(b) is the zoom of a restricted region of (a) capturing the important and hidden features of the resonance diagram (i.e. the peak and the depth) within a short range of the phase difference φ.

Figure 4 :
Figure 4: (a) Dependence of Q on g for different values of time-delay (ζ = 0.10, 0.20, 0.60, 1.50 and 3.0), showing significant enhancement of the response with the attainment of a peak value of Q ≥ 70, over a wide range of amplitudes of the fast acoustic field for φ = 1.17π while keeping other parameters values fixed.(b) The dependence of Q on φ for different values of the time-delay (ζ = 0.50,0.60,0.70, 0.80, 0.904, and 1.00) at g = 5, and other parameters fixed and illustrating different resonance pattern for smaller values of g, with the attainment of a plateau of semi-circular shape -indicating the control impact of ζ.The continuous curves represent the numerically-computed Q from Eqn. (13) using Eqn.(14), while the analytically calculated Q were obtained using Eqn.(38) and Eqn.(14), are indicated by marker points with broken lines.

Figure 5 :
Figure 5: (a) The response Q against the amplitude of the fast acoustic field with varying time-delay ζ, and φ = 3 2 π and other parameters fixed showing significant enhancement at lower values of ζ, and gradual suppression with increase in ζ; (b) The response Q against the amplitude of the fast acoustic field with varying φ, and ζ = 1.2.Continuous curves represent the numerically-computed Q from Eqn.(13) using Eqn.(14), while the analytically calculated Q from Eqn. (38) and Eqn.(14), are indicated by marker points with broken lines.

Figure 6 :
Figure 6: (a) The dependence of Q on the time delay parameter ζ with φ = 1.17π showing that maximizing the amplitude of the high-frequency component of the acoustic waves suppresses the system's response amplitude, with peaks occurring at lower values of ζ (ζ ≤ 0.5).Continuous curves represent the numerically-computed Q from Eqn.(13) using Eqn.(14), while the analytically calculated Q from Eqn. (38) and Eqn.(14), are indicated by marker points with broken lines; (b) Three-dimensional plot showing the dependence of the response amplitude Q on g and the delay parameter ζ, with the appearance of multiple VR peaks as g becomes progressively large (typically g ≥ 5.0), for φ = 1.17π and Ω = 10ω.Two resonance peaks are evident for g ≥ 500 and ζ ≥ 2.0.