domain: phys spi

text:
Paul Bonnefis

Javier Sierra-Ausin

David Fabre

Jacques Magnaudet

Path instability of deformable bubbles rising in Newtonian liquids: A linear study

The first stages of the path instability phenomenon affecting the buoyancy-driven motion of gas bubbles rising in weakly or moderately viscous liquids are examined using a recently developed numerical approach designed to assess the global linear stability of incompressible flows involving freely-evolving interfaces. Predictions for the critical bubble size and frequency of the most unstable mode are found to agree well with reference data obtained in ultrapure water and in several silicone oils. By varying the bubble size, stability diagrams are built for several specific fluids, revealing three distinct regimes with different bifurcation sequences. The spatial structure of the unstable modes is analysed, together with the variations of the bubble shape, position and orientation. For this purpose, displacements of the bubble surface are split into rigidbody components and volume-preserving deformations, allowing us to determine how the relative magnitude of the latter varies with the fluid properties and bubble size. Predictions obtained with freely-deformable bubbles are compared with those found by maintaining the bubble shape determined in the base state frozen during the stability analysis. This comparison reveals that deformations leave the phenomenology of the first bifurcations unchanged in low-viscosity fluids, especially water. Hence, in such fluids, bubbles behave essentially as freely-moving rigid bodies submitted to constant-force and zero-torque constraints, at the surface of which the fluid obeys a shear-free condition. In contrast, deformations change the nature of the primary bifurcation in oils slightly more viscous than water, whereas, somewhat surprisingly, they leave the near-threshold phenomenology unchanged in more viscous oils.

Introduction

A great deal of effort has been invested since the middle of the last century to elucidate the origin and predict the characteristics of the intriguing path instability phenomenon affecting the buoyancy-driven motion of millimeter-sized gas bubbles rising in weakly or moderately viscous liquids. Attempts were mostly experimental during the first decades, before numerical simulation techniques became mature enough to shed some new light on this puzzling phenomenon. Among them, linear stability analyses have kept pace of advances in the computational techniques required to determine efficiently the discretized form of the stationary solutions of the Navier-Stokes equations subjected to appropriate boundary conditions, and to solve subsequently the large-size generalized eigenvalues problems giving access to the possibly linearly unstable eigenmodes.

Early studies considered bubbles with a frozen shape held fixed in a uniform stream. Assuming an oblate spheroidal shape with an arbitrary aspect ratio (the ratio of the major and minor axes lengths, hereinafter denoted as χ), with the short axis of the spheroid aligned with the incoming flow, conditions under which the axisymmetric wake becomes unstable, the nature of the corresponding first bifurcations and the structure of the associated unstable modes could be explored [START_REF] Yang | Linear stability of the flow past a spheroidal bubble[END_REF][START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF]. It was eventually concluded that wake instability occurs within a finite range of Reynolds number when χ 2.21, in agreement with predictions provided by the numerical solution of the full Navier-Stokes equations in the same configuration [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF]. The first bifurcation is stationary, i.e. the corresponding eigenvalue is real. Beyond this bifurcation, the stationary wake exhibits a pair of counterrotating trailing vortices which are responsible for a nonzero lift force. For larger aspect ratios, increasing the Reynolds number beyond the primary threshold while keeping χ fixed or vice versa, a second bifurcation of Hopf type takes place. The wake keeps the previous planar symmetry unchanged beyond the corresponding threshold but vortices are periodically shed downstream, resulting in periodic oscillations of the lift force about a nonzero mean. This transition scenario, with the axisymmetric wake becoming unstable through a stationary bifurcation and the resulting non-axisymmetric stationary wake undergoing a Hopf bifurcation, is similar to that obeyed by the wake of rigid spheres and discs [START_REF] Natarajan | The instability of the steady flow past spheres and disks[END_REF].

Compared with the actual problem of buoyancy-driven rising bubbles, the configuration contemplated in the above studies is purely academic, as the bubble is held fixed despite the existence of a nonzero lift force. An important step towards physically realistic conditions was achieved during the last decade, starting with global linear stability studies of the flow past freely moving buoyancy/gravity-driven rigid bodies with various shapes. In this situation, the body motion is subjected to two additional constraints, since the hydrodynamic force must balance the net body weight at all times (implying that the horizontal force remains null whatever the body orientation and velocity), and the hydrodynamic torque must also remain zero, provided that the body density is uniform. Therefore, disturbances in the body velocity or orientation modify the velocity and pressure distributions in the surrounding fluid, which in turn induces changes in the force and torque acting on the body that must respect the above contraints. Two-dimensional plates and rods were considered by [START_REF] Fabre | A quasi-static approach to the stability of the path of heavy bodies falling within a viscous fluid[END_REF] and [START_REF] Assemat | The onset of unsteadiness of two-dimensional bodies falling or rising freely in a viscous fluid: a linear study[END_REF]. The former study focused on 'heavy' bodies with a density much larger than that of the fluid. An asymptotic approach could then be employed and revealed the existence of four 'aerodynamic' modes distinct from the 'fluid' (or 'wake') modes associated with wake instability (qualitatively similar modes are encountered in flight dynamics, and those evidenced in the above studies were coined after them). Two of these aerodynamic modes have real negative eigenvalues, i.e. they are always stable. In contrast, the other two are associated with a pair of complex conjugate eigenvalues, hence with a Hopf bifurcation. The corresponding oscillations are slow compared with those typical of vortex shedding in two-dimensional wakes, and their frequency varies as the inverse square root of the body-to-fluid density ratio, m * , thus becoming vanishingly small for very heavy bodies. Existence of these modes was confirmed for arbitrary m * in the global linear stability analysis (hereinafter abbreviated as GLSA) performed numerically by [START_REF] Assemat | The onset of unsteadiness of two-dimensional bodies falling or rising freely in a viscous fluid: a linear study[END_REF], and the corresponding thresholds were compared with those of the fluid modes for thin plates and square rods. Increasing the Reynolds number, it was observed that, for infinitely thin plates, the primary wake mode always becomes unstable first, no matter what the value of m * is. The same conclusion was found to hold for square rods as long as the density ratio exceeds 1.22. In contrast, for lower m * , the aerodynamic oscillatory mode is destabilized first and exhibits a frequency (corresponding to a fluttering of the falling/rising rod about the vertical) smaller than that of the wake mode by about a factor of two. Similar conclusions were previously obtained through fully-resolved numerical simulations by [START_REF] Alben | An implicit method for coupled flow-body dynamics[END_REF] with light two-dimensional ellipses, and the critical Reynolds number at which the fluttering motion sets in was found to be approximately half that corresponding to the onset of wake instability. The two-dimensional GLSA approach of [START_REF] Assemat | The onset of unsteadiness of two-dimensional bodies falling or rising freely in a viscous fluid: a linear study[END_REF] was extended to axisymmetric bodies by Tchoufag et al. (2014a) who considered the stability of rising/falling discs and short cylinders. Whatever the body aspect ratio, they found the aerodynamic low-frequency mode to be always the most unstable beyond a critical χ-dependent m * (note the difference with two-dimensional square rods for which the same happens below a critical m * ). For lower m * , the most unstable mode is either the one associated with the wake (here characterized by weak deviations of the path in the presence of sustained fluid oscillations at the back of the body), or a stationary mode in which the body follows an inclined path and its wake is made of a pair of steady counter-rotating streamwise vortices. Based on these findings, Tchoufag et al. (2014a) concluded that the first deviation of the path of such freely moving axisymmetric bodies cannot in general be predicted by considering the instability of the sole wake. Rather, the instability of their path results under most conditions from the intrinsic coupling between the body and fluid implied by the constant-force and zerotorque conditions. Tchoufag et al. (2014b) applied the GLSA approach to freely rising spheroidal bubbles with a prescribed oblateness. They observed that the vertical path becomes unstable within a finite range of Reynolds number when the bubble oblateness exceeds 2.15, which is slightly lower than the wake instability threshold determined for a fixed bubble [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF][START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF]. The unstable Reynolds number range widens as χ increases. For 2.15 < χ < 2.25, increasing the Reynolds number, i.e. the bubble size, while maintaining the oblateness fixed, the path first becomes unstable within a narrow Reynolds number range through a low-frequency mode similar to the aerodynamic mode encountered with heavy short cylinders and discs. Increasing the Reynolds number further, a stationary mode yielding an inclined path and growing much faster than the previous low-frequency mode becomes dominant. This situation prevails within an intermediate Reynolds number range, beyond which the low-frequency mode becomes dominant again. Last, the system restabilises beyond a fourth critical Reynolds number. For instance, the path of a bubble with χ = 2.22 becomes oscillatory through a low-frequency mode in the two separate ranges 310 Re 340 and 450 Re 980, while it takes a nonzero stationary inclination with respect to the vertical in the intermediate range 340 Re 450, the Reynolds number Re being based on the bubble rise speed and equivalent diameter. Obviously, the description provided by this simplified model suffers from the fact that the bubble shape is not allowed to vary with the Reynolds number, in contrast to what happens under real conditions. Nevertheless, as will become apparent later, the first half of the above scenario, i.e. the succession of a low-frequency mode becoming unstable first and a stationary mode taking over it, is relevant to the case of deformable bubbles rising in low-viscosity high-surface-tension fluids, especially water.

