3.1. Dimensionless numbers

Both dimensional and dimensionless parameters provide useful information in this problem. This paper focuses on the surfactant effect at relatively low concentrations, for which the liquid surface tension is similar to that of water. Therefore, $R$ , $\rho$ , ${\unicode{x03BC}}$ and $\gamma _c$ can be regarded as characteristic quantities in the presence of the surfactant. As mentioned in § 2, the gas dynamics inside the bubble have negligible effects, allowing us to eliminate the gas density and viscosity from the analysis.

The following dimensionless numbers can be defined based on the characteristic quantities: the Bond number $B=\rho g R^2/\gamma _c$ and the Galilei number $\textit{Ga}=\rho g^{1/2} R^{3/2}/\mu$ . As explained later, it is convenient to replace $B$ with $\mathcal{C}\equiv B^3\textit{Ga}^{-4}$ , which takes a constant value for a solvent–surfactant system. The effect of the monolayer is quantified through the dimensionless concentration $\hat {c}=c_{\infty }/c_{\textit{cmc}}$ ( $c_{\textit{cmc}}$ is the critical micelle concentration) and the set of dimensionless parameters $\{\mathcal{P}_i\}$ characterising the solvent–surfactant system.

The following dimensionless parameters must be considered to characterise the physical properties of the fast surfactant monolayer: the bulk and surface Péclet numbers, $\textit{Pe}=\ell _c v_c/\mathcal{D}$ and $\textit{Pe}_s=\ell _c v_c/\mathcal{D}_{S}$ , the dimensionless depletion length $\varLambda _d=L_d/\ell _c$ , and the Marangoni (elasticity) number Ma = $\varGamma _{\infty } R_g T/\gamma _c$ . In these expressions, $\ell _c=\mu ^2/(\rho \gamma _c$ ) and $v_c=\gamma _c/\mu$ are the (intrinsic) viscous-capillary length and velocity, respectively. As explained in § 3.2, sorption kinetics is characterised by the single parameter $\varLambda _d$ in the fast-surfactant case (see also § 2). Table 1 shows the values of the intrinsic dimensionless numbers characterising the surfactant.

Table 1. Values of the dimensionless numbers involving the physical properties of the fluids and surfactant monolayer.

We consider the following dimensionless dependent parameters: the Reynolds number ${Re}=\rho v_t R/\mu$ , where $v_t$ is the terminal velocity (the vertical velocity in the base flow); the aspect ratio $\chi =a/b$ , where $a$ and $b$ are the half-length and half-breadth of the cross-sectional shape, respectively; and $\mathcal{S}=S/R^2$ , where $S$ is the cross-sectional area. For a given solvent–surfactant system, $\mathcal{C}$ is constant and, therefore,

(3.1) \begin{equation} {Re}={Re}(\textit{Ga},\hat {c};\{\mathcal{P}_i\}), \end{equation}

(3.2) \begin{equation} \chi =\chi ({Re},\hat {c};\{\mathcal{P}_i\}), \quad \mathcal{S}=\mathcal{S}({Re},\hat {c};\{\mathcal{P}_i\}), \end{equation}

where we have considered that $\textit{Ga}=\textit{Ga}({Re},\hat {c};\{\mathcal{P}_i\})$ in (3.2).

Consider the drag force $F_D$ experienced by the bubble in a steady vertical motion. The dimensionless force $\mathcal{F}_D=F_D/(\rho g R^3)$ is $\mathcal{F}_D=4\pi /3$ . We define the drag coefficient $C_D=C_D({Re},\hat {c};\{\mathcal{P}_i\})$ as

(3.3) \begin{equation} C_D=\frac {F_D}{1/2\, \rho v_t^2 S}=\frac {2\mathcal{F}_D}{\textit{Ga}^{-2}\, {Re}^2 \mathcal{S}}. \end{equation}

3.2. Dimensionless governing equations

It is instructive to rewrite the governing equations in their dimensionless form. The equations can be made dimensionless by measuring the lengths, velocities, time, stresses, surfactant surface and volumetric concentrations, and surfactant flux in terms of the intrinsic quantities $\ell _c$ , $v_c$ , $\ell _c/v_c$ , $\rho v_c^2$ , $\varGamma _{\infty }$ , $\varGamma _{\infty }/\ell _c$ , $\varGamma _{\infty }/t_c$ , respectively. In this section, all the symbols represent the dimensionless counterpart of the previously defined dimensional quantity. The non-dimensional continuity and momentum equations are

(3.4) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}=0, \end{equation}

(3.5) \begin{equation} \frac {\mathrm{D}\boldsymbol{v}}{\mathrm{D}t}=-\!\left (\mathcal{C}+\frac {\mathrm{d}^2h}{{\rm d}t^2}\right )\, \boldsymbol{e_z}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }, \quad \boldsymbol{\sigma }=-p \boldsymbol {I}+\boldsymbol{\tau },\quad \boldsymbol{\tau }=\boldsymbol{\nabla }\boldsymbol {v}+\left (\boldsymbol{\nabla }\boldsymbol {v}\right )^T \!, \end{equation}

where $\mathcal{C}$ is a dimensionless parameter that takes a constant value for a fixed solvent–surfactant system, as mentioned in the previous section. The kinematic compatibility condition is the same as (2.5). The equilibrium of normal and tangential stresses on the interface leads to the equations

(3.6) \begin{equation} \boldsymbol{n}\boldsymbol{\cdot }{\boldsymbol \sigma }\boldsymbol{\cdot }\boldsymbol{n}=\widetilde {\gamma } \kappa ,\quad \boldsymbol{t_1}\boldsymbol{\cdot }{\boldsymbol \sigma }\boldsymbol{\cdot }\boldsymbol{n}=\tau _{\textit{Ma}}^{(1)}, \quad \boldsymbol{t_2}\boldsymbol{\cdot }{\boldsymbol \sigma }\boldsymbol{\cdot }\boldsymbol{n}=\tau _{\textit{Ma}}^{(2)}, \end{equation}

where $\widetilde {\gamma }=\gamma /\gamma _c$ , $\tau _{\textit{Ma}}^{(1)}=-{ \boldsymbol{\nabla} }_{\!S}\widetilde {\gamma }\boldsymbol{\cdot }\boldsymbol{t_1}$ and $\tau _{\textit{Ma}}^{(2)}=-{ \boldsymbol{\nabla} }_{\!S}\widetilde {\gamma }\boldsymbol{\cdot }\boldsymbol{t_2}$ . The monomer conservation equation in the bulk and the Langmuir adsorption isotherm are

