3.1. Bubble shapes

Figure 2 shows how the equilibrium shape of the bubble changes with the Bond number in three different fluids. As the Morton number keeps a constant value in each fluid, increasing $Bo$ (or $Ga$) in a given series amounts to increasing the bubble size; since $Bo\propto D^2$, the volume of the bubble grows as $Bo^{3/2}$. In each series, starting from nearly spherical shapes at low Bond number, bubbles take oblate spheroidal shapes when finite-$Bo$ effects become significant. The figure makes it clear that the properties of the carrying fluid have a major influence on the range of sizes in which this geometric approximation is valid: while the shape of a bubble with $Bo=3$ is still close to a spheroid in silicone oil DMS-T05 ($Mo= 6.2\times 10^{-7}$, see table 1 for the physical properties of the various fluids discussed below), a bubble with $Bo=0.5$ is seen to already exhibit a significant fore–aft asymmetry in water. Beyond that range, the front part of the bubble flattens increasingly as its size increases, while its rear part becomes more and more rounded. As indicated by the change in colour of the bubble contour in the figure, the transition between these two ‘regimes’ corresponds approximately to the onset of path instability. Although they do not rise in a straight line, bubbles with a mildly convex front and a pronounced rounded rear are routinely observed in experiments. In contrast, bubbles with an almost flat or even concave front such as those shown in the figure for $Bo\gtrsim 5$ in DMS-T05 or $Bo\gtrsim 0.7$ in water are not. The obvious reason is that such bubbles actually experience large shape oscillations, so that the present stationary solution makes little sense in these regimes.

Figure 2. Equilibrium shapes of bubbles with various Bond numbers in three different fluids; the Bond number is specified within each contour. Contours represent different bubbles; their centroids are arbitrarily shifted in the vertical direction for readability, making the local value of $z$ irrelevant. Blue and red contours refer to stable and unstable configurations, respectively. (a) Silicone oil DMS-T05 ($Mo=6.2\times 10^{-7}$); (b) silicone oil DMS-T02 ($Mo=1.6\times 10^{-8}$); (c) water at $20\,^\circ \text {C}$ ($Mo=2.54\times 10^{-11}$).

Table 1. Physical properties of some specific fluids considered in this study. Note that the surface tension of Galinstan may vary from $535$ to $718\ {\rm mN}\ {\rm m}^{-1}$, depending on the exact alloy composition.

3.2. Rise speed

Variations of the equilibrium rise speed with the bubble size and liquid properties are displayed in figure 3 in the form of a $Re(Bo)$ diagram. Comparison with experimental data reveals an excellent agreement in most fluids, especially in the most challenging case of water. A consistent $10\,\%$ over-estimate may be noticed in the least viscous silicone oil, DMS-T00, presumably because of some variation in the actual viscosity or surface tension of the sample used in the experiments with respect to the reference values provided by the manufacturer (Zenit & Magnaudet Reference Zenit and Magnaudet2008). In each series, the location of the incipient path instability detected in the experiments is specified with a bullet in the figure. The threshold Reynolds number is seen to vary by nearly one order of magnitude from water to the most viscous oil, DMS-T11, $9.4$ times more viscous than water. As the three intermediate series make clear, the predicted Reynolds number departs from that determined in the experiments beyond the threshold. This is no surprise, since part of the potential energy of zigzagging or spiralling bubbles is transferred to the radial and azimuthal fluid motions, yielding a decrease in the rise speed compared with the rectilinear regime.

Figure 3. Rise Reynolds number vs the Bond number in several fluids. Solid lines: present predictions; symbols: experimental data, with $\boxdot$ (blue): ultrapure water at $20\,^\circ {\rm C}$ ($Mo=2.54\times 10^{-11}$) (Duineveld Reference Duineveld1995), ${\times }$ (orange): silicone oil DMS-T00 ($Mo=1.8\times 10^{-10}$), ${\odot }$ (yellow): DMS-T02 ($Mo=1.6\times 10^{-8}$), (green): DMS-T05 ($Mo=6.2\times 10^{-7}$), and (purple): DMS-T11 ($Mo=9.9\times 10^{-6}$) (data for all four series taken from Zenit & Magnaudet Reference Zenit and Magnaudet2008). Red bullets mark the onset of path instability detected experimentally in each series.

Data reported in figure 3 may also be used to shed some light on the influence of the bubble shape and properties of the carrying fluid on the bubble rise speed. To this end, the above numerical predictions for the rise Reynolds number are replotted against the Galilei number in figure 4. This figure reveals the succession of two distinct behaviours in each fluid. First, at low enough $Ga$, all data collapse onto a master curve following the power law $Re\propto Ga^{5/3}$. This behaviour corresponds to the regime in which the bubble rise speed does not depend on surface tension, hence on the Morton number. In this regime, the bubble shape remains close to a sphere, i.e. its aspect ratio $\chi$ is close to unity. For low enough $We$, departures from sphericity are known to increase linearly with the Weber number according to the law $\chi =1+\frac {9}{64}We$ (Moore Reference Moore1965). Since $We=Re^2Mo^{1/3}Ga^{-2/3}$, the lower $Mo$ the wider the $Re$ range pertaining to this first regime, which extends up to $Re\approx 200$ in water. Bubbles become flatter as $We$, hence $Re$, increases, opposing a larger resistance to the fluid. This makes their rise speed grow more slowly with $Ga$ when the Weber number becomes of order unity. Since their frontal area depends on $We$, which itself depends on $Mo$, the Reynolds number is no longer independent of the Morton number in this second regime. To further examine the two behaviours, one has to refer to the force balance on the bubble which, in non-dimensional form, reads $C_D(Re)Re^2=\frac {4}{3}Ga^2$, with $C_D(Re)$ the usual drag coefficient. The drag coefficient of a spherical bubble is known to vary from $16Re^{-1}$ at low Reynolds number to $48Re^{-1}$ in the limit of very large Reynolds number (Batchelor Reference Batchelor1967). Hence, for nearly spherical bubbles, the force balance implies $k(Re)Re=\frac {1}{12}Ga^2$, with $k(Re\ll 1)=1$, $k(Re\gg 1)=3$. The slow increase of $k(Re)$ in between these two limits, which follows approximately the power law $k(Re)\sim Re^{1/5}$ up to Reynolds numbers of a few hundreds, is the reason why the Reynolds number grows slightly more slowly than $Ga^2$ in the first regime. Bubbles are significantly distorted in the second regime, which corresponds to the intermediate range $We={O}(1)$ ($We\approx 0.62$ at $Re=200$ in water). Assuming that their frontal area, which is proportional to $\chi ^{2/3}$, grows approximately linearly with the Weber number in that range yields $C_D(Re)\sim k(Re)Re^{-1}We\sim Re^{6/5}Mo^{1/3}Ga^{-2/3}$. The force balance then implies $Re\sim Ga^{5/6}Mo^{-5/48}$, a prediction supported by the numerical data for the various fluids in the $Ga$ range where path instability occurs.

Figure 4. Numerical predictions for the rise Reynolds number, plotted vs the Galilei number. The colour code and the meaning of the red bullets are similar to those of figure 3.

3.3. Flow structure

