2 Asymptotic form of the fluid dynamical equations

We consider a gas bubble in an infinite incompressible liquid and driven by an acoustic field. For the sake of simplicity, we assume a uniform incident acoustic pressure field and neglect the translational motion of the bubble. The bubble surface is described by the radial coordinate

where the small parameter $\unicode[STIX]{x1D716}$ stands for the ratio of shape modes amplitude to bubble radius ( $\unicode[STIX]{x1D716}\ll 1$ ). Following the mathematical formulation of Shaw (Reference Shaw2006), the dynamical equations accurate to second order and governing the spherical mode $R(t)$ and the shape modes $a_{n}(t)$ ( $n\geqslant 2$ ) are written as

and

where $R_{0}$ is the static bubble radius, $\unicode[STIX]{x1D70E}$ is the surface tension, $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D708}$ denotes the density and kinematic viscosity of the liquid, respectively, $\unicode[STIX]{x1D6FE}$ is the polytropic exponent of the gas, $p_{b0}$ is the equilibrium gas pressure inside the bubble defined as

where $p_{0}$ is the static liquid pressure, and $p_{\infty }$ is the instantaneous liquid pressure written in the isotropic acoustic forcing assumption as

where $p_{a}$ is the acoustic pressure amplitude and $\unicode[STIX]{x1D714}$ is the angular frequency. We must emphasise that viscous dissipation in the liquid has been added in a ad hoc manner in (2.2) and (2.3) by introducing the classical viscous terms to the lowest order. The viscous effects in the shape modes dynamical equation are taken into account through the well-known functions obtained in the boundary layer approximation

where $\unicode[STIX]{x1D6FF}=\sqrt{\unicode[STIX]{x1D708}/\unicode[STIX]{x1D714}}$ is the viscous boundary layer thickness (Brenner et al. Reference Brenner, Lohse and Dupont1995). Note that a rigorous mathematical treatment of the liquid viscosity would lead to nonlinear viscous terms at the higher orders in (2.2) and (2.3) (Shaw Reference Shaw2009). However, as will be specified in the following, we will consider viscosity as a small perturbation so that such higher-order terms should be neglected in the further asymptotic expansion of (2.2) and (2.3) anyway. Finally, we also specify that the quantities denoted $I_{a}$ , $M_{a}$ , $M_{b}$ , $M_{c}$ , $M_{d}$ and $G_{d}$ in (2.3) involve integrals of triple products of Legendre polynomials and their derivatives. Their definitions can be found in the work of Shaw (Reference Shaw2006) and those retained in the later derivations are recalled in appendix A. These integrals can be also expressed by means of the Clebsch–Gordan coefficients (Tsamopoulos & Brown Reference Tsamopoulos and Brown1983; Doinikov Reference Doinikov2004). It is worth noting that (2.2) and (2.3) are equivalent to those derived by Doinikov (Reference Doinikov2004, equations (2.23) and (2.25)) under the assumptions of isotropic pressure forcing and negligible translational motion.

Let us now introduce new non-dimensional coordinates $x(t)$ and $s_{n}(t)$ defined as

and the new time scale $\unicode[STIX]{x1D70F}$ defined as

In the following, we assume that $x$ , $p_{a}$ and the squared ratio $(\unicode[STIX]{x1D6FF}/R_{0})^{2}$ measuring the importance of viscous damping are all small quantities of the order of $\unicode[STIX]{x1D716}$ . Thus performing a power series expansion of $x$ and neglecting terms of third and higher orders in $\unicode[STIX]{x1D716}$ , equations (2.2) and (2.3) become

and

where the dot now denotes time derivative with respect to the variable $\unicode[STIX]{x1D70F}$ . Several non-dimensional parameters appear in (2.9) and (2.10) which are defined as