(3.7) \begin{equation} \frac {\partial c}{\partial t}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }c=\textit{Pe}^{-1}\boldsymbol{\nabla }^2 c, \quad \frac {\varGamma }{1-\varGamma }-\varLambda _d c_s=0. \end{equation}

The surfactant surface advection–diffusion equation with $\mathcal{D}_S=\text{const.}$ reads

(3.8) \begin{equation} \frac {\partial \varGamma }{\partial t}+{ \boldsymbol{\nabla} }_{\!S}\boldsymbol{\cdot }(\varGamma \boldsymbol { v}_S)+\varGamma ({ \boldsymbol{\nabla} }_{\!S}\boldsymbol{\cdot }\boldsymbol{n})(\boldsymbol{ v}\boldsymbol{\cdot }\boldsymbol{n})=\textit{Pe}_s^{-1}{ {\nabla} }_S^2 \varGamma +\textit{Pe}^{-1}\left . { \boldsymbol{\nabla }} c\right |_n.\end{equation}

The Langmuir equation of state is

(3.9) \begin{equation} \widetilde {\gamma }=1+\textit{Ma}\,\ln \left (1-\varGamma \right )\!. \end{equation}

The dimensionless counterparts of the boundary conditions must also be considered, including the condition $c=c_{\infty }$ at the outer surface. Finally, the bubble volume equation becomes

(3.10) \begin{equation} \frac {4}{3}\mathcal{C}^{-1}\,\textit{Ga}^2=\int _0^{1} \frac {r_i^2}{[1+({\rm d}r_i/{\rm d}z)^2]^{1/2}}\, {\rm d}s. \end{equation}

The concentration $c_{\infty }$ measured in terms of $\varGamma _{\infty }/\ell _c$ is replaced by $c_{\infty }/c_{\textit{cmc}}$ , which is more convenient experimentally. Under this consideration, the above-mentioned equations involve the same dimensionless numbers $\{\mathcal{C}$ , $\varLambda _d$ , $\textit{Ma}$ , $\textit{Pe}$ , $\textit{Pe}_S$ ; $\textit{Ga}$ , $c/c_{\textit{cmc}}\}$ defined in the previous section. In our analysis, the values of $\{\mathcal{C}$ , $\varLambda _d$ , $\textit{Ma}$ , $\textit{Pe}$ , $\textit{Pe}_S\}$ are fixed, while $\{\textit{Ga}$ , $c/c_{\textit{cmc}}\}$ are varied.

The ‘fast kinetics’ limit considered in this work can be derived from the general dimensionless equations. In these equations, the adsorption and desorption fluxes are

(3.11) \begin{equation} \mathcal{J}_{a}=\textit{Bi}\, c_s(1-\varGamma ),\quad {\mathcal J}_{d}^{(\,j)}=\frac {\textit{Bi}\,\varGamma }{\varLambda _d}, \end{equation}

where $\textit{Bi}=k_a/v_c$ is the Biot number. The conservation of surfactant exchanged between the interface and the sublayer leads to

(3.12) \begin{equation} {\mathcal J}_{a}^{(\,j)}-{\mathcal J}_{d}^{(\,j)}=\frac {\textit{Bi}(1-\varGamma )}{\varLambda _d}\left (\varLambda _d c_s-\frac {\varGamma }{1-\varGamma }\right )=\textit{Pe}^{-1}\left . \boldsymbol{\nabla }c\right |_n. \end{equation}

This equation replaces (3.7)-right in the general case. The ‘fast kinetics’ limit corresponds to $\textit{Bi}\to \infty$ and $[\varLambda _d c_s-\varGamma /(1-\varGamma )]\to 0$ so that $\textit{Pe}^{-1} \boldsymbol{\nabla }c|_n$ remains finite. Therefore, ${\mathcal J}_{a}^{(\,j)},{\mathcal J}_{d}^{(\,j)}\gg \textit{Pe}^{-1} \boldsymbol{\nabla }c|_n$ although the term $\textit{Pe}^{-1} \boldsymbol{\nabla }c|_n$ is not negligible in (3.8).

5.1. Types of paths

This section compares our experimental results for Surfynol with those obtained by Rubio et al. (Reference Rubio, Vega, Cabezas, Montanero, LÓpez-Herrera and Herrada2024) for SDS. For illustration, figure 6 shows four examples of the paths followed by bubbles in the presence of Surfynol. The corresponding bubble vertical velocity $v_z$ and aspect ratio $\chi$ are plotted in figure 7 as a function of the vertical coordinate $z$ .

Figure 6. Bubble trajectory for (a) $\textit{Ga}=79$ and $c_{\infty }/c_{\textit{cmc}}=10^{-5}$ , (b) $\textit{Ga}=79$ and $c_{\infty }/c_{\textit{cmc}}=10^{-4}$ , (c) $\textit{Ga}=58$ and $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ , and (d) $\textit{Ga}=66$ and $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ . The graphs also show the projections of the trajectories onto the planes $(x,y)$ , $(x,z)$ and $(y,z)$ .

Figure 7. Bubble vertical velocity $v_z$ and aspect ratio $\chi$ as a function of the vertical coordinate $z$ for the cases in figure 6.

Figure 6(a) corresponds to a helical motion characterised by a double-threaded wake (two counter-rotating vorticities) (Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012). Panel (b) shows a zigzag motion obtained when the surfactant concentration increases for the same Galilei number (bubble radius). A vortex is shed periodically from the wake during the bubble rising (Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012). Oscillatory instabilities reduce the mean vertical velocity because energy transfers from the base flow to the unstable mode. This is observed in panels (a) and (b), in which the velocity decreases at $z\simeq 110$ mm due to the growth of the helical and zigzag motions (figure 7).

Figure 6(c) shows an oblique path with a tilt angle approximately equal to 2.5 $^{\circ }$ . This path results from the constant lift force generated by a pair of semi-infinite counter-rotating streamwise vortices (no vortex shedding occurs in this case) (Tagawa et al. Reference Tagawa, Takagi and Matsumoto2014). Oblique paths are a characteristic effect of the surfactant monolayer. They have been reported neither in experiments nor in simulations with clean water (Bonnefis et al. Reference Bonnefis, Sierra-Ausin, Fabre and Magnaudet2024). These paths typically occur when the surfactant monolayer immobilises the bubble surface (Tagawa et al. Reference Tagawa, Takagi and Matsumoto2014), which resembles the path instability of light spheres (Jenny et al. Reference Jenny, Bouchet and Dusek2003, Reference Jenny, Dusek and Bouchet2004). One expects this instability to set in when both a sufficient vertical velocity has been reached and a sufficiently long wake has developed (Jenny, Dusek & Bouchet Reference Jenny, Dusek and Bouchet2004). In the experiment of figure 6(c), this occurs at a distance from the ejector $z\gtrsim 100$ mm, much larger than that at which the terminal velocity is reached ( $z\simeq 25$ mm) (figure 7 b). A significant decrease in the aspect ratio accompanies the growth of the stationary instability. The tiny oscillation observed for $z\lesssim 100$ mm disappears as the stationary instability grows. The initial conditions select the symmetry plane.