A somewhat more sophisticated track was followed by Cano-Lozano et al. (2016a). Rather than prescribing arbitrarily an oblate spheroidal shape, they computed the actual bubble shape and rise speed corresponding to the stationary axisymmetric base state, using the open software Gerris [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF][START_REF] Popinet | The Gerris flow solver[END_REF]. Then, they introduced this shape and speed into the GLSA solver of Tchoufag et al. (2014b), keeping the shape frozen despite the possible development of flow and path instabilities. With this approach, they could determine for a wide variety of Newtonian fluids (characterized by their density, viscosity and surface tension), the critical Reynolds number, i.e. the critical bubble size, beyond which the vertical path of a bubble with a nearly realistic shape becomes unstable. They could also compare these predictions with those corresponding to 'pure' wake instability, obtained by holding the same bubble fixed in a uniform stream. For lowviscosity high-surface-tension fluids, they found that, in agreement with the conclusions of Tchoufag et al. (2014b), the system first becomes unstable through a Hopf bifurcation leading to low-frequency path oscillations. At the corresponding Reynolds number, the wake of the same bubble held fixed is still stable, which proves that the mechanism responsible for path instability is intimately linked to the constant-force and torque-free conditions that constrain the dynamics of the coupled fluid-body system. For liquids with a viscosity 5 to 10 times larger than water and a surface tension 3 -4 times smaller, they observed that path instability first arises through an oscillatory mode with a frequency 4 -5 times higher than that found in low-viscosity high-surface tension liquids. However, the picture looked dramatically different in the intermediate range corresponding to 'moderate' liquid viscosities and surface tensions. In particular, the instability thresholds of freely-moving and fixed bubbles were found to be identical in that range. Hence, the first unstable mode is stationary, similar to that observed for fixed spheroidal bubbles [START_REF] Yang | Linear stability of the flow past a spheroidal bubble[END_REF][START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF][START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF]. When compared with reference experiments carried out in ultrapure water [START_REF] Duineveld | The rise velocity and shape of bubbles in pure water at high Reynolds number[END_REF][START_REF] De Vries | Path and wake of a rising bubble[END_REF] or silicone oils [START_REF] Zenit | Path instability of rising spheroidal air bubbles[END_REF][START_REF] Sato | Deformation, wake and collision of rising bubbles[END_REF]), these predictions experience mixed successes. In low-viscosity high-surface-tension fluids, they properly capture the low-frequency oscillatory nature of the incipient path instability, and slightly over-predict the size of the smallest bubble whose path becomes unstable. The same conclusion holds for the 'high'-frequency path oscillations observed in oils 5 -10 times more viscous than water. In contrast, for oils with an intermediate viscosity, experimental observations provide evidence that the first instability of the system keeps a similar oscillatory nature, at odds with the above predictions.

This blatant disagreement points to the role of transient bubble deformations not accounted for in the above studies. It provided one of the main motivations that decided us to undertake the development of a more general linear stability approach capable of dealing with free gas-liquid boundaries, thanks to which the fate of time-dependent deformations of the bubble-fluid interface may be predicted together with that of flow and path disturbances. The corresponding method was successfully developed and validated by [START_REF] Bonnefis | Etude des instabilités de sillage, de forme et de trajectoire de bulles par une approche de stabilité linéaire globale[END_REF]; other groups have pursued similar objectives in parallel [START_REF] Zhou | Marginal stability curve of a deformable bubble[END_REF][START_REF] Herrada | Path instability of an air bubble rising in water[END_REF]. These efforts may be considered as the third and final step on the long route of the linear stability analysis of the flow past rising bubbles. To summarize, in the first step the bubble shape was frozen, its centroid was assumed to remain fixed and only the stability of the wake was assessed [START_REF] Yang | Linear stability of the flow past a spheroidal bubble[END_REF][START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF]). In the next step, the bubble shape was still frozen but the bubble centroid was allowed to move freely under the effect of buoyancy, being only constrained by the constant-force and zero-torque conditions (Tchoufag et al. 2014b;Cano-Lozano et al. 2016a). Most of the investigations performed in these first two steps assumed the bubble to keep an exact oblate spheroidal shape, except the last of them in which this shape was obtained as part of the stationary solution of the full Navier-Stokes equations. The time-dependent adjustment of the bubble shape to its possibly non-straight unsteady motion was the last ingredient missing to reach a fully realistic description of the bubbleflow interaction. This is what the third step achieved by [START_REF] Bonnefis | Etude des instabilités de sillage, de forme et de trajectoire de bulles par une approche de stabilité linéaire globale[END_REF] enables. The purpose of this paper is to present and discuss the predictions provided by this approach across a wide range of fluids, from low-viscosity high-surface-tension liquid metals to oils with a viscosity typically ten times larger and a surface tension four times less than those of water, and to compare them with reference data when available. Predictions obtained in the specific case of water were already summarized by [START_REF] Bonnefis | When, how, and why the path of an air bubble rising in pure water becomes unstable[END_REF] and are discussed in more detail here. The numerical approach, already described by [START_REF] Sierra-Ausin | Dynamics of a gas bubble in a straining flow: Deformation, oscillations, self-propulsion[END_REF], is summarized in § 2. The characteristics of the base state are presented in § 3, before the neutral curves resulting from the linear stability analysis are discussed in § 4, and the nature and sequencing of the first unstable modes are examined in § 5. Then, the relative magnitude of deformations of the bubble-fluid interface with respect to the lateral drift of the bubble centroid are discussed in § 6. Section 7 summarizes the main findings of the study and discusses the respective roles of wake, bubble deformation and fluid-bubble hydrodynamic couplings in the different regimes encountered by varying the physical properties of the carrying liquid.