The flow structure in the base state reveals several interesting features. Figure 5 displays the azimuthal vorticity (left half of each panel) and pressure (right half) and distributions past a bubble rising in three different fluids in the $Bo$-range where the path instability threshold will later be shown to lie. With no surprise, the near-surface pressure distribution reaches its minimum very close to the equatorial plane, defined as the horizontal plane in which the longest axis of the bubble lies. This is a mere consequence of the increase of the fluid velocity along the interface from the apex of the bubble down to its equatorial plane. The norm of the azimuthal vorticity also reaches its maximum in that plane because this quantity is proportional to the product of the interface curvature in the vertical diametrical plane and the relative fluid velocity along the interface (Batchelor Reference Batchelor1967), both of which are maximum there. This vorticity results directly from the shear-free condition obeyed by the fluid at the interface and is responsible for the existence of a boundary layer that encapsulates the bubble and turns into a wake downstream of it (Moore Reference Moore1963, Reference Moore1965). The wake structure is found to depend dramatically on the Bond number and, for a given $Bo$, on the fluid properties (compare panels (a,f), both for $Bo=3$). In water, all streamlines are open for $Bo=0.4$ (panel g) and they remain so up to $Bo\approx 0.5$. Beyond this point, a standing eddy exists and grows sharply with the Bond number, its length becoming of the same order as that of the long bubble axis for $Bo=0.8$ (panel i). Qualitatively similar trends are observed in the more viscous fluids. However, for a given $Bo$, the normalized rise speed $U=u_b/(gD)^{1/2}$ is significantly smaller, and so is the Weber number $We\equiv \rho u_b^2D/\gamma =Bo\,U^2$. Indeed, keeping $Bo$ unchanged, the Weber numbers in two different fluids, 1 and 2, obey the relation $We_2/We_1=(Re_2/Re_1)^2(Mo_2/Mo_1)^{1/2}$. Hence, in the nearly five times more viscous silicone oil DMS-T05 for instance, the Weber number of a bubble with $Bo=0.5$ is eight times smaller than in water, leaving it nearly spherical. Consequently, compared with water, much larger $Bo$ are required for the tangential fluid velocity and interface curvature to be large enough for a standing eddy to develop at the back of the bubble. In the case of DMS-T05, the required conditions are achieved only for $Bo\gtrsim 2.4$.

Figure 5. Base flow past bubbles rising in three different liquids: (a–c) DMST05 ($Mo=6.2\times 10^{-7}$), with, from left to right, $Bo=3, 5$ and $7$; (d–f) DMST02 ($Mo=1.6\times 10^{-8}$), with $Bo=1$, $2$ and $3$; (g–i) water at $20\,^\circ {\rm C}$ ($Mo=2.54\times 10^{-11}$), with $Bo=0.4$, $0.6$ and $0.8$. The left and right halves of each panel display the azimuthal vorticity and pressure distributions, respectively. The thin lines are the streamlines in the reference frame rising with the bubble.

The above information regarding the existence of a standing eddy in the various liquids is summarized in figure 6. This plot shows the critical curve separating the upper region of the $(Bo,Ga)$ plane where a standing eddy exists, from the lower region in which no such structure takes place; the two insets display an example of each situation. Present predictions are seen to be in excellent agreement with those of the fully resolved axisymmetric simulations carried out with Gerris by Cano-Lozano et al. (Reference Cano-Lozano, Bohorquez and Martinez-Bazán2013). It will become clear in the next section that, in low-$Mo$ fluids, path instability arises for $(Ga,Bo)$ pairs located below this critical curve. This finding is of special importance in that it establishes that the existence of a standing eddy is not a prerequisite for the path of a freely rising bubble to become unstable, unlike the case of a fixed bubble (Tchoufag et al. Reference Tchoufag, Magnaudet and Fabre2013). More generally, it is well established that the initial ‘seed’ leading to global wake instability past bluff bodies held fixed in a uniform stream is the growth of disturbances in the core of the standing eddy (Chomaz Reference Chomaz2005). That path instability may happen in the absence of such a near-wake structure proves that, in the relevant $Mo$-range, the mechanism that governs this instability is not driven by wake instability per se. This confirms the findings of Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martinez-Bazán2016a) who established that, for $Mo\lesssim 10^{-9}$, the path of bubbles with a frozen but realistic shape becomes unstable while their wake is still stable, the low-frequency ‘aerodynamic’ mode being destabilized first.

Figure 6. Critical curve corresponding to the onset of a recirculating region at the back of the bubble in the $(Bo,Ga)$ plane. Yellow line: present results; black dashed line: predictions of Cano-Lozano, Bohorquez & Martinez-Bazán (Reference Cano-Lozano, Bohorquez and Martinez-Bazán2013). The left and right insets show some streamlines around a bubble with $Bo=0.3$ in water and a bubble with $Bo=6.5$ in DMS-T05, respectively. The thin dotted lines correspond to constant values of the Morton number, i.e. to a given liquid; $Mo$ values are specified along the upper and right sides of the figure.

4.1. Critical bubble size

Using the procedure described in § 2, we selected a set of $Mo$ values, some of which correspond to specific liquids. For each of these values, we increased the bubble size, i.e. the Bond and Galilei numbers, until the real part of one of the eigenvalues associated with the non-axisymmetric modes $m=\pm 1$ changes sign, which determines the threshold of path instability. The resulting neutral curve, obtained by linking these thresholds, is shown in figure 7. Together with figure 8, this is presumably the most important result of the present study with respect to the original physical problem. For each fluid characterized by a Morton number in the range $10^{-13}\unicode{x2013}10^{-5}$, this curve readily provides the size of the smallest bubble whose vertical path becomes linearly unstable. For instance, the path of a bubble rising in water at a temperature of $20\,^\circ {\rm C}$ becomes unstable at $Bo=0.463$ and $Ga=250$, yielding a critical diameter $D\approx 1.854$ mm (Bonnefis Reference Bonnefis2019; Bonnefis et al. Reference Bonnefis, Fabre and Magnaudet2023). Similarly, the critical Bond numbers for Galinstan ($Mo=1.4\times 10^{-13}$) and DMS-T11 ($Mo=9.9\times 10^{-6}$) are $0.19$ and $6.85$, respectively. Therefore, the critical bubble sizes for these two fluids located close to the bounds of the $Mo$-interval considered here are $1.48$ mm and $3.875$ mm, respectively.

Figure 7. Neutral curve in the $(Bo,Ga)$ plane, obtained by connecting the (orange) symbols at which the threshold was determined. The pale and dark grey lines are the two branches of the neutral curve obtained by considering the same base state but keeping the bubble shape frozen in the stability analysis; the thin yellow dashed line is the critical curve of figure 6 beyond which a standing eddy exists at the back of the bubble in the base state. The thin dotted lines are iso-Morton-number lines corresponding to a given liquid; $Mo$ values are specified along the upper and right sides of the figure. In a given fluid, open (respectively closed) symbols indicate stable (respectively unstable) paths observed in the simulations of Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazán, Magnaudet and Tchoufag2016b) (, $\blacklozenge$), in the experiments of De Vries (Reference De Vries2002) (${\odot }$, ${\bullet }$, purple) and Duineveld (Reference Duineveld1995) (${\boxdot }$, ${\blacksquare }$ purple) with ultrapure water (performed at temperatures of $28\,^\circ {\rm C}$ and $20\,^\circ {\rm C}$, respectively), and in those of Zenit & Magnaudet (Reference Zenit and Magnaudet2008) (,${\blacktriangle }$ purple) and Sato (Reference Sato2009) (,${\blacktriangledown }$ purple) with various silicone oils.