Finally, figure 6(d) shows the transition from a zigzag to an oscillatory oblique path. The path tilt partially suppresses the bubble oscillations, slightly increasing the vertical velocity (figure 7 b). The discrete Fourier transform of $v_z(t)$ exhibits a peak at the frequency $\omega \simeq 15$ Hz. As shown in § 6.1, this frequency approximately corresponds to that of the unstable oscillatory mode in the secondary instability obtained in the global stability analysis. Oscillatory oblique paths are also observed when a light solid sphere rises in water (Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012).

The local velocity profiles exhibit the so-called ‘overshooting’ phenomenon: the velocity reaches a maximum and subsequently decreases to its terminal value (figure 7). The non-monotonous behaviour of the bubble aspect ratio accompanies this phenomenon. Axisymmetric transient simulations have allowed us to explain the overshooting with non-fast surfactants in terms of the extra drag resulting from the late growth of the dynamic surfactant layer (Fernández-Martínez et al. Reference Fernández-Martínez, Cabezas, López-Herrera, Herrada and Montanero2025). Interestingly, the same phenomenon occurs with Surfynol despite its much faster sorption kinetics. This suggests that surfactant diffusion prevents the full growth of the dynamic surfactant layer during the bubble acceleration.

Figure 8 shows experimental results for SDS and Triton X-100 obtained by Fernández-Martínez et al. (Reference Fernández-Martínez, Cabezas, López-Herrera, Herrada and Montanero2025) for the same values $\textit{Ga}=66$ and $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ as those in figure 7(b). The time interval (length) over which the overshooting occurs increases as the desorption rate decreases. This explains why the overshooting extends over a much longer time with Triton. The comparison between these results and those of Surfynol indicates that Surfynol desorption takes place much faster (overshooting occurs for $z\lesssim 20$ mm in this case). Bubble rising can be used as a testbed to evaluate qualitatively the desorption rate of a fast surfactant.

Figure 8. Bubble velocity $v_z$ and aspect ratio $\chi$ as a function of the vertical position $z$ of the centre of gravity for $\textit{Ga}=66$ and $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ (Fernández-Martínez et al. Reference Fernández-Martínez, Cabezas, López-Herrera, Herrada and Montanero2025). The solid and dashed lines correspond to the stable and oscillatory parts of the bubble trajectory, respectively.

The experiment with $\textit{Ga}=43$ and $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ exhibits a slightly oblique path (figure 9). The tilt angle is approximately 0.4 $^{\circ }$ , smaller than the threshold (0.5 $^{\circ }$ ) adopted to categorise the path as oblique. Unlike in figure 6(c), the path is inclined right after the bubble ejection and the aspect ratio does not significantly change during the bubble rising. As shown in § 6.1, the global stability analysis predicts a stable path for these experimental conditions. The bubble radius ( $R=0.57$ mm) is the smallest in our experiments. The bubble rising is the slowest and the surfactant monolayer is practically formed when the maximum velocity is reached, despite the relatively large surfactant concentration $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ . Consequently, the overshooting phenomenon is hardly observed (figure 9).

Figure 9. (a) Bubble velocity $v_z$ and aspect ratio $\chi$ as a function of the vertical position $z$ of the centre of gravity for $\textit{Ga}=43$ and $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ , and (b) its corresponding trajectory.

5.2. Transition from clean to immobilised interface

Consider the Reynolds number (the terminal velocity) measured for a fixed value of the Galilei number (bubble radius) in the presence of SDS (figure 10 a). The surfactant produces small effects for low surfactant concentrations. In this regime, the Reynolds number is slightly smaller than that corresponding to the surfactant-free case. As the surfactant concentration increases, the Reynolds number decreases to a value approximately equal to that corresponding to a solid sphere. In this regime, the surfactant distributes smoothly along the interface, which is immobilised by the Marangoni stress (Rubio et al. Reference Rubio, Vega, Cabezas, Montanero, LÓpez-Herrera and Herrada2024).

Figure 10. Reynolds number Re as a function of the surfactant concentration $c_{\infty }/c_{\textit{cmc}}$ for (a) SDS and (b) Surfynol. The solid and open symbols correspond to stable and unstable realisations, respectively. The labels S, O and S/O indicate whether the instability is stationary, oscillatory or a combination of both. The error bars in the SDS data correspond to the standard deviation and show the high reproducibility of our experimental method.

One can identify an intermediate regime corresponding to surfactant concentrations between those of nearly clean bubble and practically immobilised interface. In this intermediate regime, SDS is absorbed in the bubble front and transported by convection towards the rear, where it is desorbed. Convection and surface diffusion compete within a short interface portion, forming a surfactant surface boundary layer (Rubio et al. Reference Rubio, Vega, Cabezas, Montanero, LÓpez-Herrera and Herrada2024). Only the interface behind that layer is immobilised (the stagnant cap model). The size of the immobilised interface portion increases with the surfactant concentration until it occupies the entire bubble surface (the immobilised interface regime).

In the case of Surfynol, the Reynolds numbers corresponding to the immobilised interface regime are reached at much smaller values of $c_{\infty }/c_{\textit{cmc}}$ (figure 10 b). This can be explained in terms of the surfactant behaviour at equilibrium. The Surfynol depletion length is much larger than that of SDS, which implies that more Surfynol is expected to be adsorbed during the bubble rising (figure 4). In addition, Surfynol is stronger than SDS for small surface coverages (figure 5). Both effects increase the Marangoni stress and, consequently, the drag force.