Problem formulation and numerical approach

We assume that a single gas bubble with volume V and equivalent diameter D = ( 6 π V) 1/3 rises through a large expanse of a Newtonian liquid with density ρ, viscosity µ and surface tension γ under the effect of gravity g. The liquid is at rest at infinity and the disturbance flow induced by the bubble ascent is assumed incompressible. The bubble is initially axisymmetric and is released with its symmetry axis aligned with gravity. Considering that the density and viscosity of the gas enclosed in the bubble are negligibly small compared with those of the surrounding liquid, the problem depends on two dimensionless parameters, namely the Bond number Bo = ρgD 2 /γ and the Galilei number Ga = ρ(gD 3 ) 1/2 /µ. These characteristic numbers compare the gravitational force ρgD 3 driving the bubble ascent with the capillary force γD and the viscous force µ(gD 3 ) 1/2 , respectively. Whatever the bubble size, a given liquid placed in a given gravitational environment is entirely characterized by the value of the so-called Morton number M o = Bo 3 Ga -4 = gµ 4 /(ργ 3 ). By defining the visco-gravitational diffusion length scale l µ = [µ 2 /(ρ 2 g)] 1/3 (the length over which viscosity diffuses a disturbance within a [µ/(ρg 2 )] 1/3 -long period of time) and the capillary length scale l γ = [γ/(ρg)] 1/2 , it is seen that l µ /l γ = M o 1/6 . Hence, the Morton number characterizes the ratio of these two length scales. The bubble diameter, D, the gravitational time, (D/g) 1/2 , and their ratio, (gD) 1/2 , will be used throughout the paper to make lengths, frequencies and velocities dimensionless. Once the bubble rise speed u b is known, the non-dimensional rise speed, U = u b /(gD) 1/2 allows the bubble Reynolds number, Re = Ga U , and Weber number, W e = Bo U 2 , to be evaluated. The stability of the coupled fluid-bubble motion is examined with the help of the Linearized Arbitrary Lagrangian-Eulerian approach, hereinafter abbreviated as L-ALE, developed by [START_REF] Bonnefis | Etude des instabilités de sillage, de forme et de trajectoire de bulles par une approche de stabilité linéaire globale[END_REF] and described by [START_REF] Sierra-Ausin | Dynamics of a gas bubble in a straining flow: Deformation, oscillations, self-propulsion[END_REF]. Here, we only summarize the main steps of this approach, referring readers to the above references for more details. The key idea is to solve the governing equations of the problem, which hold on a time-deforming domain Ω(t) since the bubble-fluid interface Γ b (t) experiences time-dependent deformations, on a prescribed reference domain Ω 0 bounded internally by a fixed surface Γ 0 representing the reference position of the interface (figure 1). In addition to Γ 0 , Ω 0 is bounded by a fixed external surface Γ ∞ on which suitable far-field conditions are imposed, and Γ 0 and Γ ∞ are connected by the vertical axis, Γ z , about which the base state exhibits an axial symmetry. The L-ALE approach is grounded on the explicit building of the bijective mapping connecting the instantaneous position x of a given geometrical point in the deforming domain to its position x 0 in Ω 0 . In other words, this approach relies on the determination of the infinitesimal displacement field ξ(x 0 ) such that x = x 0 + ξ(x 0 ) everywhere in Ω 0 ∪ Γ 0 ∪ Γ z ∪ Γ ∞ . On Γ 0 , the unit normal of which is n 0 , the normal component of the displacement, ξ • n 0 , must equal the displacement of the bubble-fluid interface, η. In contrast, ξ may be arbitrarily chosen anywhere else. To ensure a smooth distribution of the displacement field, we assume that ξ obeys the Cauchy-Navier equation of elastostatics with unit Lamé coefficients in Ω 0 , together with a suitable symmetry condition on Γ z and a Dirichlet condition ξ = 0 on Γ ∞ .

Knowing ξ(x 0 ) and its gradients everywhere, the time and space derivatives involved in the governing equations written in Ω(t) may be expressed in Ω 0 and linearized consistently with respect to ξ, together with the variations of the unit normal and local mean curvature of Γ b (t). Governing equations in Ω 0 were detailed in [START_REF] Sierra-Ausin | Dynamics of a gas bubble in a straining flow: Deformation, oscillations, self-propulsion[END_REF]. They consist of the continuity and Navier-Stokes equations, while the kinematic (no-penetration) and dynamic boundary conditions must be enforced on Γ 0 . Owing to the negligible viscosity of the gas enclosed in the bubble, the tangential component of the dynamic boundary condition reduces to a shear-free condition that must be satisfied on Γ 0 by the gradients of the liquid velocity field, u. Since the gas density is also negligible, the momentum balance implies that the pressure is uniform (but possibly time-dependent) within the bubble. Consequently, the normal component of the dynamic boundary condition on Γ 0 expresses the fact that the difference between the uniform pressure p b inside the bubble and the local pressure p in the liquid (from which the hydrostatic component has been subtracted) balances the difference between the local capillary pressure and the normal viscous stress in the liquid. Last, two integral conditions have to be added to ensure that the volume of the bubble does not change over time, and that the vertical position of its geometrical centroid does not vary once the problem has been expressed in a reference frame translating with the bubble rise speed, u b .

A key feature of the L-ALE approach is that the entire set of equations governing the evolution of the state vector q = [u, p, η, ξ, u b , p b ] T ( T denoting the transpose) is solved simultaneously (such an approach is frequently referred to as 'monolithic').