Figure 8. Variation of the reduced frequency $St=\lambda _i D/(2{\rm \pi} u_b)$ at the threshold vs the Morton number. The yellow curve was obtained by connecting the (orange) symbols at which the frequency was determined. The pale and dark grey lines correspond to the predictions along the two branches of the neutral curve predicted by the FSA using the same base state (for readability, only five of the computed points are highlighted with a marker on the dark grey line). The vertical arrow at $Mo=2.7\times 10^{-7}$ indicates the frequency jump associated with the switching from the stationary mode to the oscillatory mode predicted by FSA. $\blacklozenge$ (black): numerical predictions of Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazán, Magnaudet and Tchoufag2016b); ${\blacksquare }$ (purple): experimental data of Duineveld (Reference Duineveld1994) for $Bo=0.54$ in ultrapure water; ${\blacktriangle }$ (purple): experimental data of Zenit & Magnaudet (Reference Zenit and Magnaudet2009) for $Bo=3.92$ in DMS-T05.

In figure 7 we also report some experimental data obtained under controlled conditions, ensuring that the carrying fluid obeys a shear-free condition at the gas–liquid interface. The reference data of Duineveld (Reference Duineveld1995) in ultrapure water at $20\,^\circ {\rm C}$ yield a critical bubble diameter close to $1.82$ mm, which differs by less than $2\,\%$ from the present prediction. De Vries (Reference De Vries2002) also determined the onset of path instability in a tank of ultrapure water with an average temperature of $28\,^\circ {\rm C}$ ($Mo=1.13\times 10^{-11}$). However, a $1.1\,^\circ {\rm C}\ {\rm m}^{-1}$ stabilizing temperature gradient was established to visualize the wake thanks to optical index gradients. The influence of this stratification on path instability is unknown and might be at the origin of the $8\,\%$ difference between the experimentally determined critical bubble size and the present prediction. Comparisons with data obtained with four different silicone oils (Zenit & Magnaudet Reference Zenit and Magnaudet2008) and with three other silicone oils in which photochromic dye was added to visualize the wake (Sato Reference Sato2009) prove satisfactory in the sense that the neutral curve lies generally within the ($Bo,Ga$)-interval whose lower (respectively upper) bound corresponds to the largest (respectively smallest) bubble for which a stable (respectively unstable) path could be identified experimentally. The only two exceptions correspond to a probable outlier in Sato's data for $Mo=2.35\times 10^{-7}$, and to DMS-T11 (the most viscous oil) in which the path of a bubble $4\,\%$ smaller than the critical size predicted here was found to be unstable. Numerical data obtained by Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazán, Magnaudet and Tchoufag2016b) with Gerris in the transition region (but not necessarily in the close vicinity of the threshold) are also reported. These data are seen to be in good agreement with present predictions for fluids with $Mo\gtrsim 10^{-8}$. In contrast, as the authors anticipated in their conclusions, these simulations underestimate the critical bubble size in fluids having a lower Morton number. The lower $Mo$ the larger the underestimate, which makes the critical diameter in water under-predicted by 15 %–20 %. Since these simulations made use of an adaptive grid with $128$ cells per bubble diameter on both sides of the interface, this significant underestimate helps appreciate how demanding fully resolved simulations of bubbly flows in low-viscosity fluids are.

In figure 7 we also reported the critical curve of figure 6 corresponding to the onset of a recirculating region at the back of the bubble in the base flow (yellow dashed line). The yellow solid and dashed lines are seen to cross each other at $Mo\approx 10^{-10}$. For lower $Mo$, the critical bubble size at which the path destabilizes first is such that no standing eddy exists in the base state. For the reasons discussed in the previous section, this leads to the conclusion that wake instability cannot be responsible for path instability in such low-$Mo$ fluids. The pale and dark grey lines in the figure also provide interesting insight into the role of time-dependent bubble deformations. To obtain the neutral curve resulting from these two lines, we considered the same base state as above, and imported the corresponding interface shape and rise speed in the stability code used by Tchoufag et al. (Reference Tchoufag, Magnaudet and Fabre2014b) and Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martinez-Bazán2016a). In this code, the interface shape is kept frozen, i.e. the bubble can only move as a rigid body. A similar neutral curve was computed by Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martinez-Bazán2016a) (black line in their figure 7), using a base state obtained with Gerris, instead of the present global Newton method. Their neutral curve is similar to that based on the two grey lines in figure 7 but the critical $Ga$ they found at a given $Bo$ stands consistently below that determined here, with differences ranging from $5\,\%$ for $Bo\geqslant 5$ to $15\,\%$ for $Bo=0.5$. Again, the reasons for these differences are to be found in the limited spatial resolution of the Gerris-based computations in high-$Ga$ low-$Bo$ configurations.

Since the yellow and grey lines in figure 7 are based on strictly identical base states, time-dependent deformations not accounted for in the stability analysis used to produce the second of them are entirely responsible for the differences observed in the fate of the imposed disturbances. As the comparison of the thresholds predicted by the two approaches in a given fluid evidences, deformations always lower the threshold, i.e. they promote path instability. However, for fluids with Morton numbers less than ${\approx }5\times 10^{-7}$, the lower $Mo$ the weaker this influence. For instance, the ‘frozen-shape’ approximation, hereinafter referred to as FSA, over-predicts the critical bubble size by nearly $23\,\%$ for $Mo=10^{-7}$, but this overestimate is reduced by a factor of three in the case of Galinstan. The FSA actually predicts that the neutral curve involves two distinct modes intersecting at $Mo\approx 2.8\times 10^{-7}$. In the narrow range $4\times 10^{-7}\lesssim Mo\lesssim 6.5\times 10^{-7}$, the mode that becomes most unstable beyond this intersection (dark grey line) exhibits a S shape, implying a destabilization–restabilization scenario. In contrast, the neutral curve resulting from the stability analysis carried out with freely deforming bubbles remains single valued throughout the whole range of Morton number. For larger $Mo$, the FSA prediction based on this second mode is also quite good, with for instance a $16\,\%$ over-prediction of the critical bubble size in DMS-T11.

At this point it is also worth mentioning that Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martinez-Bazán2016a) established that the neutral curve predicted by FSA and that corresponding to the onset of ‘pure’ wake instability for a bubble with the same frozen shape held fixed in a uniform stream coincide in the range $10^{-9}\lesssim Mo\lesssim 6\times 10^{-7}$. Outside of this range, the threshold values of the Bond and Galilei numbers corresponding to the onset of wake instability in a given fluid are consistently larger than those at which FSA predicts the occurrence of path instability. For $Mo\lesssim 10^{-9}$, the difference between the two thresholds increases gradually as $Mo$ is decreased. For $Mo\gtrsim 6\times 10^{-7}$, the threshold of wake instability coincides with the upper neutral branch captured by FSA (here the pale grey line). In other words, depending on the Morton number, the threshold corresponding to the onset of wake instability past a fixed bubble is either identical to or larger than that of the path instability predicted by FSA. These findings, combined with the fact that figure 7 reveals that FSA overestimates the threshold of path instability whatever $Mo$, establish that the wake of a fixed bubble with a frozen shape is always stable at the actual onset of path instability. Hence, path instability cannot be explained on the sole basis of the mechanism responsible for wake instability.

4.2. Frequency at threshold