The transition from the clean bubble behaviour to the immobilised interface regime occurs within a narrower interval of $c_{\infty }/c_{\textit{cmc}}$ in the Surfynol case. For instance, consider the experiments with $\textit{Ga}=66$ . The Reynolds number decreases approximately from the clean surface value ${Re}=258.7$ at $c_{\infty }/c_{\textit{cmc}}=10^{-5}$ to the immobilised interface value ${Re}=116.6$ at $c_{\infty }/c_{\textit{cmc}}=5\times 10^{-5}$ . Conversely, several experimental realisations with intermediate Reynolds numbers were found in the presence of SDS (figure 10). This difference can be explained as follows. The Surfynol surface concentration is practically at equilibrium with the volumetric one at any interface point. Therefore, a surfactant surface boundary layer entails large gradients of the volumetric concentration next to the interface in the streamwise direction. However, strong surfactant convection parallel to the bubble surface reduces these gradients. Therefore, the surfactant surface boundary layer and stagnant caps, characteristic of the above-mentioned intermediate regime, appear under more stringent parameter conditions with fast surfactants.

Figure 11. (a) Aspect ratio $\chi$ and (b) Reynolds number Re as a function of the Galilei number Ga for SDS (grey symbols) and Surfynol (blue symbols). The solid and open symbols correspond to stable and unstable realisations, respectively. The labels S, O and S/O indicate whether the instability is stationary, oscillatory or a combination of both. The error bars correspond to the standard deviation and show the high reproducibility of our experimental method. The upper and lower solid lines are the function Re(Ga) for a clean bubble (Herrada & Eggers Reference Herrada and Eggers2023) and a solid hollow sphere, respectively.

Figure 11 also shows the difference between SDS and Surfynol discussed previously. A sharp transition occurs at the Surfynol concentration $c_{\infty }/c_{\textit{cmc}}\simeq 5\times 10^{-5}$ . For $c_{\infty }/c_{\textit{cmc}}\leq 10^{-5}$ , the bubble shape significantly deviates from the sphere and the Reynolds number values are relatively close to those of the clean bubble. Both $\chi$ and Re exhibit a non-monotonous dependency on Ga. For $c_{\infty }/c_{\textit{cmc}}\geq 5\times 10^{-5}$ , the bubble adopts a quasi-spherical shape ( $\chi \lesssim 1.1$ ) and the Reynolds number sharply decreases. In this concentration regime, both $\chi$ and Re increase linearly with Ga.

The upper and lower solid lines in figure 11(b) are the function Re(Ga) for a clean bubble (Herrada & Eggers Reference Herrada and Eggers2023) and a solid hollow sphere (Clift, Grace & Weber Reference Clift, Grace and Weber1978), respectively. For Surfynol concentrations $c_{\infty }/c_{\textit{cmc}}\geq 5\times 10^{-5}$ , the Reynolds number falls below the solid line corresponding to a solid sphere. This can be explained by the drag force increase due to the small bubble deformation ( $\chi \simeq 1.1$ ) and the velocity decrease caused by the oscillatory instability.

Figure 12 shows the projection of the experimental data on the parameter plane (Re, $C_D^*$ ). Here, $C_D^*$ is the normalised drag coefficient defined as

(5.1) \begin{equation} C_D^*=\frac {C_D-C_{D,\textit{clean}}}{C_{D,\textit{rigid}}-C_{D,\textit{clean}}}, \end{equation}

where $C_{D,\textit{clean}}$ and $C_{D,\textit{rigid}}$ are the drag coefficients of a spherical bubble with the free-slip condition (clean interface) (Mei, Klausner & Lawrence Reference Mei, Klausner and Lawrence1994) and no-slip condition (rigid interface) (Clift et al. Reference Clift, Grace and Weber1978), respectively. It must be noted that the effect of the bubble shape on $C_D$ , $C_{D,\textit{clean}}$ and $C_{D,\textit{rigid}}$ is much smaller than the influence of the interface condition for the range of Reynolds numbers considered in our analysis (Tagawa et al. Reference Tagawa, Takagi and Matsumoto2014). Therefore, the coefficient $C_D^*$ essentially indicates the degree of slip on the bubble surface. The free-slip and no-slip conditions correspond to $C_D^*=0$ and 1, respectively.

Figure 12. Normalised drag coefficient $C_D^*$ as a function of the Reynolds number Re. The solid and open symbols correspond to stable and unstable realisations, respectively. The labels S, O and S/O indicate whether the stability is stationary, oscillatory or a combination of both. The arrows indicate the critical Reynolds numbers mentioned in the text.

Consider the experiments with Surfynol. For the higher surfactant concentrations, $C_D^*$ is significantly larger than unity. This effect is smaller without oscillations (solid symbols) and is explained by the less aerodynamic shape adopted by the bubble ( $\chi \gt 1$ ). The bubble oscillation reduces the terminal velocity, increasing the drag coefficient up to 1.5. The magnitude of this effect depends on the type and amplitude of the oscillation. The same phenomenon is observed with SDS.

5.3. Path instability transition

Now, we focus on path instability. In figures 10–12, we distinguish between the oscillatory (O) instability leading to either helical or zigzag motion and the stationary (S) instability leading to an oblique path with a tilt angle greater than 0.5 $^{\circ }$ . There are several experimental realisations where the oscillatory and stationary instabilities coexist. These realisations are denoted with the symbol O/S.

Significant differences exist between the path stability in the presence of SDS and Surfynol. In the case of SDS, oscillatory instability arises for Reynolds numbers larger than those in the immobilised interface regime. This does not occur in the Surfynol case, where unstable realisations can be found only at Reynolds numbers similar to that of the immobilised interface case (figure 11).

All the experimental realisations with relatively large concentrations ( $c_{\infty }/c_{\textit{cmc}}\geq 10^{-3}$ ) of SDS are unstable. Conversely, path instability arises at a critical Galilei number for the Surfynol concentration $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ (figure 10). As explained in § 6, this is a critical characteristic of Surfynol because it allows us to study the instability transition from the global stability analysis.

In the SDS case, the stability character always exhibits a monotonous dependence with respect to the surfactant concentration: for a fixed value of Ga (bubble radius), instability arises above a critical surfactant concentration. Interestingly, this does not occur with Surfynol for $\textit{Ga}=53$ (figure 10). In this case, the path becomes unstable at $c_{\infty }/c_{\textit{cmc}}=5\times 10^{-5}$ and stability is recovered for larger concentrations.