Starting from an arbitrary initial solution q 0 = [u 0 , p 0 , 0, 0, u b0 , p b0 ] T (here assumed to be axisymmetric) and dropping time derivatives in the Navier-Stokes equations and in the no-penetration condition, a Newton algorithm is used to obtain the stationary base state through a series of global iterations in which the state vector and the reference domain are iteratively updated. For this purpose, we express the problem in weak form and make use of the finite element software FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF] to build and invert the various matrices involved. The volume fields (u, p, ξ) are discretized on the finite element basis suitable for each of them, while the normal displacement of the interface, η, is discretized on the local Fourier basis. Provided that q 0 is chosen 'not too far' from the stationary solution (i.e. within its basin of attraction), the above algorithm converges quadratically in a few iterations. Moreover, compared with time-marching algorithms, it has the decisive advantage that unstable steady solutions may also be captured, allowing the features of the corresponding branches (if any) of the bifurcation diagram to be studied in detail; see, e.g., [START_REF] Sierra-Ausin | Dynamics of a gas bubble in a straining flow: Deformation, oscillations, self-propulsion[END_REF].

Once the state vector and the shape of the bubble-fluid interface defining the axisymmetric base state have been obtained, the linear stability of this solution is assessed in the Galilean reference frame translating with the corresponding rise speed. For this purpose, we consider a cylindrical coordinate system (r, θ, z) whose axis r = 0 and cross-sectional plane z = 0 correspond to the symmetry axis and vertical position of the bubble geometrical centroid in the base state, respectively. Then we compute the complex eigenvalues λ = λ r + iλ i associated with disturbances of the form q (r, θ, z, t) = [(û, p, η, ξ)e imθ , pb ] T e λt , with m denoting the azimuthal wave number and θ the azimuthal angle, the hatted amplitudes (û, p, η, ξ) depending on both r and z (owing to the choice of the reference frame, ûb ≡ 0, so that the instantaneous position of the bubble centroid no longer necessarily coincides with that of the reference frame). In what follows, we shall be primarily interested in the family of non-axisymmetric modes |m| = 1 which are those associated with the first deviations of the bubble path from the vertical, and in modes m = 0 and |m| = 2 associated with the oscillations of the bubble shape. The eigenvector q is obtained via the Newton algorithm already employed to determine the base state. The set of equations to be solved is similar, except that the momentum equation and the no-penetration condition now involve terms proportional to λ since the sought solution is time-dependent; the reference state vector and domain to be considered are those corresponding to the base state. The complex eigenvalues are finally obtained using a Krylov-Shur projection method available in the SLEPc library.

Obviously, the accuracy with which the base state and the eigenvalues are determined depends critically on the finite element triangulation of Ω 0 , especially along the bubblefluid interface Γ 0 and within the boundary layer and wake regions. The position of the fictitious outer boundary Γ ∞ with respect to the bubble is also important, the computational domain having to be large enough to avoid spurious confinement effects. It must also be kept in mind that, in low-Morton-number fluids such as water (M o ≈ 2.54×10 -11 under standard conditions) and even more in liquid metals, such as Galinstan (M o ≈ 1.4 × 10 -13 ), the onset of path instability takes place in regimes where the bubble Reynolds number Re = ρu b D/µ is of the order of 10 3 (≈ 670 in water, ≈ 2100 in Galinstan, see figure 7 below). To ensure that the flow is still fully resolved in the thin boundary layers encountered in such regimes, we design an initial grid in which 18 nodes are distributed across a 5 × Re -1/2 -thick layer surrounding the bubble. Along the bubble surface, the number of nodes N D over a D-long arc length is approximately 7.2 × Re 1/2 , with at least N D = 64 when Re becomes less than 80. Since u b , hence Re, is not known a priori, once an approximate steady state has been computed on the initial grid, a new grid is generated following an adaptive mesh refinement procedure. This new grid 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 mN m -1 , depending on the exact alloy composition.