The most spectacular qualitative difference between the predictions of the present approach and those of FSA is seen in figure 8. This figure displays the Strouhal number, or reduced frequency, corresponding to the first unstable mode at the onset of the instability, computed as $St=\lambda _iD/(2{\rm \pi} u_b)$, with $\lambda _i$ the imaginary part of the unstable eigenvalue at the threshold. As the yellow curve evidences, the present approach predicts that the most unstable mode of the system is always oscillatory, in line with experimental observations. The reduced frequency increases continuously with the Morton number, ranging from $St\approx 0.02$ for $Mo=10^{-13}$ to $St\approx 0.12$ for $Mo=10^{-5}$. The FSA predicts that the first non-vertical path of the bubble is associated with ‘low-frequency’ oscillations in low-$Mo$ liquids ($Mo\lesssim 7\times 10^{-11}$) and with ‘high-frequency’ oscillations in liquids with $Mo\gtrsim 3\times 10^{-7}$, in agreement with the conclusions of Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martinez-Bazán2016a). Present findings with deformable bubbles are qualitatively consistent with these predictions in both ranges. The reduced frequencies determined in both approaches for low-$Mo$ liquids are even in fairly good quantitative agreement. This confirms that deformations do not bring major changes in the dynamics of isolated bubbles rising in such liquids, which extends to all low-$Mo$ fluids the conclusions drawn by Bonnefis et al. (Reference Bonnefis, Fabre and Magnaudet2023) in the specific case of water. More surprisingly, the agreement is still reasonable in fluids with $Mo\gtrsim 3\times 10^{-7}$, where the FSA neutral curve (now corresponding to the dark grey line) underestimates the reduced frequency by only roughly $20\,\%$. In contrast, the two approaches dramatically disagree in the intermediate $Mo$-range, where FSA predicts that the most unstable mode (pale grey line) is stationary. Hence, one has to conclude that, in the linear framework adopted here, time-dependent deformations promote the instability of the oscillatory mode and succeed in making it more unstable than the stationary mode in this intermediate range, while only the latter mode may become unstable if these deformations are ignored. Some experimental and numerical data are also reported in figure 8. Most of them were obtained significantly above the threshold, which explains why the corresponding frequencies lie above present predictions (for instance, the bubble size is $8\,\%$ beyond the threshold value in the two experimental determinations). Only the two numerical determinations of Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazán, Magnaudet and Tchoufag2016b) in silicone oils DMS-T02 ($Mo=1.6\times 10^{-8}$) and DMS-T11 ($Mo=9.9\times 10^{-6}$) correspond to near-threshold conditions according to figure 7. The corresponding two $St$ values are seen to agree very well with present predictions.

Having determined the normalized rise speed in the base state, we can re-plot the neutral curve of figure 7 in two different forms, to make the variations of the critical Reynolds number, $Re=Ga\,U$, and Weber number, $We=Bo\,U^2$, with the liquid properties apparent. Figure 9(a) shows that the critical Reynolds number decreases sharply as $Mo$ increases, from ${\approx }2070$ for $Mo=10^{-13}$ to $70$ for $Mo=10^{-5}$, via $Re=668$ for water at $20\,^\circ {\rm C}$. Conversely, the critical Weber number (figure 7b) is seen to increase gradually from $3.20$ for $Mo=10^{-13}$ to $6.36$ for $Mo=10^{-5}$, via $We=3.30$ for water. The above values for water are to be compared with those reported by Duineveld (Reference Duineveld1995), namely $Re=658$, $We=3.275$, from which they only differ by $1.5\,\%$ and $0.75\,\%$, respectively.

Figure 9. Influence of the fluid properties on (a) the critical Reynolds number and (b) the critical Weber number.

5.1. Silicone oil DMS-T05 ($Mo=6.2\times 10^{-7}$)

5.1.1. First unstable eigenvalues: thresholds and nature of the associated modes

Figure 10 displays the variations of the first five eigenvalues whose real part becomes positive as the Bond number, hence the size, of a bubble rising in silicone oil DMS-T05 is increased. All but one of these eigenvalues have a non-zero imaginary part, indicating that the corresponding instabilities arise through Hopf bifurcations. The threshold of path instability is encountered at $Bo=3.4$. This first unstable mode, associated with azimuthal wavenumbers $|m|=1$, has a dimensionless frequency $\lambda _i/2{\rm \pi} \approx 0.115$ at the threshold (throughout this section, eigenvalues are normalized using the gravitational time scale $(D/g)^{1/2}$). As panel (c) suggests, this mode is essentially associated with a lateral drift of the bubble centroid accompanied by some rigid-body rotation. Time-dependent deformations are small: it will be shown in § 6 that their maximum magnitude is barely $3\,\%$ of that of the horizontal displacement of the bubble centroid. For this reason, we refer to this mode as ‘rigid’.

Figure 10. Variation of the first five unstable eigenvalues with the Bond number for bubbles rising in DMS-T05 ($Mo=6.2\times 10^{-7}$). (a) Growth rate ($\lambda _r$), normalized by the gravitational time $(D/g)^{1/2}$; (b) radian frequency ($\lambda _i$) normalized similarly; (c) equilibrium shapes at the threshold (labels refer to the transition points marked with bullets in a). In (a,b), dashed (respectively solid) lines represent the stable (respectively unstable) part of the various branches. In (c), the bubble contour in the vertical diametrical plane is shown at the threshold of each of the five successive modes in the base state (black solid line) and, with an arbitrary amplitude of the disturbance, for $t=0$ (red dashed line) and $t={\rm \pi} /(2\lambda _i)$ (blue dashed line); bubble centroids are arbitrarily shifted in the $z$ direction for readability.

We continue to track this mode and other possible ones beyond this primary threshold by computing the base flow, still assumed axisymmetric, at the relevant Bond number. The growth rate of the rigid oscillatory mode continues to increase with $Bo$, and a second mode becomes unstable at $Bo\approx 5.5$ (point $B$ in figure 10a). This mode is axisymmetric and oscillatory, with a dimensionless frequency approximately $5.5$ times higher than that of the rigid mode. Examination of the corresponding successive contours in panel (c) reveals that this mode is associated with shape oscillations known as $l=2,m=0$, or $(2,0)$, in the terminology of spherical harmonics, i.e. the associated interface displacements leave the interface position unchanged at two angular positions located on both sides of the equatorial plane. We shall show in § 6 that these deformations have a magnitude of the same order as the vertical displacement of the bubble centroid with respect to its reference position in the base state. For these reasons, we refer to this mode and those sharing similar properties as ‘shape’ modes. By further increasing $Bo$, a third unstable mode is encountered at point $C$ ($Bo\approx 5.9$). Its growth rate increases strongly beyond the threshold so that this mode would quickly become dominant, were its amplitude similar to that of the primary oscillatory mode. Panel (b) reveals that this mode is stationary and, since its azimuthal wavenumber is $|m|=1$, it is also asymmetric. These two properties imply that this mode induces an inclined path of the bubble centroid. Moreover, panel (c) indicates that the corresponding interface deformations are weak compared with those associated with the previous ‘shape’ mode, so that this stationary mode can be considered ‘rigid’, just like the primary mode. Therefore, two rigid modes coexist in the system beyond $Bo\approx 5.9$, but one is oscillatory while the other is stationary. To the best of our knowledge, steady oblique bubble paths have never been reported in experiments, suggesting that this second rigid mode never becomes dominant. Nonlinearities not being taken into account here, i.e. the fact that the actual base flow at point $C$ is no longer axisymmetric due to the nonlinear corrections brought by the primary oscillatory mode and is presumably even not stationary owing to the shape oscillations associated with the $(2,0)$ shape mode, are likely the reason for this.