and the viscous dissipation functions are expressed at the lowest order following Discarding the quadratic terms involving shape modes in (2.9) leads to the equation derived by Prosperetti (Reference Prosperetti1974, equation (4)) reduced to second order, except for the difference in the chosen time scale. In these equations, we recognise the non-dimensional forms of the resonant frequency of the spherical oscillations $\unicode[STIX]{x1D714}_{0}$ and the resonant frequency of the shape modes $\unicode[STIX]{x1D714}_{0,n}$ (Lamb Reference Lamb1932).

3 First-order solutions

We now look for solutions of (2.9) and (2.10) in the form of asymptotic expansions

Substituting (3.1) and (3.2) into (2.9) and (2.10) and equating the terms of same order in $\unicode[STIX]{x1D716}$ leads to the set of equations

To ensure the validity of the asymptotic expansions (3.1) and (3.2), suitable mathematical forms of the solutions at the lowest order must be considered depending on the resonance conditions of the system. For instance, considering a shape mode $n$ close to a parametric resonance so that $\unicode[STIX]{x1D714}_{0,n}\simeq k$ , where $k$ is an integer, a suitable solution of (3.3b ) must be chosen as

In the following, we limit our analysis to the first parametric resonance of the shape mode $n$ , thus setting $k=1$ . On the other hand, we consider non-resonant spherical oscillations of the bubble, sufficiently far away from the fundamental or harmonic resonances, so that $\unicode[STIX]{x1D714}_{0}\neq 2k$ . In the non-resonant region, the solution of (3.3a ) is written as

by implicitly assuming that the steady state of the spherical oscillations is already established at the initial time. Different definitions of $\unicode[STIX]{x1D703}_{x}$ should be considered in the vicinity of resonances (Prosperetti Reference Prosperetti1974).

Substituting (3.5) and (3.4) into (3.3c ) and (3.3d ) produces secular terms that must be cancelled to ensure uniformly valid solutions. Depending on their respective phases $\unicode[STIX]{x1D703}_{i}$ and $\unicode[STIX]{x1D703}_{j}$ , the product of two shape modes can also produce secular terms as soon as $\unicode[STIX]{x1D703}_{i}-\unicode[STIX]{x1D703}_{j}\sim \unicode[STIX]{x1D703}_{n}$ or $\unicode[STIX]{x1D703}_{i}+\unicode[STIX]{x1D703}_{j}\sim \unicode[STIX]{x1D703}_{n}$ . Therefore, closing the formulation requires us to consider the other modes of integers $m\neq n$ and their respective dynamics. All these secondary modes are considered as far from their own parametric resonances, so that $\unicode[STIX]{x1D714}_{0,m}\neq k$ . Contrary to the parametrically forced mode $n$ , the modes $m\neq n$ are thus assumed to oscillate at the frequency of the acoustic forcing at the lowest order, and we write

Naturally, a complete solution for a non-resonant mode $m$ should include a freely oscillating part, however this part would exponentially decrease due to viscous damping anyway, so (3.6) is sufficient to describe the steady state of the oscillations. The solution at first order for $m\neq n$ must also be solution of (3.3d ).

Following the Krylov–Bogoliubov asymptotic method (Nayfeh Reference Nayfeh2008), allowing the amplitudes $X$ , $S_{n}$ , $S_{m}$ and phases $\unicode[STIX]{x1D719}_{x}$ , $\unicode[STIX]{x1D719}_{n}$ , $\unicode[STIX]{x1D719}_{m}$ to be slowly varying functions at the time scale $\unicode[STIX]{x1D70F}$ and cancelling the secular terms appearing in (3.3c ) and (3.3d ) leads to the first-order differential equations

governing the amplitudes and phases of the zero-order solutions, where the coefficients $\unicode[STIX]{x1D6FD}_{n}$ , $\unicode[STIX]{x1D6FE}_{n}$ , $W_{nm}$ and $O_{mn}$ are summarised in appendix A. The solutions of (3.3c ) and (3.3d ) free of secular terms are then given by