Liquid ρ (kg m -3 ) µ (mPa s) γ (mN m -1 ) M o = gµ 4 /(
is designed based on the previous steady-state solution, with some constraints ensuring that the boundary layer remains sufficiently resolved. This procedure is repeated a few times until the solution becomes grid-independent. Extensive tests were carried out to evaluate the sensitivity of the base state (appreciated for instance through the variations of Re, χ, λ r and λ i ) to factors such as the domain size, distance from the bubble to the downstream boundary and density of nodes along the bubble surface. For instance, in the case of a bubble with Bo = 0.6 rising in water, it was found that relative variations of these indicators for the non-axisymmetric modes m = ±1 are all less than 0.1% provided that N D 70 and the distance from the bubble centroid to the downstream end of the domain is at least 30D, the upstream end and the lateral surface being both located 15D away from this centroid (Bonnefis 2019).

Base state

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 ∝ 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 (M o = 6.2 × 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 5 in DMS-T05 or Bo 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.

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 [START_REF] Zenit | Path instability of rising spheroidal air bubbles[END_REF]. In each series, the location of the incipient path instability detected in the experiments is specified in the figure. The threshold Reynolds (data for all four series taken from [START_REF] Zenit | Path instability of rising spheroidal air bubbles[END_REF]). Red bullets mark the onset of path instability detected experimentally in each series.

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. 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 on a master curve following the power law Re ∝ 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 χ is close to unity. For low enough W e, departures from sphericity are known to increase linearly with the Weber number according to the law χ = 1+ 9

64 W e [START_REF] Moore | The velocity of rise of distorted gas bubbles in a liquid of small viscosity[END_REF]. Since W e = Re 2 M o 1/3 Ga -2/3 , the lower M o the wider the Re range pertaining to this first regime, which extends up to Re ≈ 200 in water. Bubbles become flatter as W e, 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 W e which itself depends on M o, the Reynolds number is no longer independent of the Morton number in this second regime. To further examine the two behaviours, on has to refer to the force balance on the bubble which, in non-dimensional form, reads C D (Re)Re 2 = 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 [START_REF] Batchelor | An Introduction to Fluid Dynamics[END_REF]. Hence, for nearly-spherical bubbles, the force balance implies k(Re)Re = 1 12 Ga 2 , with k(Re 1) = 1, k(Re 1) = 3. The slow increase of k(Re) in between these two limits, which follows approximately the power law k(Re) ∼ 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 W e = O(1) (W e ≈ 0.62 at Re = 200 in water). Assuming that their frontal area, which is proportional to χ 2/3 , grows approximately linearly with the Weber number in that range yields C D (Re) ∼ k(Re)Re -1 W e ∼ Re 6/5 M o 1/3 Ga -2/3 . The force balance then implies Re ∼ Ga 5/6 M o -5/48 , a prediction supported by the numerical data for the various fluids in the Ga range where path instability occurs.

(a) (b) (c) (d) (e) (f) (g) (h) (i)

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 [START_REF] Batchelor | An Introduction to Fluid Dynamics[END_REF], both of which are maximum there. This vorticity results directly from the shearfree 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 [START_REF] Moore | The boundary layer on a spherical gas bubble[END_REF][START_REF] Moore | The velocity of rise of distorted gas bubbles in a liquid of small viscosity[END_REF]. The wake structure is found to depend dramatically on the Bond number and, for a given Bo, on the fluid properties (compare subfigures (a) and (f ), both for Bo = 3). In water, all streamlines are open for Bo = 0.4 (subfigure (g)) and they remain so up to Bo ≈ 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 (subfigure (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 W e ≡ ρu 2 b D/γ = Bo U 2 . Indeed, keeping Bo unchanged, the Weber numbers in two different fluids, 1 and 2, obey the relation W e 2 /W e 1 = (Re 2 /Re 1 ) 2 (M o 2 /M o 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 2.4.

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 flow 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 [START_REF] Cano-Lozano | Wake instability of a fixed axisymmetric bubble of realistic shape[END_REF]. It will become clear in the next section that, in low-M o 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 [START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF]. 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 [START_REF] Chomaz | Global instabilities in spatially developing flows: Non-normality and nonlinearity[END_REF]. That path instability may happen in the absence of such a near-wake structure proves that, in the relevant M o-range, the mechanism that governs this instability is not driven by wake instability per se. This confirms the findings of Cano-Lozano et al. (2016a) who established that, for M o 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 destabilised first.

Neutral curves

Critical bubble size

Using the procedure described in § 2, we selected a set of M o 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 = ±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 -10 -5 , this curve readily provides the size of the smallest bubble whose vertical path becomes unstable. For instance, the path of a bubble rising in water at a temperature of 20 • C becomes unstable at Bo = 0.463 and Ga = 250, yielding a critical diameter D ≈ 1.854 mm [START_REF] Bonnefis | Etude des instabilités de sillage, de forme et de trajectoire de bulles par une approche de stabilité linéaire globale[END_REF][START_REF] Bonnefis | When, how, and why the path of an air bubble rising in pure water becomes unstable[END_REF]. Similarly, the critical Bond numbers for Galinstan (M o = 1.4 × 10 -13 ) and DMS-T11 (M o = 9.9 × 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 M o-interval considered here are 1.48 mm and 3.875 mm, respectively.

In figure 7 we also reported 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 [START_REF] Duineveld | The rise velocity and shape of bubbles in pure water at high Reynolds number[END_REF] in ultrapure water at 20 • C yield a critical bubble diameter close to 1.82 mm, which differs by less than 2% from the present prediction. De Vries (2002) also determined the onset of path instability in a tank of ultrapure water with an average temperature of 28 • C (M o = 1.13 × 10 -11 ). However, a 1.1 • C m -1 stabilizing temperature gradient was established to visualize the wake thanks to optical index gradients. 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 in four different silicone oils [START_REF] Zenit | Path instability of rising spheroidal air bubbles[END_REF] and in three other silicone oils in which photochromic dye was added to visualize the wake [START_REF] Sato | Deformation, wake and collision of rising bubbles[END_REF] prove satisfactory in the sense that the neutral curve lies generally within the (Bo, Ga)-interval whose lower (resp. upper) bound corresponds to the largest (resp. smallest) bubble for which a stable (resp. unstable) path could be identified experimentally. The only two exceptions correspond to a probable outlier in Sato's data for M o = 2.35 × 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. (2016b) 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 M o 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 M o 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 M o ≈ 10 -10 . For lower M o, 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-M o 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. (2014b) and Cano-Lozano et al. (2016a). 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. (2016a) (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 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 evidence, deformations always lower the threshold, i.e. they promote path instability. However, for fluids with Morton numbers less than ≈ 5 × 10 -7 , the lower M o 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 M o = 10 -7 , but this overestimate is reduced by a factor of three in the case of Galinstan. FSA actually predicts that the neutral curve involves two distinct modes intersecting at M o ≈ 2.8×10 -7 . In the narrow range 4×10 -7 M o 6.5×10 -7 , the mode that becomes most unstable beyond this intersection (dark grey line) exhibits a S shape implying a destabilisation-restabilisation 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 M o, 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. (2016a) 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

M o 6 × 10 -7 . Out 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 M o 10 -9 , the difference between the two thresholds increases gradually as M o is decreased. For M o 6 × 10 -7 , the threshold of wake instability The yellow curve was obtained by connecting the 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 'frozen-shape' approximation (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 M o = 2.7 × 10 -7 indicates the frequency jump associated with the switching from the stationary mode to the oscillatory mode predicted by FSA. : numerical predictions of Cano-Lozano et al. (2016b); : experimental data of [START_REF] Duineveld | Bouncing and coalescence of two bubbles in water[END_REF] for Bo = 0.54 in ultrapure water; : experimental data of [START_REF] Zenit | Measurements of the streamwise vorticity in the wake of an oscillating bubble[END_REF] for Bo = 3.92 in DMS-T05. 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 path instability predicted by FSA. These findings, combined with the fact that figure 7 reveals that FSA overestimates the threshold of path instability whatever M o, 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.

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 = λ i D/(2πu b ), with λ 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 ≈ 0.02 for M o = 10 -13 to St ≈ 0.12 for M o = 10 -5 . FSA predicts that the first non-vertical path of the bubble is associated with 'low-frequency' oscillations in low-M o liquids (M o 7 × 10 -11 ) and with 'high-frequency' oscillations in liquids with M o 3×10 -7 , in agreement with the conclusions of Cano-Lozano et al. (2016a). Present findings with deformable bubbles are qualitatively consistent with these predictions in both ranges. The reduced frequencies determined in both approaches for low-M o 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-M o fluids the conclusions drawn by [START_REF] Bonnefis | When, how, and why the path of an air bubble rising in pure water becomes unstable[END_REF] in the specific case of water. More surprisingly, the agreement is still reasonable in fluids with M o 3 × 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 M o-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. (2016b) in silicone oils DMS-T02 (M o = 1.6 × 10 -8 ) and DMS-T11 (M o = 9.9 × 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, W e = Bo U 2 , with the liquid properties apparent. Figure 9(a) shows that the critical Reynolds number decreases sharply as M o increases, from ≈ 2070 for M o = 10 -13 to 70 for M o = 10 -5 , via Re = 668 for water at 20 • C. Conversely, the critical Weber number (figure 7(b)) is seen to increase gradually from 3.20 for M o = 10 -13 to 6.36 for M o = 10 -5 , via W e = 3.30 for water. The above values for water are to be compared with those reported by [START_REF] Duineveld | The rise velocity and shape of bubbles in pure water at high Reynolds number[END_REF], namely Re = 658, W e = 3.275, from which they only differ by 1.5% and 0.75%, respectively.

Unstable modes: order of occurrence and physical nature

In this section, we examine the variations of the most unstable eigenvalues of the coupled bubble-fluid system as the size of the bubble (here measured through the Bond number) is varied in the vicinity of the primary threshold. We also analyse the spatial structure of the associated global modes. It will soon become apparent that the order in which the unstable modes follow one another as the Bond number increases, and even in some cases the way the nature of a given mode varies, depend dramatically on the range of fluid properties under consideration, the complexity of the picture increasing as the Morton number is decreased. For this reason, we discuss these features in descending order of M o.

5.1. Silicone oil DMS-T05 (M o = 6.2× 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 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 nonzero 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 λ i /2π ≈ 0.115 at the threshold (throughout this section, eigenvalues are normalized using the gravitational time scale (D/g) 1/2 ). As subfigure (c) suggests, this mode is essentially associated with a lateral drift of the bubble centroid accompanied by some rigid-body rotation. Timedependent 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'.

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 ≈ 5.5 (point B in figure 10(a)). 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 subfigure (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 ≈ 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. Subfigure (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, subfigure (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 ≈ 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 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 ≈ 6.6), and restabilizes itself within the narrow interval 8.65 Bo 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 [START_REF] Meiron | On the stability of gas bubbles rising in an inviscid fluid[END_REF]). 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 e ±2iθ , this (2, 2) mode does not induce any lateral displacement of the bubble. Last, at point E (Bo ≈ 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 subfigure (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 .

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 τ = π 2 λ -1 i . 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 subfigure (b); in a three-dimensional representation, this structure would look like a pair of crescents connected through their horns in the vertical plane θ = ±π/2, with fluid particles rotating in the same direction in every azimuthal plane. As subfigures (cd) make clear, this structure is actually the first of a series that develop further downstream in the wake, with alternating positive and negative directions of rotation.

Silicone oil DMS-T02 (M

o = 1.6 × 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 subfigure (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.

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 ≈ 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 M o-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. We finally consider the case of bubbles rising in pure water maintained at a temperature of 20 • 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 λ i /2π = 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 ≈ 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 [START_REF] Bergholtz | Exceptional topology of non-Hermitian systems[END_REF]). 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 [START_REF] Dangelmayr | The Takens-Bogdanov bifurcation with O(2)symmetry[END_REF]. 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.

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 λ i /2π = 0.162. This mode (identified as rigid osc. II in the legend of figure 13(a)) is associated with an azimuthal wavenumber |m| = 1, exhibits low-frequency oscillations (subfigure (b)), and is not accompanied by significant shape oscillations (subfigure (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 ≈ 0.88 (point E). As subfigures (bc) make clear, this eigenvalue corresponds to the 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. 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 0.4, followed by the mode (2, 0) and then the mode (2, 1). This is in line with the conclusions obtained by [START_REF] Meiron | On the stability of gas bubbles rising in an inviscid fluid[END_REF] 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 that 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-M o-range, say for M o 10 -9 . For instance, the same first five modes become unstable in the same order in the case of DMS-T00 (M o = 1.8 × 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.

Spatial structure of rigid and shape modes

Figure 14 displays the structure of the three successive rigid modes encountered in the sequence analysed above. Subfigures (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 lowfrequency 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. Subfigures (c) and (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 subfigure (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. inclination of the body path resulting from the constant lift force generated by a pair of semi-infinite counter-rotating streamwise vortices.

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 subfigures (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 subfigures (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 subfigures (ab), 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 = λ i /(2πU ), not directly by λ 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 1 and Bo -1 Bo 1; see equation (7-3) and the associated figure in [START_REF] Clift | Bubbles, Drops and Particles[END_REF]).

6. Respective magnitudes of rigid-body motions and time-dependent deformations

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 clarify 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 θ = (0, π) in the cylindrical (r, θ, z) system defined in § 2. At any position x 0 on the undisturbed bubble-fluid interface, the normal displacement of the interface, η(x 0 ), may be decomposed into a rigid motion inducing a lateral (resp. vertical) displacement Tx (resp. Tz ) about the x (resp. z) axis and an inclination ψ resulting from a rotation about the y axis, augmented by a local volume-preserving deformation ζ(x 0 ) in the normal direction. Projecting this decomposition onto the local normal unit vector n 0 , one has η(x 0 ) = Tx (e x • n 0 ) + Tz (e z • n 0 ) + ψ(x 0 × n 0 ) • e y + ζ(x 0 ) , (6.1) with e x , e y , e z the unit vectors in the x, y, z directions, respectively. Odd modes m = ±1 do not induce any translation in the vertical direction, so that Tz = 0 for such modes. Conversely, the lateral displacement Tx and the inclination ψ of the bubble remain null for even modes m = -2, 0, 2. To determine the nonzero displacement and inclination associated with a given mode, we make use of a least-squares fitting technique. That is, we define the local functional F Tx, Tz, ψ

(x 0 ) ≡ η -[ Tx (e x • n 0 ) + Tz (e z • n 0 ) + ψ(x 0 × n 0 ) • e y ]
2 , integrate it over the entire bubble contour and minimize the corresponding residue with respect to the nonzero components of the triplet ( Tx , Tz , ψ). 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 η(x 0 ) by Tx for odd modes and Tz 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 ( 2019)), and focus on the relative variations of the different contributions and their physical origin. 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 η vary 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×10 -2 and experiences little variation up to Bo ≈ 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 [START_REF] Mougin | Wake-induced forces and torques on a zigzagging/spiralling bubble[END_REF], 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 nonzero imaginary component of the deformation (dashed green line in figure 16(a)). The relative magnitude of the deformation does not exceed 3 × 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-M o liquids. These findings also validate the 'rigid' qualification employed for the primary lowfrequency 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 (resp. imaginary) part of the displacements associated with the bubble inclination (resp. deformation) is typically one order of magnitude smaller than its imaginary (resp. 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. These two weak components increase sharply with Bo and become dominant for Bo 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 0.59. This behaviour indicates that the bubble now exhibits a nonzero 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 nonzero 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 Bo 0.6.

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

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 nonoscillatory 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) and (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 to 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 16(a)), 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 × 10 -2 at Bo = 0.8.

Figure 16(b) shows how the relative magnitude of the three contributions to η 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 M o explored here, while those associated with the bubble inclination are almost in quadrature. Again, we attribute the non-smooth variations of the weak real (resp. imaginary) component of the inclination-induced (resp. deformation-induced) contribution observed for M o 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 ηdef max 0 , follow the approximate power law ηdef max 0 ∼ M o 0.23 , increasing by two orders of magnitude from Galinstan (η def max 0 ≈ 10 -3 ) to DMS-T11 (η def max 0 ≈ 7 × 10 -2 ). The relative inclination-induced displacements only grow by one order of magnitude, from 4 × 10 -2 to 5 × 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.

Summary and final discussion

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 ≈ 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 M o 10 -7 . Conversely, this mode gives birth to two stationary modes if the Morton number is low enough (M o 10 -9 ). In the intermediate range 10 -9 M o 10 -7 , this mode remains oscillatory even far from the primary threshold but its growth rate varies nonmonotonically with the bubble size, and it may even re-stabilise 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 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 rigidbody 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 × 10 -3 for water. It exceeds 1% only for liquids with M o 1.5 × 10 -8 , reaching 7% for DMS-T11 (M o = 9.9 × 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 frozen-shape approximation (FSA) by Cano-Lozano et al. (2016a). 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 M o ≈ 10 -13 to 23% at M o ≈ 6 × 10 -7 , before decreasing sharply to nearly 15% for larger M o. For M o 7 × 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 M o 3×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 × 10 -11 M o 3 × 10 -7 .

Discussion

To better understand the above issue and gain some insight into the subtle balances involved in the various regimes encountered by varying the fluid properties, it is useful to put the first stages of the bifurcation diagrams obtained with deformable bubbles in perspective with those determined by Tchoufag et al. (2014b) who assumed a perfectly oblate spheroidal bubble with a frozen shape. By increasing the geometrical aspect ratio χ in the range [2.21, 2.30] they identified three successive phenomenologies summarized in their figures 1 and 2. For 2.21 χ 2.23, the system first becomes unstable through a low-frequency mode which, by increasing the bubble size (i.e. the Reynolds number in their case), turns into a pair of stationary modes, one of which is quickly damped as the bubble size is increased further while the other gets strongly amplified (their figure 1). Present findings summarized in figure 13 indicate that this is exactly what happens with deformable bubbles rising in water and more generally in low-M o liquids. That the scenario relevant to frozen spheroidal bubbles remains unchanged in the case of freely-deforming bubbles rising in low-M o liquids strongly supports the view that time-dependent deformations do not play any causal role in the mechanisms governing the first stages of path instability in such fluids, although they lower the critical bubble size by 18% for water and increase the reduced frequency by nearly 30%, as may be deduced from the two neutral curves in figures 7 and 8 for M o = 2.54 × 10 -11 . The wake alone is not the source of the instability, since it is still stable at the FSA threshold (Cano-Lozano et al. 2016a). Therefore, as already stated by [START_REF] Bonnefis | When, how, and why the path of an air bubble rising in pure water becomes unstable[END_REF] in the case of water, the key explanation of path instability in low-Morton-number fluids lies in the specific physical processes accounted for in FSA, namely the degrees of freedom allowed by the rigid-body buoyancy-driven motion of the bubble, and the back reaction imposed by the fluid to this motion by the constant-force and zero-torque constraints.

Increasing the bubble aspect ratio to the range 2.23 χ 2.25, Tchoufag et al. (2014b) observed that, although it is still present, the low-frequency mode is no longer unstable, its growth rate keeping weak negative values. In the approximation they used, path instability then arises through a stationary bifurcation, similar to what the FSA neutral curves (pale grey lines) in figures 7 and 8 indicate for 7×10 -11 M o 3×10 -7 . The corresponding stationary mode is that associated with the wake instability observed when the bubble is held fixed in a uniform stream. Hence, one has to conclude that the additional degrees of freedom offered by the free motion of a frozen-shape bubble are not sufficient to change qualitatively the overall dynamics of the system in that range. The yellow neutral curve St = f (M o) in figure 8 proves that, despite their weak relative magnitude (from 0.4% for DMS-T00 to 2% for M o = 3 × 10 -7 according to figure 16(b)), deformations are able to restore the instability of the oscillatory mode and to make it more unstable than the stationary mode, keeping the near-threshold path instability phenomenology unchanged with respect to the one prevailing in the low-M o-range. Nevertheless, expanding the growth rate in the form λ r ≈ α(M o)(Bo -Bo c (M o)), with Bo c the threshold Bond number, and considering variations of α(M o) from water (figure 13(a)) to DMS-T02 (figure 12(a)) via DMS-T00 (not shown), one finds α(M o) ∼ M o -0.15 . Hence, the larger M o the slower the increase of the growth rate of the oscillatory mode with the distance to the threshold, despite the fact that the magnitude of surface deformations at the threshold increases as M o 0.23 (figure 16(b)). This confirms that the damping mechanism predicted by FSA, be it with the schematic spheroidal shape prescribed by Tchoufag et al. (2014b) or with the more realistic shape corresponding to the actual base state, is still active in this intermediate M o-range in the case of deformable bubbles. The reasons for this damping are not obvious. For a non-deforming body, the eigenvalues of the system depend directly on the drag and torque associated with a small edgewise translation or a slow rotation about the major axis perpendicular to the plane of motion [START_REF] Fabre | A quasi-static approach to the stability of the path of heavy bodies falling within a viscous fluid[END_REF]. We believe that when the bubble geometrical anisotropy is varied, small changes in the structure of the base flow, especially in the near wake, make these loads vary in a non-monotonic manner with χ, hence with M o, which makes eventually the real part of the eigenvalues associated with the oscillatory mode be positive in some ranges of the parameter space but slightly negative in others.

Last, for χ 2.25, Tchoufag et al. (2014b) found the first unstable mode to be oscillatory again. Similar to what is noticed on the FSA neutral curve in figure 8, the corresponding oscillations are much faster than those observed for 2.21 χ 2.23, the two oscillatory modes predicted by FSA in the two subdomains being distinct (pale and dark grey lines in figures 7 and 8). Deformation effects totally smooth out this abrupt change by allowing the frequency to increase gradually throughout the intermediate M o-range discussed above. Increasing the Reynolds number at a given χ, i.e. the bubble size in a given fluid, Tchoufag et al. (2014b) found that the next mode that becomes unstable is the stationary mode associated with wake instability. The FSA prediction obtained with the actual base flow confirms this conclusion. For instance, the dark grey line in figure 7 shows that, for DMS-T05, the 'high-frequency' mode becomes unstable at Bo ≈ 5.0. Then the stationary mode (pale grey line) becomes unstable in turn at Bo = 5.85, i.e. for a bubble only 8% larger than at the primary threshold. This sequence subsists when the bubble deforms freely, but the gap between the two thresholds increases considerably, with the oscillatory and stationary modes becoming unstable at Bo = 3.4 and Bo = 5.95, respectively. Moreover, the influence of deformation manifests itself through the emergence of the (2, 0) axisymmetric shape mode that becomes unstable in between these two thresholds. The above comparisons establish that the role of bubble deformations on the most unstable non-axisymmetric mode responsible for path instability deeply differs among the three M o-ranges identified by FSA. This role is limited to mere, albeit significant, quantitative changes in the critical bubble size and oscillations frequency in the low-and 'high'-M o regimes, whereas it induces a complete change in the nature of the most unstable mode in the intermediate-M o regime.

A weakly nonlinear approach inspired by those developed by [START_REF] Fabre | The steady oblique path of buoyancy-driven disks and spheres[END_REF] and [START_REF] Tchoufag | Weakly nonlinear model with exact coefficients for the fluttering and spiraling motion of buoyancy-driven bodies[END_REF] for buoyancy/gravity-driven rigid bodies is desirable to extend present predictions to slightly supercritical bubble sizes. By incorporating nonlinear corrections to the base state resulting from the primary non-axisymmetric unstable mode, this extension may in particular allow the disappearance of the stationary mode at any Morton number to be rationalized. Such a weakly nonlinear approach would also settle the question of the supercritical or subcritical nature of the primary bifurcation in each Morton number range, which the present linear approach cannot answer. Another more modest but much cheaper approach worth exploring is the extension to bubbles (and more generally to rigid bodies with a moderate or low body-to-fluid density ratio) of the 'aerodynamic' reduced-order model derived by [START_REF] Fabre | A quasi-static approach to the stability of the path of heavy bodies falling within a viscous fluid[END_REF] for heavy freely-falling two-dimensional bodies. This model should be able to predict accurately the most unstable or least stable eigenvalues of the system for bubbles with a frozen shape, provided that the drag and torque coefficients associated with a small edgewise translation and a slow rotation of the bubble are computed in the relevant base flow. This is the route to be followed to quantify the influence of base flow variations resulting from slight changes in the interface geometry on the hydrodynamic reaction experienced by bubbles close to the threshold of path instability.

Declaration of Interests. The authors report no conflict of interest.

Figure 1 :

 1 Figure 1: Flow domain and geometrical transformation involved in the L-ALE approach. The black line represents the current bubble-fluid interface Γ b (t) bounding internally the physical fluid domain Ω(t), while the blue line corresponds to the initial interface Γ b0 bounding the reference domain Ω 0 . Both domains exhibit a rotational symmetry about the Γ z axis.

Figure 2 :

 2 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 (M o = 6.2× 10 -7 ); (b): silicone oil DMS-T02 (M o = 1.6 × 10 -8 ); (c): water at 20 • C (M o = 2.54 × 10 -11 ).

Figure 3 :

 3 Figure 3: Rise Reynolds number vs the Bond number in several fluids. Solid lines: present predictions; symbols: experimental data, with : ultrapure water at 20 • C (M o = 2.54 × 10 -11 ) (Duineveld 1995), : silicone oil DMS-T00 (M o = 1.8 × 10 -10 ), : DMS-T02 (M o = 1.6×10 -8 ), : DMS-T05 (M o = 6.2×10 -7 ), and : DMS-T11 (M o = 9.9×10 -6 )(data for all four series taken from[START_REF] Zenit | Path instability of rising spheroidal air bubbles[END_REF]). Red bullets mark the onset of path instability detected experimentally in each series.

Figure 4 :

 4 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.

Figure 5 :

 5 Figure 5: Base flow past bubbles rising in three different liquids. (ac): DMST05 (M o = 6.2• 10 -7 ), with, from left to right, Bo = 3, 5 and 7; (df ): DMST02 (M o = 1.6 × 10 -8 ), with Bo = 1, 2 and 3; (gi): water at 20 • C (M o = 2.54 • 10 -11 ), with Bo = 0.4, 0.6 and 0.8. The left and right halves of each subfigure display the azimuthal vorticity and pressure distributions, respectively. The thin lines are the streamlines in the reference frame rising with the bubble.

Figure 6 :

 6 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 et al. (2013). The left and right inserts 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; M o values are specified along the upper and right sides of the figure.

Figure 7 :

 7 Figure 7: Neutral curve in the (Bo, Ga) plane, obtained by connecting the 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; M o values are specified along the upper and right sides of the figure. In a given fluid, open (resp. closed) symbols indicate stable (resp. unstable) paths observed in the simulations of Cano-Lozano et al. (2016b) ( , ), in the experiments of De Vries (2002) ( , ) and Duineveld (1995) ( , ) in ultrapure water (performed at temperatures of 28 • C and 20 • C, respectively), and in those of Zenit & Magnaudet (2008) ( , ) and Sato (2009) ( , ) in various silicone oils.

Figure 8 :

 8 Figure8: Variation of the reduced frequency St = λ i D/(2πu b ) at the threshold vs the Morton number. The yellow curve was obtained by connecting the 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 'frozen-shape' approximation (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 M o = 2.7 × 10 -7 indicates the frequency jump associated with the switching from the stationary mode to the oscillatory mode predicted by FSA. : numerical predictions ofCano-Lozano et al. (2016b); : experimental data of[START_REF] Duineveld | Bouncing and coalescence of two bubbles in water[END_REF] for Bo = 0.54 in ultrapure water; : experimental data of[START_REF] Zenit | Measurements of the streamwise vorticity in the wake of an oscillating bubble[END_REF] for Bo = 3.92 in DMS-T05.

Figure 9 :

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

Figure 10 :

 10 Figure 10: Variation of the first five unstable eigenvalues with the Bond number for bubbles rising in DMS-T05 (M o = 6.2 × 10 -7 ). (a): growth rate (λ r ), normalized by the gravitational time (D/g) 1/2 ; (b): radian frequency (λ i ) normalized similarly; (c): equilibrium shapes at the threshold (labels refer to the transition points marked with bullets in (a)). In (a) -(b), dashed (resp. solid) lines represent the stable (resp. 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 = π/(2λ i ) (blue dashed line); bubble centroids are arbitrarily shifted in the z direction for readability.

Figure 11 :

 11 Figure 11: Structure of the first unstable mode past a bubble rising in DMS-T05 at Bo = 3.7. (a) and (c): real part of the mode near the bubble and further downstream in the wake, respectively; (b) and (d): same for the imaginary part. The left and right halves in (ab) display the pressure and azimuthal vorticity iso-levels, respectively; those in (cd) 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 (η), respectively.

Figure 12 :

 12 Figure 12: Same as figure 10 for bubbles rising in DMS-T02 (M o = 1.6 × 10 -8 ).

Figure 13 :

 13 Figure 13: Same as figure 10 for bubbles rising in water at 20 • C (M o = 2.54 × 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.

Figure 14 :

 14 Figure 14: Structure of the rigid modes past a bubble rising in water. Top: primary lowfrequency mode at Bo = 0.48 in the bubble vicinity, with the real and imaginary parts shown in (a) and (b), respectively. The left and right halves of the two subfigures 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. Bottom: 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 subfigure displays the real part of the azimuthal vorticity iso-levels; the left half of (c) and (e) (resp. (d)) displays the real (resp. imaginary) part of the azimuthal velocity iso-levels.

Figure 15 :

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

Figure 16 Figure 16 :

 1616 Figure16(a) considers the case of the most amplified rigid mode for a bubble rising in water at 20 • C, i.e. the low-frequency oscillating mode up to Bo ≈ 0.64, followed by the

  ργ 3 )

	Galinstan Water (28 • C, 1 atm) Water (20 • C, 1 atm) DMS-T00 DMS-T02 DMS-T05 DMS-T11	6440 996 1000 761 873 918 935	2.40 0.80 1.00 0.49 1.75 4.59 9.35	718 71.0 72.8 15.9 18.7 19.7 20.1	1.4 × 10 -13 1.13 × 10 -11 2.54 × 10 -11 1.8 × 10 -10 1.6 × 10 -8 6.2 × 10 -7 9.9 × 10 -6

† Email address for correspondence: jmagnaud@imft.fr

Author ORCIDs. P. Bonnefis, https://orcid.org/0000-0003-2825-2974; J. Sierra-Ausin, https://orcid.org/0000-0001-6036-5093; D. Fabre, https://orcid.org/0000-0003-1875-6386; J. Magnaudet, https://orcid.org/0000-0002-6166-4877.