A fourth mode becomes unstable at point $D$ ($Bo\approx 6.6$), and restabilizes itself within the narrow interval $8.65\lesssim Bo\lesssim 9.45$ before becoming unstable again at larger $Bo$. This mode has much in common with the $(2,0)$ shape mode, except that it is not axisymmetric. Rather, it is associated with azimuthal wavenumbers $|m|=2$ and, in the $Bo$-range under consideration, its frequency is nearly half that of the $(2,0)$ mode. That the frequencies of two modes with the same $l$ but different $m$ differ in this context is no surprise, since spherical harmonics are no longer the eigenmodes with which a bubble oscillates when its undisturbed shape is not spherical (Meiron Reference Meiron1989). Here, the Bond number is large, and it is only in the small-$Bo$ limit that the two frequencies are expected to get close to each other. Due to the symmetries preserved by angular variations of the form ${\rm e}^{\pm 2{\rm i}\theta }$, this $(2,2)$ mode does not induce any lateral displacement of the bubble. Last, at point $E$ ($Bo\approx 9.45$), a third shape mode with a frequency close to that of the $(2,2)$ mode becomes unstable. Deformations associated with this mode share the same characteristics as those induced by modes $(2,0)$ and $(2,2)$ on both sides of the bubble equatorial plane. In contrast, as panel (c) makes clear, the two halves of the bubble deform in an asymmetric manner at every instant of time in the vertical diametrical plane, one half becoming more pointed while the other becomes more rounded. These characteristics identify this mode as the $(2,1)$ shape mode. Unlike those associated with $m=0$ and $|m|=2$, this mode may contribute to the lateral drift of the bubble at large enough Bond numbers, owing to its azimuthal asymmetry. It must be stressed that the full sequence described above is kept unchanged when the Morton number is increased up to $10^{-5}$.

5.1.2. Spatial structure of the first unstable mode

Figure 11 shows the spatial structure of the first unstable mode slightly above the threshold. The real and imaginary parts of a given mode may be thought of as the associated disturbance at two different instants of time a quarter of a period apart, i.e. shifted by an interval $\tau =({{\rm \pi} }/{2})\lambda _i^{-1}$. Therefore, if the bubble is zigzagging in a vertical plane (as it does if the two disturbances associated with modes $m=+1$ and $m=-1$ have the same amplitude), the real part corresponds to the time instant by which the bubble reaches its maximal lateral excursion (red contours in figures 11(a) and 10(c)), while the imaginary part corresponds to that by which it crosses the midline of its path (blue contours in figures 11(b) and 10(c)). The streamlines pattern evidences the antisymmetric nature of the disturbance flow. An alternation of positive and negative vorticity zones may be noticed in the wake, both in every horizontal plane (varying $r$) and at every radial position from the axis of the base flow (varying $z$). This pattern is indicative of a vortex shedding process. A specific structure in which fluid particles rotate clockwise on both sides of the axis of the base flow is visible near the bottom of panel (b); in a three-dimensional representation, this structure would look like a pair of crescents connected through their horns in the vertical plane $\theta =\pm {\rm \pi}/2$, with fluid particles rotating in the same direction in every azimuthal plane. As panels (c,d) make clear, this structure is actually the first of a series that develops further downstream in the wake, with alternating positive and negative directions of rotation.

Figure 11. Structure of the first unstable mode past a bubble rising in DMS-T05 at $Bo=3.7$. (a,c) Real part of the mode near the bubble and further downstream in the wake, respectively; (b,d) same for the imaginary part. The left and right halves in (a,b) display the pressure and azimuthal vorticity iso-levels, respectively; those in (c,d) display the azimuthal velocity and vorticity iso-levels, respectively; some streamlines defined in the reference frame of the base configuration are also shown. Black, red and blue bubble contours correspond to the base state, and the real and imaginary parts of the interface disturbance ($\hat \eta$), respectively.

5.2. Silicone oil DMS-T02 ($Mo=1.6\times 10^{-8}$)

We now move to silicone oil DMS-T02 which is $2.6$ times less viscous than DMS-T05 and has almost the same surface tension. Figure 12 shows how the first five eigenvalues yielding unstable modes vary with the Bond number in this case. The beginning of the sequence is similar to that described for DMS-T05. That is, path instability first arises through an oscillatory ‘rigid’ mode that becomes unstable at $Bo=1.56$ (point $A$ in panel a). Then the shape mode $(2,0)$ becomes unstable at point $B$ ($Bo=2.2$), followed by the stationary rigid mode at point $C$ ($Bo=2.43$). Further increasing $Bo$, the $(2,1)$ and $(2,2)$ shape modes become unstable at points $E$ ($Bo=3.89$) and $D$ ($Bo=4.12$), respectively. Note that, compared with the case of DMS-T05, these last two modes emerge in reverse order.

Figure 12. Same as figure 10 for bubbles rising in DMS-T02 ($Mo=1.6\times 10^{-8}$).

Nevertheless, the main originality of the present bifurcation diagram compared with that of DMS-T05 is that the growth rate of the primary mode starts to decrease slightly beyond point $B$, and becomes even very weakly negative within a tiny interval around $Bo\approx 2.55$, before turning positive again and growing continuously up to the upper bound of the explored Bond number range. This non-monotonic behaviour is to be related to the fact that DMS-T02 stands in the middle of the $Mo$-range where FSA predicts the emergence of path instability through a stationary bifurcation instead of a Hopf bifurcation, being unable to detect that the oscillatory rigid mode becomes actually unstable first. This failure is presumably closely connected with the fact that, although interface deformations succeed in making the oscillatory mode unstable first, the growth rate of this mode remains quite low throughout the interval bounded by points $A$ and $C$. We shall come back to this point in § 7.

5.3. Water at $20\,^\circ {\rm C}$ ($Mo=2.54\times 10^{-11}$)

5.3.1. First unstable eigenvalues: thresholds and nature of the associated modes

We finally consider the case of bubbles rising in pure water maintained at a temperature of $20\,^\circ {\rm C}$. Variations of the eigenvalues leading to the first five unstable modes are shown in figure 13. Again, the sequence begins with the emergence of an oscillatory rigid mode that becomes unstable at $Bo=0.463$ with a dimensionless frequency $\lambda _i/2{\rm \pi} =0.088$, followed by the $(2,0)$ shape mode at $Bo=0.525$. However, the next steps differ deeply from those encountered in the previous two cases. The primary mode grows continuously until $Bo\approx 0.64$. There, this mode splits into two separate branches, both of which are stationary. Therefore, at the point (marked with a square in the figure) where the three branches meet, the pair of complex eigenvalues associated with the oscillatory mode turns into a pair of real eigenvalues. Such a codimension 2 singular point is referred to as ‘exceptional’ in the theory of non-Hermitian Hamiltonian systems (Bergholtz, Budich & Kuns Reference Bergholtz, Budich and Kuns2021). In the present case, this exceptional point, at which the two eigenvalues change their nature, corresponds to a Takens–Bogdanov bifurcation with ${O}(2)$-symmetry (Dangelmayr & Knobloch Reference Dangelmayr and Knobloch1987). Beyond this point, the upper stationary branch exhibits a rapid and continuous growth. Conversely, the growth rate of the lower branch decreases sharply with the distance to the exceptional point, until it becomes negative at $Bo=0.685$ (point $C$ in the figure). Similar to the cases of DMS-T02 and DMS-T05, and undoubtedly for the same reasons, inclined paths corresponding to the above stationary modes have not been observed in experiments carried out in water, irrespective of the bubble size.

Figure 13. Same as figure 10 for bubbles rising in water at $20\,^\circ {\rm C}$ ($Mo=2.54\times 10^{-11}$). In (a), the square indicates the codimension-two ‘exceptional’ point at which the pair of complex conjugate eigenvalues associated with the first unstable mode turns into a pair of real eigenvalues.