where the coefficients $\unicode[STIX]{x1D6FC}_{n}$ , $\unicode[STIX]{x1D6FF}_{n}$ and $J_{a_{nij}}$ are also summarised in appendix A.

Figure 1.

Oscillations of the spherical mode and the seven first shape modes for a bubble of radius $R_{0}=44.8~\unicode[STIX]{x03BC}\text{m}$ driven above the first parametric resonance of the $n=2$ mode ( $\unicode[STIX]{x1D714}/(2\unicode[STIX]{x03C0})=31.25~\text{kHz}$ , $P_{a}=22.25~\text{kPa}$ , $\unicode[STIX]{x1D714}_{0,2}\sim 1$ , corresponding to the point $\triangledown$ in figure 5). Numerical simulations (blue curves) are compared with first-order approximate solutions (orange curves). (a) Presents the results over long times, from transient oscillations to the establishment of the steady-state regime. (b) Focuses on few acoustic periods during the steady state.

Figure 2.

Same as figure 1, for a bubble of radius $R_{0}=66.9~\unicode[STIX]{x03BC}\text{m}$ driven above the first parametric resonance of the $n=3$ mode ( $\unicode[STIX]{x1D714}/(2\unicode[STIX]{x03C0})=31.25~\text{kHz}$ , $P_{a}=12.69~\text{kPa}$ , $\unicode[STIX]{x1D714}_{0,3}\sim 1$ , corresponding to the point ▫ in figure 6).

Figure 3.

Same as figure 1, for a bubble of radius $R_{0}=90.72~\unicode[STIX]{x03BC}\text{m}$ driven above the first parametric resonance of the $n=4$ mode ( $\unicode[STIX]{x1D714}/(2\unicode[STIX]{x03C0})=31.25~\text{kHz}$ , $P_{a}=3.62~\text{kPa}$ , $\unicode[STIX]{x1D714}_{0,4}\sim 0.95$ , corresponding to the point ▿ in figure 7).

Figure 4.

Modal amplitudes including the spherical ( $0$ ) and translational ( $1$ ) vibrational modes, for two bubbles driven at $\unicode[STIX]{x1D714}/(2\unicode[STIX]{x03C0})=10~\text{kHz}$ . (a) $R_{0}=144~\unicode[STIX]{x03BC}\text{m}$ , $P_{a}=13~\text{kPa}$ , $\unicode[STIX]{x1D714}_{0,3}\sim 0.99$ . (b) $R_{0}=189~\unicode[STIX]{x03BC}\text{m}$ , $P_{a}=8~\text{kPa}$ , $\unicode[STIX]{x1D714}_{0,4}\sim 0.99$ . Circles (○) correspond to results from the Krylov–Bogoliubov perturbation method, while diamonds (♢) correspond to amplitudes in steady state obtained from the numerical simulations presented in the paper of Shaw (Reference Shaw2017). For the spherical mode, we represent the difference between the amplitude before ( $-$ ) and after (symbol) the establishment of the surface oscillations. Sustaining shape oscillations generally goes with a decreasing of the amplitude of the spherical oscillations.

The approximate solutions (3.5)–(3.6) and (3.8)–(3.10) are compared to the direct numerical integration of (2.2) and (2.3) in figures 1–3, for the parametric excitation of the first three shape modes $n=2$ , $n=3$ and $n=4$ , respectively. We consider air bubbles of typical radii around $R_{0}\simeq 50{-}90~\unicode[STIX]{x03BC}\text{m}$ in water, driven by an external forcing at frequency $\unicode[STIX]{x1D714}/(2\unicode[STIX]{x03C0})=31.25~\text{kHz}$ , following experimental conditions available in the literature (Guédra et al. Reference Guédra, Cleve, Mauger, Blanc-Benon and Inserra2017). Eight modes are retained in the computation. The initial conditions for all shape modes are set to $s_{n}(0)=s_{m}(0)=10^{-2}$ and ${\dot{s}}_{n}(0)={\dot{s}}_{m}(0)=0$ , while the steady state of the spherical oscillations is already established at initial time as mentioned previously (a simulation of the Rayleigh–Plesset equation has been run prior to the calculation to set the initial conditions for $x$ and ${\dot{x}}$ ). The simulations are then performed over 100 periods, corresponding to a total time of 3.2 ms.