The critical Reynolds number for the path stability of a solid sphere is approximately 102.5 (Jenny et al. Reference Jenny, Dusek and Bouchet2004). The instability corresponds to a non-oscillating growth of small non-axisymmetric perturbations, which gives rise to a stationary oblique path (S in our notation). A secondary Hopf bifurcation occurs at ${Re}=112$ , leading to an oblique oscillatory path (O/S in our notation). Finally, a low-frequency zigzag (O in our notation) arises at ${Re}=122.5$ (Jenny et al. Reference Jenny, Dusek and Bouchet2004). The above-mentioned critical Reynolds numbers are indicated by arrows at the line $C_D^*=1$ in figure 12. The same sequence of instabilities is observed with Surfynol for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ and $\textit{Ga}\geq 53$ (see supplementary movie 1). The surfactant practically immobilises the interface at this concentration, which explains why the sequence of instabilities is the same as that of a solid sphere (Jenny et al. Reference Jenny, Dusek and Bouchet2004). We analyse this sequence of instabilities from global linear stability in § 6.

6.1. Validation of the numerical solution

Table 3 compares the terminal velocity, aspect ratio and the stability character obtained from the global linear stability with their counterparts in the experiment. There is a good agreement for the three Galilei numbers. The numerical terminal velocities are slightly larger than the experimental values. The bifurcation sequence predicted by the global stability analysis is found in the experiments. The path is stable for $\textit{Ga}=53$ . A stationary instability arises at $\textit{Ga}=58$ . A secondary oscillatory Hopf bifurcation occurs at Ga $\simeq 66$ . The stationary and oscillatory instabilities coexist at $\textit{Ga}=66$ . This corresponds to the paths observed in the experiments: a straight path for $\textit{Ga}=53$ , an oblique path for $\textit{Ga}=58$ and oscillations around an oblique path for $\textit{Ga}=66$ . The oscillation frequency predicted by the linear stability analysis for $\textit{Ga}=66$ is 22 Hz, whereas the corresponding experimental value is approximately 15 Hz.

Table 3. Experimental and numerical results for different Galilei numbers around its critical value for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ . The table shows the terminal velocity $v_t$ , aspect ratio $\chi$ and stability character, both from the experiment (Exp.) and from the global stability analysis (GSA).

6.2. Base flow

We will return to the global stability analysis in § 6.3. Now, we examine the base flow for $\textit{Ga}=53$ , 58 and 66. Consider the marginally stable flow obtained for the critical Galilei number $\textit{Ga}=58$ . Figure 13 shows the flow streamlines and the volumetric surfactant concentration. As shown later, the surfactant monolayer practically immobilises the interface. The viscous boundary layer remains attached to the bubble front and separates from the bubble surface behind the bubble equator. Two recirculation cells arise in the wake. The surfactant concentration in the bubble front is lower than the upstream value  $c_{\infty }$ . The concentration increases and reaches its maximum value approximately at the boundary layer separation point (indicated by a circle in figure 13). As shown later, the surfactant is adsorbed at the bubble front, while desorption occurs beyond a certain point in front of the bubble’s equator. The surfactant excess is transferred back to the liquid by diffusion.

Figure 13. (a) Streamlines and (b) surfactant volumetric concentration $c/c_{\infty }$ for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ and $\textit{Ga}=58$ . The red circle indicates the viscous boundary layer separation point.

Figure 14 shows the distribution of the relevant quantities over the bubble surface. As the sublayer and the interface are at equilibrium, both the bulk (figure 14 a) and the surface concentration (figure 14 b) show the same behaviour. Both concentrations show a minimum at the front stagnation point and increase smoothly downstream. Their maximum values are reached behind the bubble equator ( $\alpha /\pi =0.63$ ), not at the rear of the bubble. This occurs due to the reverse flow in the wake, which convects the surfactant from the south pole to the equator.

Figure 14. (a) Volumetric surfactant concentration $c_s$ evaluated at the free surface, (b) surfactant surface concentration $\varGamma$ , (c) surfactant fluxes, (d) surface velocity $v_s$ , (e) Marangoni stress $\tau _{\textit{Ma}}^{(1)}$ , tangential viscous stress $\tau _{nt1}$ , and normal viscous stress $\tau _{nn}$ , and ( f) pressure $p$ as a function of the polar angle $\alpha$ for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ .

Figure 14(c) shows the values of the fluxes in the surfactant conservation equation (2.10); i.e. the flux of surfactant exchanged with the sublayer, $J=\mathcal{D} \boldsymbol{\nabla }c|_n$ , the flux convected over the interface, $J_c={\boldsymbol {\nabla} }_{\!S}\boldsymbol{\cdot }(\varGamma \boldsymbol {v}_s)$ , and the flux difussed across the interface, $J_{Ds}={\boldsymbol {\nabla} }_{\!S}\boldsymbol{\cdot }(\mathcal{D}_S{\boldsymbol{\nabla} }_S\varGamma )$ . In all the cases, positive (negative) values imply that the surfactant is gained (lost) in the interface element. Due to uneven surfactant distribution, surface tension gradients cause Marangoni stresses that almost immobilise the interface. However, surfactant convection plays a crucial role. Convection pushes the surfactant away from the bubble front. Surface diffusion does not compensate for convection. Therefore, the surfactant is transferred from the sublayer in this part of the bubble.

The velocity $v_s$ (figure 14 d) and the surfactant gradient ${\rm d}\varGamma /{\rm d}s$ decrease in the bubble south hemisphere, resulting in a reduction of the convection flux magnitude in that region (figure 14 c). This flux becomes positive, implying that convection accumulates surfactant in this part of the bubble. Most of the surfactant gained by convection is transferred back to the liquid. For this reason, the sublayer concentration grows downstream along the interface (figure 13 b), and the interface and sublayer remain at equilibrium.

Convection plays a key role in maintaining the variation of surfactant concentration over the interface, even though the surface velocity is three orders of magnitude smaller than the terminal velocity. The bulk and surface diffusion collaborate to reduce the surfactant gradient. Since the equilibrium equation (2.9) links $c_s$ and $\varGamma$ , both mechanisms transfer surfactant in the same direction over most of the interface. This explains why $J$ and $J_{Ds}$ have the same sign for most values of $\alpha$ (figure 14 c). This also explains why $J$ vanishes near the inflection point ${\nabla} _s^2 \varGamma =0$ ( $J_{Ds}=0$ ), implying that $J_c$ vanishes near that point too.

The Marangoni stress almost immobilises the interface. The tiny values taken by the surface velocity (figure 14 d) are crucial to understanding the surfactant behaviour. In the absence of surfactant convection (a fully immobilised interface), diffusion would render surfactant concentration uniform, eliminating the Marangoni stress responsible for interface immobilisation.