Next, figure 13 indicates that the mode $(2,2)$ becomes unstable at $Bo=0.68$ (point $D$ in the figure) and then grows slowly with the Bond number. Then, a fourth mode becomes unstable at $Bo=0.809$ (point $A'$) with a dimensionless frequency $\lambda _i/2{\rm \pi} = 0.162$. This mode (identified as rigid osc. II in the legend of figure 13a) is associated with an azimuthal wavenumber $|m|=1$, exhibits low-frequency oscillations (panel b), and is not accompanied by significant shape oscillations (panel c). These characteristics point to another rigid mode. Hence, similar to the right part of figures 10 and 12, two rigid modes coexist in the case of water beyond point $A'$. In all three cases, the most amplified of them is stationary (yellow curve in figure 13) while the other (red curve) exhibits slow oscillations. Finally, figure 13 indicates that a fifth eigenvalue crosses the real axis at $Bo\approx 0.88$ (point $E$). As panels (b,c) make clear, this eigenvalue corresponds to the asymmetric shape mode $(2,1)$ already encountered with DMS-T02 and DMS-T05. One may notice that the mode $(2,1)$ becomes unstable at a larger $Bo$ than the mode $(2,2)$ in the cases of water and DMS-T05, while the reverse is true for DMS-T02. It is also worth noting that, among the three shape modes revealed by figure 13, the mode $(2,2)$ is the one with the lowest frequency in the range of Bond numbers relevant here, say $Bo\gtrsim 0.4$, followed by the mode $(2,0)$ and then the mode $(2,1)$. This is in line with the conclusions obtained by Meiron (Reference Meiron1989) assuming an inviscid potential flow. In contrast, figures 10 and 12 show that modes $(2,1)$ and $(2,2)$ have nearly equal frequencies, both significantly lower than that of the mode $(2,0)$, in DMS-T05 and DMS-T02. This suggests that vortical effects not accounted for in the above reference significantly modify the dynamics of shape oscillations of strongly deformed bubbles, even in low-viscosity fluids. Indeed, although the viscosity of DMS-T02 is less than twice that of water, its surface tension is four times lower, so that the relevant Bond numbers in DMS-T02 are typically larger than $3.0$, while those in water are less than $1.0$.

To finish with figure 13, it is important to stress that the succession of modes it reveals is the canonical sequence encountered throughout the low-$Mo$-range, say for $Mo\lesssim 10^{-9}$. For instance, the same first five modes become unstable in the same order in the case of DMS-T00 ($Mo=1.8\times 10^{-10}$), with the low-frequency rigid mode splitting into two stationary modes at an exceptional point, just as in figure 13. Only the thresholds and frequencies differ, with points $A$–$B$ then located at $Bo=0.66$ and $0.805$, followed by points $C$ and $D$ almost coincident at $Bo=0.98$, and finally points $A'$ and $E$ located at $Bo=1.135$ and $1.42$, respectively.

5.3.2. Spatial structure of rigid and shape modes

Figure 14 displays the structure of the three successive rigid modes encountered in the sequence analysed above. Panels (a,b) look qualitatively similar to their counterparts in figure 11. However, the ‘double-crescent’ structure identified in the bottom part of the latter is not present here. The reason is that the Strouhal number of the primary low-frequency mode is approximately three times lower in water (see figure 8). For this reason, the vertical distance separating vortices shed in the wake is much larger, so that the double-crescent structure closest to the bubble is located further downstream. Panels (c,e) compare the wake structure of the primary and secondary oscillatory rigid modes. They evidence the higher vortex shedding frequency of the latter, in line with the nearly twice as large dimensionless frequency noticed at point $A'$ in figure 13 compared with that at point $A$. In the case of the stationary mode displayed in panel (d), the azimuthal vorticity and velocity disturbances are seen to keep a constant sign all along the wake (the real part of the azimuthal velocity is null in this case, which is why the imaginary part is shown in the figure). This is the generic hallmark of wakes associated with stationary inclined paths of axisymmetric bodies; e.g. Tchoufag et al. (Reference Tchoufag, Fabre and Magnaudet2014a, Reference Tchoufag, Magnaudet and Fabreb). No vortex shedding takes place in that situation, the inclination of the body path resulting from the constant lift force generated by a pair of semi-infinite counter-rotating streamwise vortices.

Figure 14. Structure of the rigid modes past a bubble rising in water. (a,b) Primary low-frequency mode at $Bo=0.48$ in the bubble vicinity, with the real and imaginary parts shown in (a,b), respectively. The left and right halves of the two panels display the pressure and azimuthal vorticity iso-levels, respectively; some streamlines defined in the reference frame of the base configuration are also shown. Black, red and blue bubble contours correspond to the base state, and the real and imaginary parts of the interface disturbance, respectively. (c–e) Wake structure of the successive rigid modes. (c) Primary low-frequency mode at $Bo=0.48$; (d) most amplified stationary mode at $Bo=0.75$; (e) secondary oscillating mode at $Bo=0.85$. The right half of each panel displays the real part of the azimuthal vorticity iso-levels; the left half of (c,e) (respectively d) displays the real (respectively imaginary) part of the azimuthal velocity iso-levels.

Figure 15 reveals the structure of the $(2,0)$ and $(2,1)$ shape modes at $Bo=0.54$ and $Bo=0.9$, respectively. The axial symmetry of the flow field is obvious in panels (a,b), with a succession of toroidal eddies of alternate signs shed downstream of the bubble. The red and blue contours evidence the fact that the bubble alternatively ‘inflates’ and ‘deflates’ along its axis while simultaneously shortening or lengthening in its equatorial plane. The streamlines in panels (c,d) make the antisymmetry of the mode $(2,1)$ just as obvious. The wake is dominated by a series of double-crescent structures qualitatively similar to those identified in figure 11. The wavelength separating two consecutive structures is slightly shorter than that of the toroidal eddies in panels (a,b), although the dimensionless frequency of the mode $(2,1)$ at $Bo=0.9$ is approximately $30\,\%$ lower than that of the mode $(2,0)$ at $Bo=0.54$ according to figure 13(b). However, these structures are advected downstream by the base flow, so that their wavelength is determined by the Strouhal number $St=\lambda _i/(2{\rm \pi} U)$, not directly by $\lambda _i$. Therefore, what the comparison of the wavelengths for the two modes indicates is that the normalized rise speed $U$ decreases by slightly more than $30\,\%$ from $Bo=0.54$ to $Bo=0.9$, owing to the progressive flattening of the bubble (in the regime where bubbles keep an approximate oblate spheroidal shape, $U$ is known to vary approximately as $Bo^{-1/2}$ for $Bo\ll 1$ and $Bo^{-1}$ for $Bo\gg 1$; see equation (7-3) and the associated figure in Clift, Grace & Weber Reference Clift, Grace and Weber1978).

Figure 15. Structure of the first axisymmetric and asymmetric shape modes past a bubble rising in water. (a,b) Axisymmetric $(2,0)$ mode at $Bo=0.54$, with the real and imaginary parts shown in (a,b), respectively. (c,d) Same for the asymmetric $(2,1)$ mode at $Bo=0.90$. For caption see figures 11 and 14.

6.1. Decomposition technique

In the previous section, we identified ‘rigid’ and ‘shape’ modes qualitatively, stating that the former leave the bubble shape almost unchanged while the latter are associated with significant time-dependent shape oscillations. However, this distinction deserves a more quantitative basis. Moreover, determining how large the deformations at the threshold of path instability are, compared with the lateral drift of the bubble, is key to clarifying the relative role of the various possible causes of instability present in the system. In order to quantify time-dependent deformations, we define a $(x,y,z)$ Cartesian coordinate system, the origin of which stands at the centroid of the bubble in the base state. We assume that the bubble performs planar zigzags in the vertical $(x,z)$ plane corresponding to the plane $\theta =(0,{\rm \pi} )$ in the cylindrical $(r,\theta,z)$ system defined in § 2. At any position ${\boldsymbol {x}}_0$ on the undisturbed bubble–fluid interface, the normal displacement of the interface, $\hat {\eta }({\boldsymbol {x}}_0)$, may be decomposed into a rigid motion inducing a lateral (respectively vertical) displacement $\hat {T}_x$ (respectively $\hat {T}_z$) about the $x$ (respectively $z$) axis and an inclination $\hat {\psi }$ resulting from a rotation about the $y$ axis, augmented by a local volume-preserving deformation $\hat {\zeta }({\boldsymbol {x}}_0)$ in the normal direction. Projecting this decomposition onto the local normal unit vector ${\boldsymbol {n}}_0$, one has

(6.1)\begin{equation} \hat{\eta}({\boldsymbol{x}}_0)=\hat{T}_x({\boldsymbol{e}}_x\boldsymbol{\cdot}{\boldsymbol{n}}_0)+ \hat{T}_z({\boldsymbol{e}}_z\boldsymbol{\cdot}{\boldsymbol{n}}_0)+\hat{\psi}({\boldsymbol{x}}_0\times{\boldsymbol{n}}_0) \boldsymbol{\cdot}{\boldsymbol{e}}_y+\hat{\zeta}({\boldsymbol{x}}_0), \end{equation}

with ${\boldsymbol {e}}_x,{\boldsymbol {e}}_y,{\boldsymbol {e}}_z$ the unit vectors in the $x,y,z$ directions, respectively. Odd modes $m=\pm 1$ do not induce any translation in the vertical direction, so that $\hat {T}_z=0$ for such modes. Conversely, the lateral displacement $\hat {T}_x$ and the inclination $\hat {\psi }$ of the bubble remain null for even modes $m=-2,0,2$. To determine the non-zero displacement and inclination associated with a given mode, we make use of a least-squares fitting technique. That is, we define the local functional $\mathcal {F}_{\hat {T}_x,\hat {T}_z,\hat {\psi }}({\boldsymbol {x}}_0)\equiv \{\hat {\eta }-[\hat {T}_x ({\boldsymbol {e}}_x\boldsymbol {\cdot }{\boldsymbol {n}}_0)+\hat {T}_z({\boldsymbol {e}}_z \boldsymbol {\cdot }{\boldsymbol {n}}_0)+\hat {\psi }({\boldsymbol {x}}_0\times {\boldsymbol {n}}_0) \boldsymbol {\cdot }{\boldsymbol {e}}_y]\}^2$, integrate it over the entire bubble contour and minimize the corresponding residue with respect to the non-zero components of the triplet $(\hat {T}_x,\hat {T}_z,\hat \psi )$. Given the linear framework of the present study, the magnitude of the displacements is arbitrary, and only their relative magnitudes with respect to a common reference makes sense. This is why we normalize the local displacement $\hat {\eta }({\boldsymbol {x}}_0)$ by $\hat {T}_x$ for odd modes and $\hat {T}_z$ for even modes, in order to make the translational displacement of the bubble centroid unity in all cases.

Thanks to this procedure, the evolution of the bubble shape and the relative magnitude of the interface displacements induced by the rigid-body rotation and volume-preserving deformations may be tracked along each branch of the bifurcation diagram. To save space, we do not discuss shape evolutions here (examples are provided by Bonnefis Reference Bonnefis2019), and focus on the relative variations of the different contributions and their physical origin.

6.2. Variations of the inclination/deformation-induced contributions and path styles

Figure 16(a) considers the case of the most amplified rigid mode for a bubble rising in water at $20\,^\circ {\rm C}$, i.e. the low-frequency oscillating mode up to $Bo\approx 0.64$, followed by the upper branch of the stationary mode for higher $Bo$ (figure 13). It shows how the maximum over the bubble contour of the three contributions to $\hat {\eta }$ varies with the Bond number. In the vicinity of the threshold ($Bo=0.463$), the inclination-induced displacement is seen to be essentially in quadrature with the lateral drift, which is the expected behaviour in a planar zigzagging motion if the bubble axis remains aligned with the path (i.e. its inclination is maximum when the bubble reaches the centreline of its zigzagging path). The relative magnitude of the maximum normal displacement resulting from this inclination is $7\times 10^{-2}$ and experiences little variation up to $Bo\approx 0.55$. At the threshold, the maximum of the volume-preserving deformation is almost in phase with the lateral drift. This is because the rise speed reaches a maximum at each extremity of the zigzag (Mougin & Magnaudet Reference Mougin and Magnaudet2006), resulting in some flattening of the bubble. Conversely, the rise speed is minimum when the bubble crosses the centreline of the zigzag, which slightly reduces its oblateness, yielding the small non-zero imaginary component of the deformation (dashed green line in figure 16a). The relative magnitude of the deformation does not exceed $3\times 10^{-3}$ in this $Bo$ range. This provides the confirmation that time-dependent deformations are extremely weak at the onset of path instability for bubbles rising in low-$Mo$ liquids. These findings also validate the ‘rigid’ qualification employed for the primary low-frequency mode throughout § 5. Conversely, the relative magnitude of deformations in modes $(2,0)$ and $(2,1)$ (not shown) is almost unity at their respective thresholds. The real (respectively imaginary) part of the displacements associated with the bubble inclination (respectively deformation) is typically one order of magnitude smaller than its imaginary (respectively real) counterpart. These weak components exhibit abrupt variations close to the threshold (and for one of them around $Bo=0.6$), presumably because such low-magnitude quantities are only determined with a limited accuracy by the least-squares technique.

Figure 16. Decomposition of the normal displacement of the bubble surface into a horizontal uniform translation, an inclination-induced displacement resulting from a rigid-body rotation, and a volume-preserving deformation. (a) Variation with $Bo$ for the most amplified rigid mode $|m|=1$ in the case of a bubble rising in water at $20\,^\circ {\rm C}$; the yellow vertical line indicates the threshold of path instability; (b) variation with the fluid properties at the path instability threshold. Each line figures the maximum of the corresponding component over the bubble contour; for each component, the solid and dashed lines refer to the real and imaginary parts, respectively.

These two weak components increase sharply with $Bo$ and become dominant for $Bo\gtrsim 0.55$, marking a change in the style of the bubble ascent. In particular, the real part of the displacement associated with the bubble inclination exceeds the imaginary part for $Bo\geqslant 0.59$. This behaviour indicates that the bubble now exhibits a non-zero inclination at the extremities of its zigzagging path, so that its minor axis is no longer aligned with its velocity. In between two successive extremities of the zigzag, the bubble now glides along its path rather than facing it. This makes the part of its surface pointing in the direction of the lateral drift slightly more pointed, while the opposite part becomes slightly more rounded. This is the reason why the imaginary part of the deformation is significant in this second style of zigzagging motion. Conversely, owing to the non-zero bubble inclination at the extremities of the zigzag, the rise speed is slightly lower than in the previous zigzag style. This mitigates the flattening of the bubble, thus reducing the real part of the volume-preserving deformation, as observed on the solid green line of figure 16(a) in the range $0.5\lesssim Bo\lesssim 0.6$.

At $Bo=0.64$, the imaginary parts of the rotation- and deformation-induced contributions fall to zero, a consequence of the change in nature of the primary rigid mode at the exceptional point revealed by figure 13. Beyond this point, the rigid mode is stationary, yielding a non-oscillatory inclined path. When the Bond number increases in the range 0.6–0.8, the bubble gradually flattens and the size of the standing eddy at its back grows significantly, as panels (h,i) of figure 5 show. This makes the strength of the trailing vortices that set in at the path instability threshold increase with $Bo$, which translates directly into an increase of the lift force acting on the bubble, hence an increase of its inclination with the Bond number. This scenario is confirmed by the solid red line in figure 16(a), which shows that the maximum displacement associated with the bubble inclination increases by a factor of four from $Bo=0.64$ to $Bo=0.8$. The larger the bubble inclination, the lower its rise speed, which in turn makes its oblateness reduce compared with that in the base state. This is the origin of the sharp increase of the real part of the volume-preserving deformation (solid green line in figure 16a), which gains two orders of magnitude between $Bo=0.6$ and $Bo=0.8$. Nevertheless, deformations remain small compared with the lateral drift, reaching only a relative magnitude of $3\times 10^{-2}$ at $Bo=0.8$.

Figure 16(b) shows how the relative magnitude of the three contributions to $\hat {\eta }$ vary with the Morton number at the path instability threshold, i.e. along the neutral curve of figure 7. Since the corresponding unstable mode is always oscillatory, one expects the relative magnitude of the different contributions to behave qualitatively in the same way as in the case of water in the vicinity of the primary threshold. Indeed, in line with the behaviours noticed above, deformation-induced displacements are essentially in phase with the bubble lateral drift throughout the eight decades of $Mo$ explored here, while those associated with the bubble inclination are almost in quadrature. Again, we attribute the non-smooth variations of the weak real (respectively imaginary) component of the inclination-induced (respectively deformation-induced) contribution observed for $Mo\lesssim 10^{-10}$ to the limitations of the minimization technique. In contrast, the dominant component of these two contributions is seen to increase smoothly with the Morton number. The relative deformation-induced displacements, say $\hat {\eta }^{def}_{max}|_0$, follow the approximate power law $\hat {\eta }^{def}_{max}|_0\sim Mo^{0.23}$, increasing by two orders of magnitude from Galinstan ($\hat {\eta }^{def}_{max}|_0\approx 10^{-3}$) to DMS-T11 ($\hat {\eta }^{def}_{max}|_0\approx 7\times 10^{-2}$). The relative inclination-induced displacements only grow by one order of magnitude, from $4\times 10^{-2}$ to $5\times 10^{-1}$, in between the same two limits. Nevertheless, the fact that these displacements reach $50\,\%$ of the lateral drift for DMS-T11 (to be compared with $7\,\%$ in the case of water) indicates that bubbles wobble significantly when their path becomes unstable in liquids having a Morton number four to five orders of magnitude higher than that of water.

7.1. Main findings

Thanks to an innovative numerical approach allowing the linear stability of the viscous flow past freely rising and deforming bubbles to be assessed accurately, even when the rise Reynolds number is large, we explored the conditions under which the path of such bubbles deviates from a straight vertical line and the physical nature of the underlying instability modes over a wide range of properties of the carrying liquids. Present results extend over eight orders of magnitude of the Morton number. Over this range, the critical Bond number increases by a factor of $35$ from Galinstan to silicone oil DMS-T11 (yielding an increase by a factor of $2.6$ of the critical bubble size), whereas the critical Reynolds number decreases by a factor of $30$, from ${\approx}\,2100$ to $70$, in between the same two liquids. We showed that path instability always arises through the destabilization of a non-axisymmetric oscillatory mode associated with weak time-dependent changes of the bubble shape. Increasing the bubble size in a given fluid, the nature of this primary ‘rigid’ mode stays unchanged if the Morton number is large enough, say $Mo\gtrsim 10^{-7}$. Conversely, this mode gives birth to two stationary modes if the Morton number is low enough ($Mo\lesssim 10^{-9}$). In the intermediate range $10^{-9}\lesssim Mo\lesssim 10^{-7}$, this mode remains oscillatory even far from the primary threshold but its growth rate varies non-monotonically with the bubble size, and it may even re-stabilize within a narrow size range. A second rigid mode, which turns to be stationary for intermediate and ‘high’ Morton numbers but oscillatory for low Morton numbers also arises further away from the primary threshold. ‘Shape’ modes distinct from the previous rigid modes in that they exhibit much larger oscillations of the bubble surface also become unstable as the bubble size is increased.

Comparison of present predictions with reference experimental and numerical data revealed a very good agreement, both on the threshold (most often measured through the Bond or Galilei number, or alternatively the rise Reynolds number) and the reduced frequency of path oscillations (measured through a Strouhal number). This agreement fully supports the view that the path instability of real gas bubbles results from a linear mechanism associated with a Hopf bifurcation. In contrast, oblique paths which are the hallmark of stationary modes have not been reported in experiments, nor in fully resolved simulations. We see their absence as an indication that the changes introduced in the base flow by the primary non-axisymmetric oscillatory mode are significant enough to alter the next steps of the actual bifurcation sequence compared with the predictions of the linear approach.

To better quantify time-dependent bubble deformations and their influence on the path instability mechanisms and characteristics, we employed two additional and complementary tools. First, a least-squares minimization technique allowed us to split the displacements of the bubble–fluid interface into contributions associated with a rigid-body motion and volume-preserving deformations. By doing so, we could among other things establish how the amplitude of the deformations varies with the Morton number at the onset of path instability. It turned out that the ratio of the maximum deformation to the maximum lateral drift of the bubble centroid is always much smaller than unity. In particular, this ratio is approximately $10^{-3}$ for Galinstan and $3\times 10^{-3}$ for water. It exceeds $1\,\%$ only for liquids with $Mo\gtrsim 1.5\times 10^{-8}$, reaching $7\,\%$ for DMS-T11 ($Mo=9.9\times 10^{-6}$). Second, we took advantage of the base states computed with the present numerical technique to update the predictions obtained in the framework of the FSA by Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martinez-Bazán2016a). With this procedure, any difference in the nature of the primary bifurcation or/and the corresponding threshold may be ascribed to effects of time-dependent bubble deformations. We found that, despite their weak magnitude, these deformations lower significantly the threshold whatever the Morton number. Hence, a general conclusion of the present study is that deformations always promote path instability in the linear framework. Differences between the predicted critical size of deformable and non-deformable bubbles increase gradually from $8\,\%$ at $Mo\approx 10^{-13}$ to $23\,\%$ at $Mo\approx 6\times 10^{-7}$, before decreasing sharply to nearly $15\,\%$ for larger $Mo$. For $Mo\lesssim 7\times 10^{-11}$, FSA correctly predicts that the most unstable mode is a ‘low-frequency’ one with $St={O}(0.01)$ at the threshold. Similarly, it rightly predicts that path instability arises through a ‘high-frequency’ mode with $St={O}(0.1)$ at the threshold for $Mo\gtrsim 3\times 10^{-7}$. The predicted reduced frequencies lie $20$ to $30\,\%$ below those obtained with deformable bubbles, but their variations with the Morton number are qualitatively similar. In contrast, FSA dramatically fails to detect that the most unstable mode is still oscillatory in the intermediate-Morton-number range $7\times 10^{-11}\lesssim Mo\lesssim 3\times 10^{-7}$.