Looking at the general dynamics shown in figures 1–3 first indicates that the equations obtained from the perturbation method reproduce the characteristic exponential growth of the instability and capture the nonlinear saturation of the shape modes. Finally, the steady-state regime with finite-amplitude oscillations is recovered. The unstable shape mode close to its first parametric resonance naturally oscillates mainly at half the driving frequency ( $\unicode[STIX]{x1D714}_{0,n}\simeq 1$ ). The approximate solutions, through the differential set of equations (3.7a–f ), also capture the exponential decay of damped non-resonant shape modes (being roughly of the type $S_{m}\propto \exp [-F_{\unicode[STIX]{x1D708},0}^{(m)}\unicode[STIX]{x1D70F}/2]$ ) as well as the possible excitation of some particular shape modes coming from nonlinear coupling. The presence or absence of shape modes is ruled by the values of $W_{mn}$ and $O_{nm}$ involved in (3.7a–f ), which can be zero for particular $(n,m)$ combinations. Comparing figures 1 to 3 notably illustrates that parametrically forced even modes are only able to excite other even modes (at least at this order of approximation), while odd modes can excite both even and odd modes (Shaw Reference Shaw2006). These results also highlight the important drift effect observed for the $s_{2}$ shape mode in figures 2 and 3 coming from the non-oscillating terms in (3.10). This effect is a typical signature of nonlinear mode coupling and has already been observed both experimentally by Guédra et al. (Reference Guédra, Inserra, Mauger and Gilles2016b ) and numerically by Shaw (Reference Shaw2017).

In figure 4, the approximate solutions are compared to a more complete theory derived by Shaw (Reference Shaw2017) including translational motion and liquid compressibility up to third order. The calculations are made for a driving frequency $\unicode[STIX]{x1D714}/(2\unicode[STIX]{x03C0})=10~\text{kHz}$ , in the vicinity of the first parametric resonance of the shape modes $n=3$ (i) and $n=4$ (ii), referring to the figures 5–6 and 9–10 presented in the paper of Shaw (Reference Shaw2017). The steady-state amplitudes are given as a function of the mode number, up to 16 modes retained for computation. Even with several assumptions made in the perturbation analysis, the general trend is satisfactory. Moreover, these results highlight the decreasing amplitude of the spherical solution (mode $0$ ) due to the triggering of non-spherical oscillations. This specific influence of shape modes on the spherical one will be further commented on in the next section.

Finally, we must also pay attention to the boundary layer approximation used to introduce viscous dissipation in the present model. This assumption should be used with caution as it discards the continuous diffusion of vorticity by viscosity (Prosperetti Reference Prosperetti1977) which leads to increasing the boundary layer on long times. The importance of vorticity effects can be inferred from the estimate of the characteristic time $R_{0}^{2}/\unicode[STIX]{x1D708}$ representative of viscous diffusion. This time ranges from 2.5 to 8 ms for the micrometre-sized bubbles considered in our simulations, which is comparable to the total computation time and cannot be strictly considered as much larger than the time required for the growth and saturation of the surface oscillations. Moreover an estimate of the Reynolds number from the calculation of the characteristic capillary time leads to values as low as $Re\simeq 85$ , which indicates that viscous effects are probably underestimated by the boundary layer approximation. For the bubbles considered in the simulations, it can in fact be expected that stronger viscous damping can participate in the saturation of the shape oscillations to a similar extent as nonlinear mode coupling. Thus a quantitative comparison with experimental data would certainly require a more rigorous mathematical treatment of the viscous terms, notably accounting for the growth of the viscous layer convected by vorticity as considered in the works of Lundgren & Mansour (Reference Lundgren and Mansour1988).