The Marangoni stress $\tau _{\textit{Ma}}^{(1)}$ is balanced by the viscous tangential stress $\tau _{nt1}$ (figure 14 e). The stress $\tau _{\textit{Ma}}^{(1)}$ and, therefore, $\tau _{nt1}$ vanish where the surfactant surface concentration is maximum ( ${\rm d}\varGamma /{\rm d}s={\rm d}\gamma /{\rm d}s=0$ ). We conclude that the boundary layer separates from the bubble surface where the concentration is maximum. The sign of $\tau _{\textit{Ma}}^{(1)}$ changes in the wake, opposing the reverse flow in that region. This prevents the interface from moving towards the bubble front except within a region very close to the bubble’s south pole (see the inset in figure 14 d). Both the low pressure in the wake (figure 14 f) and the viscous stress contribute to the drag.

Figure 14 also shows the results for a subcritical ( $\textit{Ga}=53$ ) and a supercritical ( $\textit{Ga}=66$ ) case. The bubble terminal velocity increases with Ga (the bubble radius) (figure 14 c), which enhances surfactant convection (figure 14 d). As mentioned previously, surface and bulk diffusion compensate for the surfactant depletion caused by convection in the bubble front. Therefore, the magnitude of these two mechanisms also increases with Ga (figure 14 d). This also explains why $c_s$ decreases at the front stagnant point as Ga increases (figure 14 a).

The Marangoni stress in terms of the dynamical pressure $\rho v_t^2$ decreases as Ga increases (figure 14 e), resulting in less interface immobilisation (figure 14 d). The variation of the Galilei number significantly affects neither the surfactant concentration inflection point ( $J_{Ds}=0$ and $J\simeq J_c\simeq 0$ ) (figure 14 c) nor the boundary layer separation point ( $\tau _{\textit{Ma}}^{(1)}=\tau _{nt1}=0$ ) (figure 14 e). In terms of the dynamical pressure force $\rho v_t^2\, S$ , drag decreases as Ga increases due to a decrease in viscous stress and a tiny increase in the pressure at the wake. The reverse flow near the bubble’s south pole disappears for the subcritical case $\textit{Ga}=53$ (see the inset in figure 14 d).

Overall, the base flow for the subcritical, critical and supercritical Galieli numbers hardly differ. The mechanisms responsible for the instability must be described by analysing the effect of the perturbation associated with the critical eigenmode.

6.3. Analysis of the instability

This section analyses the eigenmodes for both subcritical and supercritical Galieli numbers. For this purpose, the eigenfunctions are normalised by setting $\text{max}(\delta r_i)=1$ . As mentioned previously, the lengths, velocities and stresses are divided by $R$ , $v_t$ and $\rho v_t^2$ , respectively. Consequently, the forces and torques are divided by $\rho v_t^2 R^2$ and $\rho v_t^2 R^3$ , respectively.

Figure 15 shows the eigenvalues for the three cases considered in § 6.1. The arrows indicate the unstable eigenmodes arising as the Galilei number increases. As mentioned previously, the bubble follows a straight (stable) path for $\textit{Ga}=53$ ( $\omega _i\lt 0$ ). A stationary instability ( $\omega _r=0$ , $\omega _i\gt 0$ ) arises at $\textit{Ga}=58$ . A secondary oscillatory Hopf bifurcation ( $\omega _r\gt 0$ , $\omega _i\gt 0$ ) occurs at $\textit{Ga}=66$ .

Figure 15. Eigenvalues for $0\leq \omega _r\leq 0.5$ and $\omega _i\geq -0.1$ . The results were obtained for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ . The arrows indicate the unstable eigenmodes. The eigenvalues are made dimensioneless with the inertio-capillary time $t_{ic}=(\rho R^3/\gamma _c)^{1/2}$ .

It is interesting to compare the instability of a surfactant-laden bubble and a solid sphere. The surfactant practically immobilises the bubble surface. This explains why the bubble path suffers from the same sequence of instabilities as that of a solid sphere: a stationary instability followed by the coexistence of stationary and oscillatory instabilities. This sequence differs from that observed in clean water: an oscillatory instability followed by the coexistence of oscillatory and stationary instabilities (Bonnefis et al. Reference Bonnefis, Fabre and Magnaudet2023). Interestingly, stationary modes (oblique paths) have been reported neither in experiments nor in simulations with clean water. It has been speculated that the changes introduced in the base flow by the primary oscillatory mode prevent the stationary perturbation from growing (Bonnefis et al. Reference Bonnefis, Sierra-Ausin, Fabre and Magnaudet2024).

A fundamental difference exists between the surfactant-covered bubble and the solid sphere. The Marangoni stress responsible for the surface immobilisation is perturbed due to the surfactant concentration perturbation. Conversely, the no-slip boundary condition cannot be perturbed. Figure 16 shows the effect of fixing the surfactant concentration (and, therefore, the Marangoni stress) on the eigenvalues. In this case, all the variables are perturbed except for the surfactant concentration, implying that the Marangoni stress is that of the base flow at any instant. The growth rate of the critical eigenmode increases from 0.00259 to 0.00723, implying that the flow becomes more unstable. In other words, the perturbation of the surfactant concentration partially suppresses the critical mode growth. The perturbation of the surfactant concentration causes the perturbation of the surface tension and, therefore, of the Marangoni stress over the bubble surface. The outer viscous stress balances the Marangoni stress. This means that the surfactant-induced perturbation of the viscous stress at the bubble surface plays a stabilising role.

Figure 16. Eigenvalues for $0\leq \omega _r\leq 0.5$ and $\omega _i\geq -0.175$ . The results were obtained for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ and $\textit{Ga}=58$ . The triangles correspond to the results obtained by fixing the surfactant concentration (all the variables are perturbed except for the surfactant concentration). The eigenvalues are made dimensioneless with the inertio-capillary time $t_{ic}=(\rho R^3/\gamma _c)^{1/2}$ .

The perturbation of the bubble shape also plays an important role in the instability. Figure 17 shows how the bubble shape is perturbed by the critical mode of the unstable case $\textit{Ga}=58$ and $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ . Several parallels have been plotted to make the effect more visible. The deformation has been adapted for visualisation. The axes have been oriented so that the shape is not perturbed in the $xz$ plane and the maximum perturbation occurs in the $yz$ plane. For all quantities analysed, the perturbation is symmetric with respect to the $yz$ plane. Since the bubble is axisymmetric in the base flow, the projections of the parallels onto the $xz$ and $yz$ planes are straight lines. The perturbation produces the bubble displacement to the right (the positive $y$ axis direction) (figure 17 d) and the bubble deformation due to the rotation of the parallels around the $x$ axis (figure 17 c–d). This rotation is anticlockwise for the equator and the parallels above it. The rotation changes to clockwise slightly below the equator. As a result, the bubble surface expands in the displacement direction and shrinks in the opposite direction.