4 Steady-state solutions

In the previous section, the temporal dynamics obtained from numerical integration revealed that nonlinear mode coupling is responsible for a saturation of the instability leading to stable non-spherical oscillations. In order to obtain the corresponding steady-state solutions, the set of differential equations (3.7a–f ) is first put in a more convenient complex form by introducing the complex variables $\bar{X}=X\text{e}^{\text{i}\unicode[STIX]{x1D719}_{x}}$ , $\bar{S_{n}}=S_{n}\text{e}^{\text{i}\unicode[STIX]{x1D719}_{n}}$ and $\bar{S_{m}}=S_{m}\text{e}^{\text{i}\unicode[STIX]{x1D719}_{m}}$ to obtain

where $\ast$ denotes the conjugate, and the complex quantities $\bar{\unicode[STIX]{x1D6E5}}$ are defined by

In the steady-state regime, cancelling the first temporal derivatives in (4.1a ) and (4.1c ) directly provides the complex amplitudes of the spherical and secondary shape modes

Then substituting (4.5) and (4.6) into (4.1b ) leads to the equation

which has the trivial solution $S_{n}=0$ (with undetermined phase $\unicode[STIX]{x1D719}_{n}$ ), and solutions ruled by

where we have introduced the complex quantities

for the conciseness of the derivations.

Figure 5.

Driving amplitude as a function of the frequency $\unicode[STIX]{x1D714}_{0,n}$ , for $n=2$ . The absolute stability threshold computed from (4.13) is given by the thick curve. The conditional stability threshold computed from (4.17) is given by the dashed curve. The other curves are isolines corresponding to fixed values of $S_{n}$ in steady state: $S_{n}=0.1$ (♢), $0.2$ (▿), $0.3$ (▫) and $0.4$ (○). The dotted lines correspond to conditionally stable branches ( $A_{th}^{\prime }<A<A_{th}$ ).

We emphasise that these solutions include the onset threshold for the parametric excitation of the shape mode $n$ in the vicinity of the first resonance. A suitable expression can be simply deduced by discarding the quadratic terms in (4.5) and (4.8) to obtain the transition curves, in terms of external driving amplitude (acoustic pressure):

or in terms of spherical mode driving amplitude:

Note that (4.14) is almost equivalent to that derived by Francescutto & Nabergoj (Reference Francescutto and Nabergoj1978), except for small differences in the viscous damping terms due to the fact that these authors only considered the lowest-order viscous terms in the boundary layer approximation following the derivations of Ceschia & Nabergoj (Reference Ceschia and Nabergoj1978).

When accounting for nonlinear interaction terms, solving (4.8) and (4.9) allows us to determine the steady-state solutions above the threshold. It is worth noting that (4.8) directly provides the value of applied pressure amplitude $A$ to get a given amplitude $S_{n}$ of the parametrically forced shape mode $n$ . This is illustrated in figures 5–7 which give the stability areas of the three first shape modes $n=2$ , $n=3$ and $n=4$ , respectively, in the vicinity of the first parametric resonance. The onset threshold curves computed from (4.13) are given by the thick black curves. In each figure, the driving amplitude curves corresponding to isovalues of the shape mode amplitude $S_{n}$ ( $S_{n}=0.1,0.2,0.3,0.4$ ) and computed from (4.8) are labelled with different symbols.

Figure 6.

Same as figure 5, for $n=3$ .

Figure 7.

Same as figure 5, for $n=4$ .

These results enable us to illustrate an important feature of the steady-state solutions, namely that the amplitude curves for a finite value of shape mode amplitude cross the onset threshold curve at particular frequencies. This has two major consequences on the sustained shape oscillations above parametric threshold. First, it indicates that increasing the acoustic driving amplitude to exceed the onset threshold in the region where curves are crossing leads to a jump of the shape mode amplitude from a zero value to a finite non-zero value. Second, one can notice that this region where such a jump is appearing always corresponds to values $\unicode[STIX]{x1D714}_{0,n}>1$ , namely for bubble static radii below the parametric resonance radius. It means that for a bubble of static radius smaller than the parametric resonance radius, shape modes suddenly appear with a sharp transition at the onset threshold when increasing the acoustic pressure amplitude. Conversely, shape modes appear with a smooth transition at threshold for a bubble of static radius larger than the parametric resonance radius.

Let us now look at the steady-state solutions of (4.8) which can be expressed as:

where $h_{i}=|\bar{h}_{i}|$ and

The first requirement for these solutions to be real is that

which gives the transition curve between spherically stable and conditionally stable regions. This threshold is represented by thick dashed curves in figures 5–7.

Figure 8.

Amplitude-independent quantity $g_{0}$ as a function of the frequency $\unicode[STIX]{x1D714}_{0,n}$ , here for $n=3$ . Note that $g_{0}=0$ for $\unicode[STIX]{x1D714}_{0,n}=\unicode[STIX]{x1D6FA}^{\ast }\simeq 1$ .

The second requirement to get real steady-state solutions is that

which depends both on the sign of $g_{0}$ and on the sign of the two branches. $g_{0}$ is a frequency-dependent quantity which cancels for $\unicode[STIX]{x1D714}_{0,n}=\unicode[STIX]{x1D6FA}^{\ast }\simeq 1$ , and would equal $1$ when neglecting viscous effects. From a graphical observation of $g_{0}$ , as shown in figure 8, it directly comes that $g_{0}>0$ when $\unicode[STIX]{x1D714}_{0,n}>\unicode[STIX]{x1D6FA}^{\ast }$ while $g_{0}<0$ when $\unicode[STIX]{x1D714}_{0,n}<\unicode[STIX]{x1D6FA}^{\ast }$ . Also note that for $\unicode[STIX]{x1D714}_{0,n}=\unicode[STIX]{x1D6FA}^{\ast }$ , setting $g_{0}=0$ in (4.17) leads to $A_{th}^{\prime }=A_{th}$ . This corresponds to the meeting points of the two transition curves in figures 5–7. On the other hand, as soon as the amplitude $A$ fulfils the condition (4.17), the upper branch is always positive, while the lower one is conditionally positive, namely if

Because $h_{i}$ and $g_{0}^{2}$ are all positive quantities, equation (4.19) is only fulfilled if

which naturally corresponds to the absolute stability threshold defined by (4.13). The conditions for real steady-state solutions are summarised as follows:

1. (i) If the driving amplitude is below the conditional stability threshold defined by (4.17), thus both branches are complex and the only possible solution is the trivial one, i.e.  $S_{n}=0$ .

2. (ii) If the driving amplitude lies between the boundaries

   (4.21) $$\begin{eqnarray}\frac{1}{h_{3}}\sqrt{(h_{1}-g_{0}/h_{2})(h_{1}+g_{0}/h_{2})}<A<\frac{h_{1}}{h_{3}},\end{eqnarray}$$
   then both solutions are real only if $g_{0}>0$ ( $\unicode[STIX]{x1D714}_{0,n}>\unicode[STIX]{x1D6FA}^{\ast }$ ). In that case, the upper branch is stable while the lower one is unstable.

3. (iii) If the driving amplitude exceeds the absolute stability threshold defined by (4.13), then the upper branch is the only real solution when $g_{0}>0$ , while the lower one is the only real solution when $g_{0}<0$ .

Figure 9.