Figure 17. (a) Bubble shape in the base flow indicating the analysed parallels and meridians, (b) perturbed bubble shape, and (c, d) projection of the perturbed bubble shape onto the (c) $xz$ and (d) $yz$ planes. The results were obtained for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ and $\textit{Ga}=58$ .

Figure 18. Perturbation of the bubble shape, curvature $\delta \kappa$ , surfactant surface concentration $\delta \varGamma$ ( $-\delta \gamma$ ) and hydrostatic pressure $\delta p$ in different sections of the bubble surface. The magnitudes have been adapted for visualisation. The sphere in the upper-left corner indicates the meridian and the parallels considered in the figure. The results were obtained for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ and $\textit{Ga}=58$ .

The bubble shape perturbation alters the hydrostatic pressure distribution over the bubble shape. As shown later, the normal viscous stress associated with the perturbation practically vanishes for $\textit{Ga}=58$ . This means that the hydrostatic pressure perturbation $\delta p$ practically equals $-\delta p_c$ , where $p_c=\gamma \, \kappa$ is the capillary pressure. Therefore, the hydrostatic pressure is perturbed due to changes in the surface tension $\gamma$ and curvature  $\kappa$ , i.e. $\delta p\simeq -(\delta \gamma \kappa _0+\gamma _0\delta \kappa )$ . These two contributions are shown in figure 18. The comparison between $\delta p$ and $\delta \kappa$ suggests that $\gamma _0\, \delta \kappa$ constitutes the main contribution to $\delta p$ . In other words, the changes in the surface tension contribute to the hydrostatic pressure perturbation on the bubble surface to a lesser extent. The expansion of the bubble’s right side (figure 17) produces a reduction in the curvature, which contributes to decreasing the capillary pressure (increasing the hydrostatic pressure) (figure 18). The opposite occurs on the bubble’s left side. This effect is more significant at the bubble front where the base flow surfactant concentration is lower ( $\gamma _0$ is higher).

An important aspect of the problem is how the velocity perturbation alters the balance between the mechanisms of surfactant transport over the interface. Figure 19(a) shows the projection onto the $yz$ plane of the velocity field perturbation evaluated in that plane. The horizontal pressure gradient (figure 19 b) pushes the liquid surrounding the bubble to the left at both the front and rear of the bubble. The magnitude of the velocity perturbation (figure 19 c) increases in the front and rear of the bubble and reaches its maximum value in the wake. In an experiment, the velocity perturbation remains small in terms of the terminal velocity. However, it can be significant compared with the surface velocity because the interface is almost immobilised. Therefore, the flow perturbation can alter the delicate balance between the surfactant transport mechanisms described in § 6.2.

Figure 19. Perturbation in the $yz$ plane of (a) the velocity projection onto the $yz$ plane, (b) the hydrostatic pressure, (c) the magnitude of the velocity perturbation and (d) the surfactant volumetric concentration. Yellow (blue) corresponds to higher (lower) values of the corresponding quantity. The results were obtained for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ and $\textit{Ga}=58$ .

It is also instructive to describe the link between the perturbation of interfacial and volumetric variables. The pressure field perturbation (figure 19 b) reflects the hydrostatic pressure perturbation $\delta p$ discussed previously (figure 18). Figure 19(d) shows that the perturbation of the surfactant volumetric concentration is antisymmetric with respect to the $yz$ plane, as can be anticipated from $\delta \varGamma$ (figure 18). The perturbation of the surfactant volumetric concentration is confined within the wake.

We concluded previously that the interface curvature perturbation constitutes the main contribution to $\delta p$ at the bubble surface. However, the capillary pressure is also affected by the surface tension perturbation $\delta \gamma$ due to changes in the surfactant surface concentration. This concentration is altered by the surface expansion/contraction and the surfactant flux produced by the perturbed flow. The general view of the velocity field perturbation shows an anticlockwise rotation of the liquid around the bubble and the formation of a recirculation cell rotating in the opposite direction in the bubble rear (figure 19 a).

The perturbation of the velocity components at the bubble surface changes the surfactant surface concentration with respect to that in the base flow. Figure 20 shows the three components of the velocity field perturbation on the bubble surface. Here, $\delta v_n$ , $\delta v_{t1}$ and $\delta v_{t2}$ are the velocity components along the normal and two tangential directions to the bubble surface, respectively. We plot the results evaluated at the meridians $\theta =0$ and $\pi /2$ , as indicated in the figure. The normal component $\delta v_n$ is insignificant compared to the tangential component. The surfactant is dragged towards the left side at the bubble front ( $\alpha /\pi \lesssim 0.1$ ) and rear ( $\alpha /\pi \gtrsim 0.95$ ), as shown by the values of $\delta v_{t1}$ and $\delta v_{t2}$ . However, the surface contraction and the incoming surfactant flux increase the concentration on the left side (figure 18), reducing the surface tension and decreasing the capillary pressure against the effect of the curvature increase. Furthermore, the horizontal flow (the tangential component $\delta v_{t2}$ ) switches direction to the right due to the significant horizontal gradient in the concentration of the surfactant. At the equator, there is almost no pressure change as the effect of the curvature and surface concentration perturbation practically compensate for each other (figure 18). Interestingly, the angular direction of the horizontal flow (the sign of $\delta v_{t2}$ ) changes several times along the bubble height (the value of $\alpha$ ) for $\textit{Ga}=58$ .

Figure 20. Components $\delta v_n$ , $\delta v_{t1}$ and $\delta v_{t2}$ of the velocity perturbation at the bubble surface. The results were obtained for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ , and $\textit{Ga}=53$ and 58.

The description of the instability mechanism must ultimately rely on the forces and torques arising in the perturbed solution. The surface forces exerted by the liquid on the bubble are altered by the bubble deformation and by the perturbation of the hydrostatic pressure and viscous stresses on the bubble surface. The equivalent system of perturbed forces is a lateral force acting in the $-y$ direction and a clockwise torque ( $-x$ ). As shown later, this lateral force and torque have the same directions in the stable ( $\textit{Ga}=53$ ) and unstable ( $\textit{Ga}=58$ ) cases. Specifically, the lateral force pushes the bubble back towards the original position, while the torque reaction exerted on the liquid contributes to the rotating perturbed flow shown in figures 19 and 21 (note that the torque exerted on the liquid is anticlockwise). Therefore, the instability must be explained by a local effect.