Amplitude of the parametrically forced shape mode (a) and the spherical mode (b) in steady state as a function of the driving amplitude, here for $n=3$ . The driving amplitude is divided by the amplitude threshold. Straight (respectively dashed) lines indicate stable (respectively unstable) solutions. The different curves correspond to different frequency values: $\unicode[STIX]{x1D714}_{0,n}=0.9$ (♢, $g_{0}<0$ ),  $1$ (○, $g_{0}<0$ ),  $1.05$ (▿, $g_{0}>0$ ) and $1.1$  (▫, $g_{0}>0$ ). In (b) the amplitude of the spherical mode is divided by the driving amplitude. The inset graph in (a) shows the curves in the close vicinity of the parametric threshold and a comparison with experimental values (●) that can be found in the work of Guédra et al. (Reference Guédra, Cleve, Mauger, Blanc-Benon and Inserra2017) with $R_{0}\simeq 73~\unicode[STIX]{x03BC}\text{m}$ and $\unicode[STIX]{x1D714}_{0,3}\simeq 0.87$ .

Figure 10.

Oscillations of the spherical mode and the seven first shape modes for $n=3$ , $\unicode[STIX]{x1D714}_{0,n}=1.05$ and $A/A_{th}=0.95$ , corresponding to the green ▿ curve shown in figure 9. (a,b) Present numerical simulations obtained from two different initial conditions on the parametrically forced shape mode, $s_{3}(0)=10^{-2}$ (a) and $s_{3}(0)=0.2$ (b), illustrating the hysteretic behaviour in the region bounded by the conditional threshold and the absolute threshold. (c) Focuses on few acoustic periods during the steady state obtained for the second initial condition, providing a comparison between numerical simulations (blue curves) and first-order approximate solutions (orange curves).

An illustration of these stability conditions is provided in figure 9(a), which gives the steady-state solutions (4.15) as a function of the driving amplitude. The solutions are given here for the $n=3$ shape mode. Among the four curves, two of them correspond to $\unicode[STIX]{x1D714}_{0,n}=0.9$ and $\unicode[STIX]{x1D714}_{0,n}=1$ , so that $\unicode[STIX]{x1D714}_{0,n}<\unicode[STIX]{x1D6FA}^{\ast }$ thus $g_{0}<0$ . In this case, the only non-zero stable solution is obtained when $A>A_{th}$ , and the transition from spherically stable oscillations to non-spherically stable oscillations is continuous. Indeed, it is easy to see that the real lower branch defined by (4.15) leads to $S_{n}=0$ when $A=A_{th}$ . The two other curves correspond to $\unicode[STIX]{x1D714}_{0,n}=1.05$ and $\unicode[STIX]{x1D714}_{0,n}=1.1$ , so that $\unicode[STIX]{x1D714}_{0,n}>\unicode[STIX]{x1D6FA}^{\ast }$ thus $g_{0}>0$ . In this case, a region lies between the boundaries $A_{th}^{\prime }$ and $A_{th}$ where the spherical and non-spherical oscillations are conditionally stable. A sharp transition exists between the two oscillating regimes, and the real upper branch defined by (4.15) leads to $S_{n}=\sqrt{2g_{0}}/h_{2}$ for $A=A_{th}^{\prime }$ and $S_{n}=2\sqrt{g_{0}}/h_{2}$ for $A=A_{th}$ . The stability of both the trivial solution and the upper branch solution in this region leads to an hysteretic behaviour which indicates that stable non-spherical oscillations can be reached even for driving amplitudes below the absolute parametric threshold. To illustrate this, we report in figure 10 two direct numerical simulations of (2.2) and (2.3) corresponding to the green $\triangledown$ curve conditions shown in figure 9 and for a driving amplitude $A/A_{th}=0.95$ . The two numerical solutions are obtained for two different initial conditions on the $n=3$ shape mode. It demonstrates that two different steady states can be reached depending on the initial conditions, namely stable spherical oscillations (figure 10 a), corresponding to the damping of all shape modes (trivial solution $S_{n}=S_{m}=0$ ), and sustained shape oscillations (figure 10 b) corresponding to the upper branch shown in figure 9. In figure 10(c) is shown a comparison between the numerical solution and the first-order approximate solution in the steady-state regime, proving the consistency of the perturbation method to reproduce the complex structure of the nonlinear coupled oscillations.