Figure 21. (a) Perturbation in the $yz$ plane of the velocity projection on the $yz$ plane. (b) Magnitude of the velocity field perturbation. Yellow (blue) corresponds to higher (lower) values. The results were obtained for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ and $\textit{Ga}=53$ .

Figure 22. (a) Tangential $\delta \tau _{nt1}$ and normal $\delta \tau _{nn}$ elements of the perturbed viscous stress tensor evaluated at the meridian $\theta =\pi /2$ , and tangential element $\delta \tau _{nt2}$ evaluated at $\theta =0$ . Contributions to (b) the viscous lateral force $\delta F_{v,\delta \tau }$ and (c) torque $\delta M_{v,\delta \tau }$ (c) resulting from the perturbed viscous stress tensor. The results were obtained for $c_{\infty }/c_{\textit{cmc}}=10^{-3}$ and $\textit{Ga}=53$ and 58.

To analyse the source of instability, we calculated the contributions of the pressure ( $\delta F_p(\alpha )$ ) and viscous ( $\delta F_v(\alpha )$ ) force perturbations to the lateral force ( $\delta F(\alpha )$ ) exerted on the annular surface element of area ${\rm d}S$ (figure 22) as

(6.1) \begin{equation} \delta F=\delta F_p+\delta F_v, \quad \delta F_p=\delta F_{p,\delta p}+\delta F_{p,sh}, \quad \delta F_v=\delta F_{v,\delta \tau }+\delta F_{v,sh}, \end{equation}

where $\delta F_{p,\delta p}$ and $\delta F_{p,sh}$ stand for the contributions of the pressure and bubble shape perturbations to $\delta F_p$ . Analogously, $\delta F_{v,\delta \tau }$ and $\delta F_{v,sh}$ are the contributions of the viscous stress and bubble shape perturbations to $\delta F_v$ . We proceeded in a similar way to determine the $x$ component of the torque ( $\delta M(\alpha )$ ):

(6.2) \begin{equation} \delta M=\delta M_p+\delta M_v, \quad \delta M_p=\delta M_{p,\delta p}+\delta M_{p,sh}, \quad \delta M_v=\delta M_{v,\delta \tau }+\delta M_{v,sh}, \end{equation}

where the meaning of the subindexes is the same as that in (6.1). Figure S4 in the supplementary material shows all the results. Here, we summarise them. For the stable and unstable cases, $\delta F_{p,\delta p}$ is the main contribution (approximately 80 %) to the lateral force $\delta F$ . In addition, $\delta M_v$ represents the major contribution (approximately 80 %) to $\delta M$ . In the unstable case, there is a noticeable reduction in $\delta M$ due to a decrease in $\delta M_v$ . As mentioned previously, $\delta M_v$ results from two contributions: the bubble shape perturbation, $\delta M_{v,sh}$ , and the viscous stress perturbation, $\delta M_{v,\delta \tau }$ . Later, we analyse the critical role played by $\delta M_{v,\delta \tau }$ in the reduction of $\delta M$ .

First, we examine the components of the viscous stress perturbation on the bubble surface (figure 22 a) and their effect on the flow and lateral force exerted on the bubble (figure 22 b). Lastly, we analyse the $x$ component of the torque caused by the viscous stress perturbation, $\delta M_{v,\delta \tau }$ (figure 22 c). The tangential components balance the Marangoni stress resulting from the perturbation of the surfactant concentration. For both values of $\textit{Ga}$ , the tangential component $\delta \tau _{nt2}$ opposes the horizontal flow to the interface right side (the expanded and less loaded side). However, the component $\delta \tau _{nt1}$ switches direction and shows a smaller magnitude in the unstable case. As mentioned previously, $\delta \tau _{nn}\simeq 0$ for $\textit{Ga}=58$ , meaning that $\delta p\simeq \delta p_c$ . For $\textit{Ga}=53$ , there is significant normal stress in the southern hemisphere associated with the presence of secondary recirculation cells in the perturbation wake (figure 21 a). Figure 22(b) shows the lateral force $\delta F_{v,\delta \tau }$ resulting from the action of the viscous stress. For $\textit{Ga}=58$ , the effects of the two tangential components add up in the northern hemisphere and are counteracted in the southern one. The opposite occurs for $\textit{Ga}=53$ . In this case, the normal stress contributes to the lateral force in the southern hemisphere.

As anticipated previously, instability is attributed to a noticeable reduction in $\delta M$ due to the viscous contribution $\delta M_v$ . This contribution results from the bubble shape and viscous stress perturbations. Consider the $x$ component of the torque caused by the viscous stress perturbation, $\delta M_{v,\delta \tau }$ (figure 22 c). For the stable case $\textit{Ga}=53$ , the integral $\varPhi =\int _0^{\pi } \delta M_{v,\delta \tau }\, {\rm d}\alpha$ is negative, which implies that the viscous stress perturbation favours the anticlockwise overall rotation of the perturbed flow around the bubble (figure 21). In contrast, $\varPhi$ is positive for the unstable case, producing the reduction of the net torque. In both cases, $\delta M_{v,\delta \tau }$ opposes the recirculation cell at the rear (see, e.g. figure 21 a), limiting the magnitude of the velocity there (see, e.g. figure 21 b). However, the torque opposing the recirculation in the rear is significantly smaller in the unstable case (figure 22 c). This explains why the magnitude of the velocity field perturbation in the wake is larger than in the stable case (figure 19 c). In the unstable case, more energy is transferred to the recirculation in the rear, which can produce the tilting of the main toroidal vortex of the steady flow, leading to the oblique trajectory. A similar mechanism has been described by Natarajan & Acrivos (Reference Natarajan and Acrivos1993) and Ern et al. (Reference Ern, Risso, Fabre and Magnaudet2012) for the first bifurcation in the fixed-sphere case.

As mentioned previously, there are similarities between the instability mechanisms for a deformable bubble and a rigid sphere. A natural question is whether the interface deformation plays a relevant role in the bubble path instability. We do not have a definitive answer to this question. The torque $\delta M_{v,sh}$ associated with the perturbation of the interface location can be seen as the sum of two contributions: (i) that due to the bubble motion as a solid rigid and (ii) that due to the interface deformation shown in figure 17. We have verified that those two contributions are commensurate with each other, indicating that the interface deformation produces a significant effect on the torque. However, neither the torque $\delta M_{v,sh}$ nor its contributions change significantly when the flow becomes unstable (see figure S4 f). So, our qualitative explanation of the instability is not influenced by the interface deformation effect. Nevertheless, this deformation may affect the instability threshold value.