Figure 11.

Amplitude of the parametrically forced shape mode $S_{n}$ in steady state as a function of the frequency $\unicode[STIX]{x1D714}_{0,n}$ , for $n=3$ . The thick black curve represents the amplitude at the absolute threshold, namely $S_{n}=0$ for $g_{0}<0$ and $S_{n}=2\sqrt{g_{0}}/h_{2}$ for $g_{0}>0$ . The other curves correspond to different driving amplitude ratios above threshold: $A/A_{th}=1.1$ (▫), $1.5$ (▿) and $2$ (○).

Oscillations of the shape mode $n$ with the complex amplitude $\bar{S}_{n}$ also imply a change in the amplitude of the spherical mode which is mathematically made explicit by (4.5). Figure 9(b) shows the corresponding amplitude of the spherical oscillations in steady state as a function of the driving amplitude, for the same conditions as in the left part. The ratio $X/A$ is plotted, so that for driving amplitudes below the absolute threshold (i.e. no shape oscillations, $S_{n}=0$ ) the linear relationship between the spherical response of the bubble and the external pressure naturally leads to the constant value $1/|\bar{\unicode[STIX]{x1D6E5}}_{0}|$ . However sustained shape oscillations of finite amplitude lead to a global decrease of the amplitude of the spherical mode. This is clearly observed in figure 9 and can be also noticed on the temporal solutions shown in figures 1–3. For the same reasons as mentioned previously, an hysteretic loop is also obtained when $\unicode[STIX]{x1D714}_{0,n}>\unicode[STIX]{x1D6FA}^{\ast }$ for driving amplitudes bounded by the conditional stability and absolute stability thresholds.

Figure 12.

Amplitude of the secondary shape modes $S_{m}$ in steady state as a function of the frequency $\unicode[STIX]{x1D714}_{0,n}$ , for $n=3$ . See the caption of figure 11 for the legend.

As the final point of the discussion, we recall that in the framework of the present analysis, the saturation of the parametrically forced shape mode $n$ can only come from nonlinear interactions involving secondary shape modes $m\neq n$ . As soon as the amplitude of the shape mode $n$ in steady state is determined, the amplitudes of all the secondary modes are directly deduced from (4.6). The figures 11 and 12 show the variations of the shape mode amplitudes in the vicinity of the first parametric resonance of the $n=3$ mode, for different values of the driving amplitude at and above the absolute threshold. The threshold curves ( $A=A_{th}$ , thick black curves) remind us that $S_{n}=0$ is the only stable solution for $\unicode[STIX]{x1D714}_{0,n}<\unicode[STIX]{x1D6FA}^{\ast }$ , while $S_{n}\neq 0$ at threshold for $\unicode[STIX]{x1D714}_{0,n}>\unicode[STIX]{x1D6FA}^{\ast }$ . The figure 12 also reveals that secondary odd shape modes have vanishing amplitudes at the lowest order ( $S_{5}=S_{7}=0$ ), while these modes are after all excited through nonlinear coupling as illustrated in figures 2 or 10. First, this means that only even modes participate to the nonlinear saturation of the parametrically excited shape mode ( $n=3$ ) at the lowest order. Second, it implies that oscillations of additional odd modes thus only appear at first order, following solution (3.10), which explains why these modes mainly oscillate at the frequencies $\unicode[STIX]{x1D70F}$ or $3\unicode[STIX]{x1D70F}$ rather than at the driving frequency $2\unicode[STIX]{x1D70F}$ . This is clearly observed on the temporal solutions for $s_{5}$ and $s_{7}$ given in figures 2 and 10. It is interesting to note that the perturbation analysis allows us to capture this specific behaviour although it goes beyond the validity framework of the asymptotic approach expressed by (3.6).