Jie Zhang

Ming-Jiu Ni

Jacques Magnaudet

Three-dimensional dynamics of a pair of deformable bubbles rising initially in line. Part 1: Moderately inertial regimes

Keywords: Bubble dynamics, hydrodynamic interaction, multiphase flow

The buoyancy-driven motion of two identical gas bubbles released in line in a liquid at rest is examined with the help of highly resolved simulations, focusing on moderately inertial regimes in which the path of an isolated bubble is vertical. Assuming first an axisymmetric evolution, equilibrium configurations of the bubble pair are determined as a function of the buoyancy-to-viscous and buoyancy-to-capillary force ratios which define the Galilei (Ga) and Bond (Bo) numbers of the system, respectively. The threedimensional solutions reveal that this axisymmetric equilibrium is actually never reached. Instead, provided Bo stands below a critical Ga-dependent threshold determining the onset of coalescence, two markedly different evolutions are observed. At the lower end of the explored (Ga, Bo)-range, the tandem follows a Drafting-Kissing-Tumbling scenario which eventually yields a planar side-by-side motion. For larger Ga, the trailing bubble drifts laterally and gets out of the wake of the leading bubble, barely altering the path of the latter. In this second scenario, the late configuration is characterized by a significant inclination of the tandem ranging from 30 • to 40 • with respect to the vertical. Bubble deformation has a major influence on the evolution of the system. It controls the magnitude of vortical effects in the wake of the leading bubble, hence the strength of the attractive force acting on the trailing bubble. It also governs the strength and even the sign of the lateral force acting on this bubble, a mechanism of particular importance when the tandem is released with a small angular deviation.

Introduction

Buoyancy-driven bubbly flows are widely encountered in natural environments (e.g. breaking waves, bubbly plumes released from the floor of lakes and oceans) and engineering devices (e.g. bubble columns, ladle steel making, boiling flows in power plants). Such applications have driven fundamental studies of fluid-bubble interactions in bubbly suspensions for a long time. Early computational investigations [START_REF] Smereka | On the motion of bubbles in a periodic box[END_REF][START_REF] Sangani | Dynamic simulations of flows of bubbly liquids at large Reynolds numbers[END_REF]) assumed spherical bubble shapes and neglected any possible influence of vorticity in the liquid. The corresponding potential flow simulations predict the formation of large horizontal bubble clusters. However, subsequent laboratory experiments and simulations based on the full Navier-Stokes equations revealed a less clear-cut picture. When bubbles remain nearly spherical and viscous effects, while smaller than inertial effects, remain significant in the bulk, experiments [START_REF] Cartellier | Bubble-induced agitation and microstructure in uniform bubbly flows at small to moderate particle Reynolds numbers[END_REF] and threedimensional simulations [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 1. Low Reynolds number arrays[END_REF][START_REF] Bunner | Dynamics of homogeneous bubbly flows. Part 1. Rise velocity and microstructure of the bubbles[END_REF][START_REF] Yin | Lattice-Boltzmann simulation of finite Reynolds number buoyancydriven bubbly flows in periodic and wall-bounded domains[END_REF][START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF] show that the microstucture is governed by the pair interaction mechanism known as Drafting-Kissing-Tumbling (hereinafter abbreviated as DKT) for sedimenting solid particles [START_REF] Joseph | Nonlinear mechanics of fluidization of spheres, cylinders and disks in water[END_REF][START_REF] Fortes | Nonlinear mechanics of fluidization of beds of spherical particles[END_REF]. This is primarily a wake effect by which two particles or bubbles initially aligned vertically are first attracted toward each other, then repel in the horizontal direction when they get very close to each other, until they reach an equilibrium separation and fall/rise side by side. The process is self-repeating, since at some point each of the two bodies enters unavoidably in the wake of one of its neighbours. In more inertial regimes, experiments [START_REF] Cartellier | Bubble-induced agitation and microstructure in uniform bubbly flows at small to moderate particle Reynolds numbers[END_REF][START_REF] Zenit | Measurements of the average properties of a suspension of bubbles rising in a vertical channel[END_REF][START_REF] Figueroa-Espinoza | Clustering in high Re monodispersed bubbly flows[END_REF] and simulations [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 1. Moderate Reynolds number arrays[END_REF][START_REF] Yin | Lattice-Boltzmann simulation of finite Reynolds number buoyancydriven bubbly flows in periodic and wall-bounded domains[END_REF] with nearlyspherical bubbles reveal a clear tendency of bubbles to align horizontally. However, the corresponding clusters are less strong, i.e. the bubble distribution is less anisotropic, than predicted by potential flow theory. Simulations also considered effects of bubble deformation. In moderately inertial regimes, the results show that pairs of significantly oblate bubbles tend to align vertically, forming vertical streams [START_REF] Bunner | Effect of bubble deformation on the properties of bubbly flows[END_REF]. However this ' chimney' effect disappears in strongly inertial regimes in which bubbles tend to follow zigzagging or spiralling paths [START_REF] Esmaeeli | A direct numerical simulation study of the buoyant rise of bubbles at O(100) Reynolds number[END_REF].

This brief review highlights the fact that the microstructure of buoyancy-driven bubbly suspensions is to a large extent governed by pair interactions. In particular, the two canonical configurations in which two bubbles are released either in line or side by side are of particular relevance to obtain a better understanding of the local mechanisms at stake in freely-rising suspensions. Configurations corresponding to intermediate initial inclinations connect these two extreme geometries and were considered in the nearlyinviscid limit, both theoretically (assuming a spherical bubble shape) and experimentally, by Kok (1993a,b). Detailed experiments were carried out with bubble pairs rising side by side [START_REF] Duineveld | Bouncing and coalescence of bubble pairs rising at high Reynolds number in pure water or aqueous surfactant solutions[END_REF][START_REF] Sanada | Motion and coalescence of a pair of bubbles rising side by side[END_REF][START_REF] Kong | Hydrodynamic interaction of bubbles rising side-by-side in viscous liquids[END_REF], varying the liquid properties, bubble sizes and initial separation. Depending on flow conditions, the two bubbles were found to repel or attract each other. In the latter case, they may reach an equilibrium separation or collide, in which case they subsequently bounce or coalesce. The respective roles of irrotational and vortical effects on the sign and magnitude of the transverse interaction force were examined in the simulations of [START_REF] Legendre | Hydrodynamic interactions between two spherical bubbles rising side by side in a viscous liquid[END_REF] assuming spherical bubbles. In particular, a regime map predicting the characteristics of the final configuration as a function of the initial separation was obtained. Influence of bubble deformation, which beyond a critical oblateness makes the wake unstable, was considered numerically by [START_REF] Zhang | Vortex interactions between a pair of bubbles rising side by side in ordinary viscous liquids[END_REF]. The resulting double-threaded wakes and their interactions were found to be critically important during the collision stages. Indeed, in most cases this interaction is responsible for an extra repulsive transverse force which makes the two bubbles bounce.

Not surprisingly, the in-line configuration was first considered assuming spherical bubbles and an axisymmetric flow at all time. With an irrotational flow in the bulk supplemented with a weak boundary layer and wake past each bubble, [START_REF] Harper | On bubbles rising in line at large Reynolds numbers[END_REF] established the existence of a finite equilibrium separation of the two bubbles. He showed that this equilibrium stems from the balance between a repulsive force corresponding to the irrotational flow past the two bodies and an attractive force resulting from the influence of the boundary layer past the leading bubble on the boundary layer of the trailing bubble. His conclusions were qualitatively confirmed and extended toward lower Reynolds numbers through axisymmetric computations by [START_REF] Yuan | On the in-line motion of two spherical bubbles in a viscous fluid[END_REF] (and later by [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF] who considered arbitrary orientations of the bubble pair). This investigation prompted Harper to improve his theory by accounting for viscous diffusion in the wake of the leading bubble, allowing him to reach a better agreement with the numerical results for large Reynolds numbers [START_REF] Harper | Bubbles rising in line: why is the first approximation so bad?[END_REF]. Early experiments with nearly spherical bubbles were performed in weaker inertial regimes corresponding to Reynolds numbers lower than those considered by [START_REF] Yuan | On the in-line motion of two spherical bubbles in a viscous fluid[END_REF]. Using distilled water, [START_REF] Katz | Wake-induced relative motion of bubbles rising in line[END_REF] observed that under such conditions the two bubbles always collide and coalesce. Experiments performed in silicone oils by Watanabe & Sanada (2006) in the same regime confirmed that the bubbles collide but revealed no coalescence. In moderately inertial regimes, these authors did not observe any head-on collision, in line with the numerical findings of [START_REF] Yuan | On the in-line motion of two spherical bubbles in a viscous fluid[END_REF]. However, they found the equilibrium axisymmetric configuration to be unstable, confirming Harper's theoretical analysis [START_REF] Harper | On bubbles rising in line at large Reynolds numbers[END_REF]. The three-dimensional evolution of the bubble pair in moderately inertial regimes was explored in more detail by [START_REF] Kusuno | Experimental investigation of the motion of a pair of bubbles at intermediate Reynolds numbers[END_REF] and [START_REF] Kusuno | Lift force acting on a pair of clean bubbles rising in-line[END_REF], using ultrapure water and silicone oil, respectively. Several interaction scenarios were reported, including the DKT process and a distinct evolution (also noticed in the computations of [START_REF] Gumulya | Interaction of bubbles rising inline in quiescent liquid[END_REF]) in which the trailing bubble drifts laterally without significantly modifying the path of the leading bubble. In this case, the separation between the two bubbles remained significantly larger than their radius throughout the rise.

Focusing on the initial in-line configuration, possibly with some small angular deviation, the present investigation aims at providing a more detailed understanding of the interaction processes reviewed above. For this purpose, we carried out high-resolution three-dimensional time-dependent computations allowing a complete interplay of inertial, viscous and capillary effects over a wide range of flow regimes. The present paper reports on the first half of this investigation. It focuses on moderately inertial regimes in which each bubble, taken isolated, would follow a straight vertical path. High-inertia regimes in which isolated bubbles follow a non-straight path will be examined in the companion paper [START_REF] Zhang | Dynamics of a pair of deformable bubbles rising in line. Part 2: Highly inertial regimes[END_REF]). The entire work is based on the open source code Basilisk [START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations[END_REF] which employs the volume of fluid (VOF) approach. This method enables the bubbles to deform freely in the liquid as they rise and interact. Moreover, the adaptive mesh refinement (AMR) technique embedded in this code, supplemented with a specific refinement [START_REF] Zhang | Vortex interactions between a pair of bubbles rising side by side in ordinary viscous liquids[END_REF], makes it possible to properly capture the flow in the vicinity of the bubble surface, in the near wake, as well as within the film that forms between the two bubbles when they get very close to each other. The paper is organized as follows. The problem and the dimensionless parameters are introduced in § 2, while § 3 summarizes the numerical method. Before embarking on the discussion of numerical results, the fundamental mechanisms involved in the interaction process for spherical and deformed bubbles are reviewed in § 4. Predictions obtained by constraining the flow to remain axisymmetric are discussed in § 5. Then, fully three-dimensional evolutions are examined in § 6. Influence of bubble deformation and initial conditions on the evolution of the bubble pair is discussed in § 7. The main findings and some open issues are summarized in § 8.

Problem statement

We consider a pair of deformable gas bubbles rising in line in a large expanse of liquid. The bubbles are assumed to have the same volume V, hence the same equivalent radius R = (3V/4π) 1/3 . Initially spherical, they are released from rest near the bottom of the numerical tank, with their line of centres vertical (Y -direction) and their centres separated by a distance S 0 . The initial configuration is illustrated in figure 1(a). The three-dimensional computational domain is cubic, with a size of (240R) 3 , which makes it large enough for minimizing artificial confinement effects. The two bubbles start to rise simultaneously under the effect of buoyancy, which is somewhat different from experimental studies in which they are usually released in sequence. During the rise, deviations of the line of centres of the bubble pair from the vertical is characterized by the angle θ(t) (figure 1(b)).

The liquid and bubble motions are governed by the incompressible one-fluid Navier-Stokes equations

∇ • u = 0 , ρ(∂ t u + u • ∇u) = (ρ -ρ)g -∇p + ∇ • Σ + F s .
(2.1)

In (2.1), F s = γκδ s n stands for the capillary force, with γ the surface tension, n the unit normal to the interface, κ = -∇ • n the interface mean curvature, and δ s the Dirac function identifying the interface position. The suspending liquid and the gas within the bubbles being assumed Newtonian, the viscous stress tensor reads Σ = µ(∇u + ∇u T ), with µ the dynamic viscosity and the superscript T standing for the transpose operator.

The density ρ and viscosity µ are uniform in both the liquid and the gas and experience a jump at the interface. A free-slip condition is imposed on all four lateral boundaries, while a periodic condition is assumed to hold on the top and bottom boundaries. In order to ensure that the net momentum flux through the bottom and top planes is constant, and to prevent gravity from accelerating the flow in the vertical direction, the gravity force ρg is supplemented by a body force -ρg with ρ = f ρ g + (1 -f )ρ l , f denoting the global volume fraction of gas in the computational domain, and ρ l and ρ g the liquid and gas densities, respectively [START_REF] Bunner | Dynamics of homogeneous bubbly flows. Part 1. Rise velocity and microstructure of the bubbles[END_REF].

In addition to the gas/liquid density and viscosity ratios, usually very small, the dynamics of the system is governed by three independent dimensionless numbers, among which the dimensionless initial separation S 0 = S 0 /R. The other two control parameters we vary Ga and Bo in the range 10 < Ga < 30 and 0.02 < Bo < 2.0, respectively. In this parameter range, an isolated bubble follows a rectilinear path, i.e. path instability which is commonly observed for millimeter-size bubbles rising in water does not take place.

For such an isolated bubble, selecting Ga and Bo in the above range leads to terminal Reynolds numbers in the range 10 Re 120, depending on bubble deformation. In most of the present work, the initial separation between the two bubble centres is set to S 0 = 8, hence the initial gap between the two bubbles is 6R. In what follows, all variables are normalized using R and R/g as characteristic length and time scales, respectively. The bubble deformation will often be characterized using the aspect ratio χ = b/a, where b and a denote the length of the major and minor axes, respectively. Note that in most available studies, especially the computational works of [START_REF] Yuan | On the in-line motion of two spherical bubbles in a viscous fluid[END_REF] and [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF], the Reynolds number is based on the bubble diameter rather than the radius. When used for comparison, the corresponding results were converted accordingly.

Numerical approach

General features

The results to be discussed below are obtained by solving (2.1) with the open source flow solver Basilisk developed by Popinet [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations[END_REF]. Basilisk (see basilisk.fr) is the successor of Gerris (http://gfs.sourceforge.net) which has been widely employed over the last fifteen years in detailed explorations of interfacial flows. It makes use of Cartesian grids with a collocated discretization of the velocity and scalar fields. The temporal discretization is based on a second-order fractional step method. In particular, the Godunov-type unsplitted upwind scheme developed by [START_REF] Bell | A second-order projection method for the incompressible navier-stokes equations[END_REF] is used to discretize the advection term, and a fully implicit scheme is used to compute the viscous term. A second-order projection method is employed to ensure that the computed velocity field is divergence-free at the end of each time step. Interfaces are tracked and geometrically reconstructed by a VOF approach in which an accurate well-balanced height function method is used to calculate the interface curvature [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Popinet | Numerical models of surface tension[END_REF]). An AMR technique allows to locally refine the grid close to interfaces and high-vorticity regions, based on a wavelet decomposition of the gas volume fraction and velocity fields, respectively [START_REF] Van Hooft | Towards adaptive grids for atmospheric boundary-layer simulations[END_REF]. This strategy greatly enhances the computational efficiency while guaranteeing a high numerical accuracy in flow regions where subtle physical phenomena are likely to take place. In the present study, the spatial resolution is refined down to ∆ min = R/68 close to the bubble interface. Hence, at the highest Reynolds number considered here (Re ≈ 120), approximately 6 grid points lie within the boundary layer whose thickness is estimated to δ b ∼ Re -1/2 . Figure 2 shows how a typical grid is refined in the vicinity of the two bubbles. In addition to the near-interface zones where ∆ = ∆ min = R/68, refined regions include the boundary layer and wake of each bubble, where the local cell size is ∆ = R/17. The grid coarsens drastically beyond the region displayed in the figure, the largest cells in the far field corresponding to ∆ ≈ 7.5R. Hence the ratio between the largest and smallest cells in the whole domain is 2 9 = 512. Previous works have established the capability of Basilisk to accurately simulate the dynamics of isolated rising bubbles. Moreover, since Basilisk succeeded Gerris and essentially makes use of the same algorithms, the numerous validations performed with Gerris over the years also hold for Basilisk. For instance, Popinet (2017) reproduced successfully with Basilisk the results obtained with Gerris by Cano-Lozano et al. (2016a) concerning the transition from rectilinear to zigzagging or spiralling paths of air bubbles rising in various liquids. In § 5 we shall compare present predictions for the equilibrium distance of two nearly-spherical bubbles (Bo 1) rising in line with the findings of [START_REF] Yuan | On the in-line motion of two spherical bubbles in a viscous fluid[END_REF] and [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF]. This comparison will establish the relevance and accuracy of the present approach in the axisymmetric case. Nevertheless, this configuration will be proved to be unstable in § 6. Therefore, it is necessary to determine to which extent the fate of the three-dimensional system depends on numerical details driving the onset of this instability. This is the purpose of appendix A in which several specific tests are reported.

Most of the computations were run on a personal computer with 24 Intel R Xeon R E5-2630 v3 processors. The Intel-MPI library was used to exchange informations between the processors. A typical run extending over 50 time units took approximately 50 days (e.g. the case Ga = 30, Bo = 0.3 discussed in § 6.3). Note that, since capillary effects impose a specific time step constraint, low-Bo cases require longer computational times. For instance, to reach a given physical time, the case Ga = 30, Bo = 0.05 was 1.5 more time consuming than the previous case. Due to the AMR procedure, the grid evolves over time. For Ga = 30, Bo = 0.3, its size stabilizes at approximately 2.5 million cells during the second half of the run. When additional levels of grid refinement are introduced because the bubbles get very close to each other (see below), this size increases significantly. In such cases, e.g. Ga = 20, Bo = 0.5 discussed in § 6.4, the complete grid involves up to 4.2 million cells.

Treatment of thin films: numerical vs. physical coalescence

In the present problem it is likely that under certain conditions the two bubbles get very close to each other. When this happens in a real flow, coalescence may or may not take place, depending on the mobility of the interfaces involved and on the strength and duration of the forces that drive the two bubbles toward each other [START_REF] Chesters | The modelling of coalescence processes in fluid-liquid dispersions: a review of current understanding[END_REF][START_REF] Chan | Film drainage and coalescence between deformable drops and bubbles[END_REF]. Dealing numerically with such situations is particularly challenging, owing to the very small scales involved. Some numerical approaches, especially the Front Tracking technique [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF]; [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF]), totally prevent coalescence. The same may be achieved in VOF approaches by identifying the two bubbles with separate markers, each of them representing the local volume fraction of the corresponding body. However, this numerical option is not fully appropriate here. Indeed, for the reasons mentioned below, bubbles rising in line in a pure liquid offer one of the physical situations with the highest coalescence probability. If coalescence takes place under real conditions, it is of course desirable to track numerically the post-coalescence dynamics, i.e. the shape and path evolution of the resulting bubble, which the above option would not allow.

That bubbles rising in line in a pure liquid are prone to coalesce in a number of cases is due to the combination of two factors. First, as will be discussed in the next section, the wake of the leading bubble provides a permanent attractive force to the trailing bubble. Under a number of flow conditions, this force is strong enough to make the two bubbles come virtually in contact in the head-on configuration. Second, the mobility of the gasliquid interfaces when the bubbles are free of any contamination makes the drainage of the interstitial film several orders of magnitude faster than that of liquid-liquid or contaminated gas-liquid interfaces [START_REF] Vakarelski | Coalescence dynamics of mobile and immobile fluid interfaces[END_REF]. The coalescence process may take several distinct forms, depending on the fluid characteristics and bubble size. After the two bubbles collide, the interstitial film may rupture quickly, or the bubbles may stay almost in contact during a long time before eventually coalescing. If viscous effects are small enough, bubbles larger than a critical size bounce after the film has been drained partially and coalesce only after one or several bounces. Energetic considerations and detailed experimental data may be employed to determine which of these scenarios takes place. This issue, together with the specificities of the coalescence of 'clean' bubbles are discussed in appendix B.

In situations potentially leading to coalescence, one would ideally like to track numerically all steps of the drainage, until non-hydrodynamic effects such as the London-van der Waals force come into play and rupture the film. This typically occurs when the minimum gap between the two interfaces is of the order of 10 nm. For millimeter-size bubbles, this would require approximately ten additional grid levels beyond the one corresponding to ∆ min . Since the largest scales to be resolved in the present context are of the order of 10 cm, i.e. seven orders of magnitude larger, the corresponding computational cost would be prohibitive. At least, it is possible to track the first stages of the drainage based on a suitable grid refinement technique. Then, referring to the available knowledge summarized in appendix B, one can reasonably predict which coalescence scenario is taking place in the real system, although the grid resolution does not allow all its details to be captured. The main shortcoming of this approach is that bubbles usually merge too early in the computations, i.e. they would merge at a higher vertical position in a real flow. Nevertheless, as we discuss in appendix B, the theoretical and experimental knowledge available on the coalescence of nearly-spherical clean bubbles allows the corresponding temporal shift to be estimated in a number of cases.

To track the first steps of the drainage, we developed a specific topology-based AMR scheme to refine the grid within the gap [START_REF] Zhang | Vortex interactions between a pair of bubbles rising side by side in ordinary viscous liquids[END_REF]. The corresponding algorithm checks whether or not any cell crossed by the gas-liquid interface has at least one neighbouring cell filled with only liquid or gas. If not, the cell crossed by the interface is automatically refined. This adaptive strategy is illustrated in figure 3. At t = t 1 (panel (a)), the two interfaces are separated by a 'pure' liquid cell, but this is no longer the case in (b) where the cells standing in the yellow region all contain a nonzero gas fraction. Therefore, these cells are refined by a factor of 2 in each direction, as shown in (c).

Ultimately, for the reasons discussed above, we let the two interfaces merge numerically if the number N of successive refinements required to satisfy the above 'pure liquid neighbouring cell' criteria exceeds a prescribed value. In practice, this is achieved by merging the two markers which identify each bubble in all previous steps. The threedimensional results to be discussed later were all obtained with N = 2, so that the minimum cell size in the gap was δ min = ∆ min /2 N = R/272. For a 0.25 mm-radius bubble, this implies δ min ≈ 1 µm, which is still one hundred times larger than the typical thickness at which the interstitial film actually ruptures.

Fundamental mechanisms: role of fluid inertia, wake effects and bubble deformation

Before analyzing the computational results, a brief review of the main physical mechanisms involved in the problem is in order. First of all, it is key to keep in mind that finite-Reynolds-number interactions between two bubbles are controlled by two antagonistic effects. The first of them is due to the outer flow past the bubble pair, in which irrotational mechanisms prevail. In the potential flow approximation, exact results for the fluid kinetic energy have been established for two spherical bubbles having identical radii, from which the interaction force may be obtained [START_REF] Voinov | On the motion of two spheres in a perfect fluid[END_REF][START_REF] Voinov | Hydrodynamic interaction of bodies in an ideal incompressible liquid and their movement in inhomogeneous flows[END_REF][START_REF] Van Wijngaarden | Hydrodynamic interactions between gas bubbles in liquid[END_REF][START_REF] Miloh | Hydrodynamics of deformable contiguous spherical shapes in an incompressible inviscid fluid[END_REF][START_REF] Bentwich | On the exact solution for the two-sphere problem in axisymmetrical potential flow[END_REF][START_REF] Biesheuvel | The motion of pairs of gas bubbles in a perfect liquid[END_REF]Kok 1993a). These predictions indicate that the interaction is repulsive when the two bubbles rise in line because the fluid velocity reaches a minimum in the gap, inducing a pressure maximum there [START_REF] Harper | On bubbles rising in line at large Reynolds numbers[END_REF]. Conversely, the interaction is attractive when the bubbles rise side by side, owing to the flow acceleration (hence the pressure minimum) in the gap. The critical angle at which the interaction force changes sign depends on the separation between the two bubbles, ranging from 35 • when they are in contact to a value close to 54 • when they are far away from each other (Kok 1993a).

Finite-Reynolds-number effects manifest themselves in the generation of vorticity at the bubble surface, owing to the shear-free condition obeyed by the carrying liquid when the gas-to-liquid viscosity ratio is negligibly small and the interface is uncontaminated by surfactants. Diffusion and advection of this vorticity in the surrounding fluid results in a boundary layer and a wake past each bubble. Vortical effects at the bubble surface lower the pressure at the rear stagnation point compared to that at the front [START_REF] Kang | The drag coefficient for a spherical bubble in a uniform streaming flow[END_REF], so that the pressure at a given position along the wake axis is lower than it would be in the potential flow limit. When the bubbles rise in line, this pressure drop makes the trailing bubble (hereinafter abbreviated as TB) sucked toward the leading bubble (hereinafter abbreviated as LB). In contrast, when they rise side by side, the interaction of the two wakes results in a pressure maximum in the gap, hence a force tending to move the two bubbles away from each other.

In a given geometrical configuration of the tandem, the relative magnitude of the above two antagonistic effects depends on the Reynolds number, assuming the bubbles to keep a spherical shape. Vortical effects dominate when Re is low enough, while the interaction is expected to become close to potential flow predictions for large enough Re. For this reason, keeping the Reynolds number and the inclination of the line of centres fixed, the overall interaction force vanishes when the separation between the two bubbles takes a specific value, S e (Re, θ). The smaller Re is, the shorter (resp. larger) this equilibrium separation is for θ = 0 • (resp. θ = 90 • ). The two bubbles collide if S e is small enough, i.e S e 2 for spherical bubbles. Then, assuming that the interface is uncontaminated, they may either coalesce, bounce or stay in contact for a very long time, depending on whether or not the net attractive force is large enough to achieve the drainage of the interstitial film. Approximate models have been proposed to predict S e in the in-line configuration for spherical bubbles. These models consider that the wake of the LB, taken into account by using Oseen or high-Re far-wake velocity distributions, decreases the fluid vertical velocity 'felt' by the TB at a given position by an amount equal to the cross-sectional average of the velocity defect at that position [START_REF] Katz | Wake-induced relative motion of bubbles rising in line[END_REF][START_REF] Ramírez-Muñoz | Drag force on interacting spherical bubbles rising in-line at large Reynolds number[END_REF], 2013).

When the TB moves behind the LB with some offset from the wake axis it faces a nonuniform, asymmetric flow which may locally be considered as a shear flow. A spherical bubble rising in a linear shear flow is known to experience a shear-induced lift force [START_REF] Auton | The lift force on a spherical body in a rotational flow[END_REF][START_REF] Legendre | The lift force on a spherical bubble in a viscous linear shear flow[END_REF]. If the bubble moves faster than the fluid along the streamlines of the base flow, this sideways force deviates it toward the direction of the descending fluid. In the in-line configuration, the relative flow faced by the TB moves downwards and its velocity grows with the distance to the wake axis. Therefore the sideways force tends to move the TB out of the wake of the LB, making the in-line configuration unstable with respect to an infinitesimal lateral deviation [START_REF] Harper | On bubbles rising in line at large Reynolds numbers[END_REF]). As will be seen later, this instability plays a crucial role in the evolution of a bubble pair.

The mechanisms reviewed so far subsist of course when bubble deformation becomes significant. However their magnitude is deeply influenced by the bubble shape. In particular, a key feature of vorticity generation on a curved shear-free interface is that the resulting tangential vorticity, say ω s , is proportional to the product of the local surface curvature and tangential velocity of the fluid [START_REF] Batchelor | An introduction to Fluid Dynamics[END_REF]. This makes the magnitude of ω s increase with bubble deformation. In contrast, for a given interface shape, ω s does not depend on Re when the Reynolds number is large, unlike the more familiar case of a no-slip surface. In the limit Re 1, the maximum of ω s (normalized by the rise velocity and equivalent bubble radius) at the surface of an oblate bubble with an aspect ratio χ increases by a factor of 4 from χ = 1 (spherical bubble) to χ = 2, eventually growing as χ 8/3 when χ 1 [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF]). In the in-line configuration, this dramatic increase implies that the attraction of the TB toward the LB becomes increasingly strong as the latter deforms. Consequently, bubble deformation is expected to reduce significantly the equilibrium separation, favouring coalescence. Another consequence of deformation is its influence on the magnitude and even the sign of the sideways force acting on the TB when the axial symmetry of the in-line configuration is broken by some lateral disturbance. The corresponding mechanisms are discussed in appendix C. Red bubbles maintain a finite separation, while blue bubbles are about to coalesce.

Axisymmetric configuration

Numerical simulations of the in-line configuration based on boundary-fitted grids were reported by [START_REF] Yuan | On the in-line motion of two spherical bubbles in a viscous fluid[END_REF] and [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF] for spherical bubbles. Here in contrast, the bubbles deform freely and continuously while rising, which allows us to investigate the influence of their deformation on the interaction process. In this section, we restrict the generality of the problem by constraining the bubbles to follow a straight vertical path. For this purpose, instead of the cubic domain described earlier, we use an axisymmetric domain whose radius and height are both 240R. Bubbles are released on the symmetry axis are rise along it. Computations are carried out with gas/liquid density and viscosity ratios set to 10 -3 and 10 -2 , respectively.

Figure 4 displays the final bubble shapes and relative positions obtained through 32 computational runs covering the domain (10 Ga 30, 0.02 Bo 1.0), all with S 0 = 8. Red bubbles maintain an equilibrium distance which is finite in the computational sense, i.e. the gap that separates them at steady state exceeds the minimum cell size δ min allowed by the specific topology-based AMR treatment described in § 3.2. Conversely, blue bubbles are such that the 'final' gap is thinner than δ min , implying that coalescence is about to take place numerically. Consider a given row in the figure, i.e. a given Ga. Increasing Bo increases the ability of the bubbles to deform. Hence it reduces the final equilibrium separation distance S e between their centroids, which makes coalescence more likely to occur. For Ga = 10, increasing the Bond number from Bo = 0.02 to Bo = 0.1 makes the Weber number increase from W e = 0.18 to W e = 0.90. As a consequence, the final deformation of the LB and that of the TB rise from χ = 1.02 to χ = 1.1 and from χ = 1.01 to χ = 1.04, respectively. Although still modest, this deformation makes the bottom region of the LB significantly flatter than the top region of the TB, allowing a thin-gap region with a finite area to develop. However, this deformation is small enough for the rise speed to remain virtually unaffected, which leaves the terminal Reynolds number unchanged throughout this range of Bo (Re ≈ 21).

That the LB deforms more than the TB becomes clearer as Ga increases. This is directly due to the fact that the wake of the former reduces the pressure at the front stagnation point of the latter, a mechanism often referred to as the 'sheltering' effect. Hence the pressure difference between the front stagnation point and the bubble equatorial plane, which drives the deformation [START_REF] Moore | The rise of a gas bubble in a viscous liquid[END_REF][START_REF] Moore | The velocity of rise of distorted gas bubbles in a liquid of small viscosity[END_REF], is smaller on the TB. For Ga = 30, the final separation distance decreases from S e = 5.72 to S e = 3.63 when the Bond number increases from 0.05 to 0.15. The Weber number is now of O(1), increasing from 0.65 to 1.44. At the same time, the Reynolds number decreases from 108 to 91, due to the drag increase associated with the increasing oblateness of the two bubbles. The variation of S e with Re and W e may be obtained by considering the bubble pairs of figure 4 once they have reached their final configuration. The result is displayed in figure 5. Specific runs were carried out with Bo = 0.005 (hollow circles in figure 5(a)) to maintain the interface shape very close to a sphere, in order to compare present predictions with results available for spherical bubbles. For (Ga, Bo) = (30, 0.005), the final aspect ratio of the LB (TB) is χ = 1.023 (χ = 1.017), and the deformation is even less for lower Ga. Therefore, the corresponding equilibrium distance is expected to agree well with the correlation proposed by [START_REF] Yuan | On the in-line motion of two spherical bubbles in a viscous fluid[END_REF] As Bo increases, the equilibrium separation quickly falls below that predicted by (5.1) and a correlation in the form S e = f (Re, W e) must be sought. A first attempt aimed at recovering (5.1) in the limit W e → 0 was not successful. After several trials, we found that the entire set of present results is best approached by the fit S e (Re, W e) = 2.025 log Re -3.56 -0.98 W e -0.36 W e 2 .

(5.2)

Figure 5(a) indicates that, once the Weber number has been properly eliminated (W e = (Re/Ga) 2 Bo), all numerical data collapse onto the corresponding curves (dotted lines).

Similarly, figure 5(b) in which Re has been eliminated for the benefit of Ga and Bo confirms that numerical data all follow the W e-dependence defined by (5.2). Therefore (5.2) is seen to provide a valid prediction of the equilibrium separation distance at least in the range 20 Re 120 and 0 < W e 1.5. Actually, we also performed some computations for Ga = 40 and Ga = 50 and found that the equilibrium separation obtained in the range 0.02 Bo 0.2 is still correctly predicted by (5.2). Note that, although (5.2) does not reduce to (5.1) in the limit W e → 0, the corresponding fit (red dash-dotted line in figure 5(a)) achieves a better agreement than (5.1) with present low-Bo predictions, as well as with the results of [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF]. Nevertheless, it must be kept in mind that neither (5.1) nor the low-W e limit of (5.2) remains valid for Re 15, since both expressions predict S e < 2 at lower Re. In the low-but-finite Reynolds number range, asymptotic results for rigid spheres [START_REF] Happel | Low Reynolds Number Hydrodynamics[END_REF] may readily be transposed to bubbles using the scaling argument developed by [START_REF] Legendre | A note on the lift force on a spherical bubble or drop in a low-Reynolds-number shear flow[END_REF]. By doing so, it is concluded that the two bubbles always collide for Re 1, since the LB experiences a larger drag than the TB. Actually, there are computational indications that collision takes place as soon as Re 15.5 (Watanabe & Sanada 2006), which corresponds well to the threshold predicted by (5.2) (S e (Re, W e = 0) = 2 for Re = 15.55).

Three-dimensional configurations

Overview of the results

As mentioned earlier, it is known since [START_REF] Harper | On bubbles rising in line at large Reynolds numbers[END_REF] that the in-line configuration of two clean spherical bubbles is unstable with respect to side disturbances when the Reynolds number is large. Experimental investigations carried out under surfactant-free conditions (silicone oils) in the range 10 Re 150 confirm this prediction [START_REF] Kusuno | Experimental investigation of the motion of a pair of bubbles at intermediate Reynolds numbers[END_REF][START_REF] Kusuno | Lift force acting on a pair of clean bubbles rising in-line[END_REF]. We performed fully three-dimensional time-dependent simulations to assess this stability issue and explore its consequences. Similar to the axisymmetric case, the gas/liquid density and viscosity ratios were set to 10 -3 and 10 -2 , respectively. In agreement with the above experimental and theoretical findings, we observed that the bubble pair never maintains a straight vertical path except when the Bond number is large enough for coalescence to eventually happen. We actually identified three drastically different evolutions, depending on the value of the Galilei and Bond numbers. Beyond a critical Bond number, Bo c (Ga) which increases approximately from 0.2 for Ga = 10 to 0.5 for Ga = 30, the two bubbles collide and eventually coalesce, this 'coalescence' having to be interpreted in the light of the discussion of § 3.2 (see below). Predictions of the three-dimensional and axisymmetric simulations superimpose for Bo 0.5 when Ga 20 (and for Bo 1 when Ga = 30), indicating that the in-line configuration is stable with respect to azimuthal disturbances for sufficiently deformed bubbles. In contrast, for Bo < Bo c (Ga), the TB escapes from the wake of the LB at some point, and the two go on rising with their line of centres more or less inclined with respect to the vertical and their centroids widely separated. In such cases, we found that the interplay of the two bubbles after the three-dimensional effects set in may follow two markedly different scenarios. One is clearly a DKT-type mechanism. In this case, both bubbles deviate from their initial trajectory and eventually rise almost side by side along two straight vertical lines distinct from the initial path. In contrast, in the other scenario, which we refer to as Asymmetric Side Escape (hereinafter abbreviated as ASE), the lateral drift of the TB leaves the path of the LB almost unaffected. Hence this bubble goes on rising virtually along its initial path, while after the system has reorganized itself, the TB follows a markedly distinct vertical path. The structural differences between the configurations corresponding to the DKT and ASE scenarios are well visible in the recent observations of [START_REF] Kusuno | Lift force acting on a pair of clean bubbles rising in-line[END_REF]; see especially their figures 4 and 5.

Figure 6 displays the phase diagrams and typical paths corresponding to the above three evolutions, still for S 0 = 8. The influence of initial conditions, i.e. angular inclination and separation, on the borders of the different subdomains will be discussed in § § 7.3 and 7.4, respectively. The influence of numerical parameters, among which the grid resolution, on the coalescence threshold, i.e. on Bo c (Ga), is estimated in appendix A. These technical aspects are found to barely change Bo c (Ga) by a few percent. As the thin black lines indicate, coalescence only takes place in the present configuration in liquids with a sufficiently high Morton number, from M o ≈ 8 × 10 -7 for Ga = 10 to M o ≈ 1.5 × 10 -7 for Ga = 30. This means in particular that bubbles rising in pure water (M o ≈ 2.6 × 10 -11 ) never coalesce in the Ga-range considered here. Interestingly, figures 6(a) -(b) also reveal that under some conditions such as (Ga, Bo) = (20, 0.2) or (30, 0.2), bubbles that would coalesce if the system were constrained to remain axisymmetric actually escape from coalescence through the ASE scenario. From figure Figure 7. Snapshots of the bubble shapes and relative positions observed during the lateral escape of the TB or the pre-coalescence process. Red, green and blue pairs correspond to the ASE, DKT, and coalescence scenarios, respectively. In the red and green regions, the snapshots are taken by the time the horizontal distance between the two centroids is approximately equal to one bubble initial radius; in the blue region, each snapshot is the last in a run before numerical coalescence occurs. Since successive snapshots are separated by a finite time interval, the remaining time until coalescence may differ among the various blue pairs. 6(b), it may also be concluded that for Re ≈ 15, even an 8% departure from sphericity is sufficient to lead to coalescence. Actually, the discussion of appendix B indicates that, for such modest Reynolds numbers, viscous effects tremendously delay the drainage of the interstitial film. This is why experiments performed in this regime [START_REF] Sanada | Interactions and coalescence of bubbles in stagnant liquid[END_REF]Watanabe & Sanada 2006) reveal that, after the two bubbles collide, they merely stay 'glued' to each other and rise as a single 'dumbbell' bubble without coalescing during the time window of the observations. However, since all forces involved in the physical system, including the London-van der Waals force, are attractive, there is no doubt that such bubbles eventually coalesce. Unfortunately, the spatial resolution of the present simulations is clearly insufficient to reproduce this slow coalescence process. That is, the computations correctly predict the final state of the physical system but this state is reached too early. This underestimate of the coalescence time subsists at larger Reynolds number (see § 6.4) but, for reasons discussed in appendix B, reduces as Re increases. In figure 6(b), the maximum deformation below which the two bubbles do not coalesce is seen to increase significantly with Re. For instance, bubbles with χ ≈ 1.5 follow an ASE scenario for Re ≈ 70, but coalesce for Re ≈ 35 (in these estimates, χ and Re are evaluated from the steady-state properties of the isolated bubble corresponding to the same (Ga, Bo) set).

Figure 7 depicts the bubble shapes and relative positions during the lateral escape of the TB or the coalescence process. Considering the DKT (green) and ASE (red) regimes at a fixed Ga, the figure indicates that the larger Bo the shorter the separation during the escape stage. This may be interpreted as a stabilizing effect of the deformation, since the bubble pair is able to maintain a vertical path during a longer time (i.e. down to a shorter separation) when the Bond number increases. Then, coalescence takes place when Bo exceeds the critical threshold Bo c (Ga) which is seen to increase significantly with Ga, as already noticed in figure 6(a). In most cases, coalescence is reached through a head-on approach, i.e. without any prior escape of the TB, the axisymmetric film in the gap being gradually squeezed. However, for Ga = 30 and Bo = 0.5, an intermediate configuration corresponding to an oblique approach is noticed. These different regimes are examined in more detail in the rest of this section.

Drafting-Kissing-Tumbling

The DKT scenario followed by sedimenting particle pairs has been widely described for spheres [START_REF] Joseph | Nonlinear mechanics of fluidization of spheres, cylinders and disks in water[END_REF][START_REF] Fortes | Nonlinear mechanics of fluidization of beds of spherical particles[END_REF][START_REF] Feng | Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid. Part 1. Sedimentation[END_REF], thick disks [START_REF] Brosse | Interaction of two axisymmetric bodies falling in tandem at moderate Reynolds numbers[END_REF], and under certain conditions prolate spheroids [START_REF] Ardekani | Numerical study of the sedimentation of spheroidal particles[END_REF]. In short, owing to the sheltering effect induced by the wake of the leading body, the trailing body first catches up with it (drafting) until the two collide (kissing). Then the resulting prolate compound body becomes unstable to transverse disturbances, which makes it rotate (tumbling) in such a way that the line joining the centroids of the two individual bodies tends to become horizontal, letting them eventually fall/rise separately in a sideby-side configuration. A similar behaviour of bubbles rising in line has been reported experimentally by [START_REF] Kusuno | Lift force acting on a pair of clean bubbles rising in-line[END_REF] in the range 10 Re 25, 0.3 W e 1.1. Here we identified it for bubble pairs corresponding to Ga 12 and Bo 0.2. With (Ga, Bo) = (10, 0.02), i.e. M o = 8 × 10 -10 , the final state corresponds to Re ≈ 16 and χ ≈ 1.02, which indicates that the bubbles shape remains close to a sphere. The corresponding evolution is displayed in figure 8.

As revealed in panel (a), the tumbling process starts at t 37 when the two bubbles reach the same rise velocity (see figure 9(a)). After this process is completed (say at t = 54 in figure 8(a)), both paths have experienced large lateral deviations. Their lateral separation is about 2.3, the LB having been slightly more deviated. They stand virtually side by side and keep on rising in this configuration, although the rise velocity of the (formerly) TB is barely larger than that of the LB, inducing a tiny difference in the final altitude reached by the two bubbles. During the side-by-side rise, the lateral separation S r evolves significantly, until an equilibrium value S re 2.55 is reached. This is consistent with the findings of [START_REF] Legendre | Hydrodynamic interactions between two spherical bubbles rising side by side in a viscous liquid[END_REF] which indicate that in this configuration, the transverse force acting on a pair of spherical bubbles almost vanishes for S r = 2.5 when Re = 25 (their figure 13). For Y > 175, the centre of inertia of the bubble tandem is seen to drift slightly towards the left, a consequence of the tiny inclination of its line of centres. The sign of this drift is consistent with the conclusions of [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF] who found (their figure 9(b)) that in the range 10 Re 25, the centre of inertia drifts laterally toward the position of the higher bubble, provided 40 • θ < 90 • . The horizontal trace of the TB and LB paths in figure 8(a) indicates that the entire motion remains planar and takes place within a vertical plane. Hence, by breaking the initial axial symmetry of the flow, the DKT mechanism changes the initial one-dimensional path into a two-dimensional path, but the latter keeps track of the former since the preferential direction of the bubble motion remains unchanged. Some three-dimensional streamlines past the two bubbles are displayed in figure 8(b). No standing eddy is observed at the back of the bubbles. This is in line with the generating mechanism of vorticity on a curved shear-free surface discussed in § 4, which maintains the flow past a shear-free sphere unseparated whatever the Reynolds number [START_REF] Blanco | The structure of the axisymmetric high-reynolds number flow around an ellipsoidal bubble of fixed shape[END_REF][START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF]. At t = 38, the separation has almost reached its minimum (S ≈ 2.3) and the flow is still almost axisymmetric, except within the gap where a small left-right asymmetry is discernible in the enlarged view. In the present case, the origin of the axial symmetry breaking stands in tiny numerical asymmetries, in the first place those resulting from the pressure field returned by the multilevel Poisson solver (see appendix A). Tumbling then starts. At t = 44 and t = 48, the flow field exhibits strong left-right asymmetries which result in transverse forces that act to move the two bubbles in opposite directions. Later, say for t 54, the flow field gradually rearranges towards a left-right symmetric configuration corresponding to a sideby-side motion of the two bubbles, with an equilibrium horizontal separation S re ≈ 2.55.

The evolution of the rise speed of the two bubbles is plotted in figure 9(a). The threedimensional prediction is seen to coincide with its axisymmetric counterpart up to t ≈ 45, i.e. throughout the time period during which the bubble pair moves in straight line and even during the early stage of the tumbling process. That the three-dimensional and axisymmetric predictions virtually superimpose up to t = 37 proves the reliability of the former. Note that the red and blue curves superimpose up to t = 4, which corresponds to the very early stage during which the two bubbles rise independently. During the next stage, say 4 < t < 20, the TB goes on accelerating while the rise speed of the LB (hereinafter denoted as V LB ) stays almost constant. This corresponds to the early stage of the interaction during which the trailing bubble is sucked in the wake of the LB while the latter remains unaffected. Then, from t = 20 to t = 37, the two bubbles get close enough for the TB to modify the wake of the LB, the rise speed of which increases sharply until the two bubbles rise with the same speed. Tumbling starts at t ≈ 38, without any real 'kiss' since figure 9(c) indicates that S never fell below 2.3 at previous times. Up to t = 46, no change is noticed in the rise speed of the two bubbles, nor in their vertical separation, although the tumbling process is going on, as the sharp rise of their lateral velocities (figure 9(b)) and that of the inclination of their line of centres (figure 9(d)) confirm. The situation significantly changes within the next short stage 46 t 50, during which V LB and V T B (the rise speed of the TB) drop sharply and the line of centres rotates by more than 30 • . Not surprisingly, this rotation forces the lateral separation to increase dramatically, from S r ≈ 1 at t = 46 to S r ≈ 2.5 at t = 50. Conversely, the vertical separation is reduced to S ≈ 2 at t = 50. The tumbling motion is also responsible for the slowing down of the two bubbles, as it makes the flow around them fully three-dimensional, which enhances the dissipation in the liquid, hence the drag on each bubble. Tumbling goes on more slowly until t ≈ 85, when the side-byside configuration (θ = 90 • ) and the equilibrium horizontal separation (S re ≈ 2.55) are reached. Ultimately, θ slightly exceeds 90 • , in line with the tiny difference between the final positions of the two bubbles noticed in figure 8(a).

Asymmetric Side Escape

Computations carried out with somewhat larger values of the Galilei number (15 Ga 30) but still fairly low Bond numbers reveal a drastically different evolution of the bubble pair. In what follows, we select the case (Ga, Bo) = (30, 0.3), i.e. M o = 3.3 × 10 -8 , as typical of this regime to discuss the corresponding dynamics. With these parameters, an isolated bubble follows a rectilinear path (Cano-Lozano et al. 2016a,b) and its characteristic rise Reynolds number and aspect ratio at steady state are Re = 66 According to panel (a), the TB also departs from its original vertical path at t ≈ 11, but the LB is left virtually unaffected by this departure until t ≈ 15. At this point, the TB still stands in the wake of the LB but its equatorial plane has tilted clockwise, which distorts the flow in the gap and makes it significantly asymmetric at the back of the LB. This asymmetry induces a slight anticlockwise tilt of the LB path until t ≈ 20, whereas the TB goes on drifting laterally and escapes completely from the wake of the LB. Then the TB rotates anticlockwise in such a way that its equatorial plane becomes again horizontal, and its drift stops at t ≈ 23, when the lateral separation between the two bubbles is approximately S r ≈ 4. Due to this temporary anticlockwise rotation, the TB slightly drifts back until t ≈ 28. The motion of the bubble pair has remained planar up to this point. However, at t ≈ 28, they both start to drift slightly out of their previous plane of rise. Then, each bubble goes on rising along a straight but slightly inclined path.

Both paths being tilted in the same direction but the angle being larger in the case of the LB, the lateral separation weakly but consistently increases over time. As the TB spent part of its potential energy to move out from the wake of the LB during the time period 11 t 23, it never catches up, and the LB remains significantly ahead of the TB in the final configuration.

Figure 11 describes the evolution of the same four characteristics as in figure 9 for the above bubble pair. According to the axisymmetric prediction reported in panel (c), the vertical separation decreases monotonically and the two bubbles enter the initial stage of coalescence at t ≈ 23. However, the three-dimensional evolution reveals a totally different evolution beyond t ≈ 11, although the rise speed of both bubbles and their vertical separation do not exhibit discernible differences with the axisymmetric results up to t ≈ 16. At this time, V T B starts to drop sharply and equals V LB at t ≈ 18. In the meantime, the horizontal velocity of the TB has grown tremendously. Its maximum value, about 40% of V T B , is reached at t ≈ 19. The horizontal velocity of the LB reaches its maximum at the same time. However the ratio of the two maxima is less than 4%, which confirms that the LB is only marginally disturbed by the lateral drift of the TB. This situation dramatically differs from that depicted in figure 9(b), where the two V r maxima have almost the same magnitude. Moreover, in the present ASE scenario, the maximum of V r for the TB is typically three times larger than the maxima encountered during the DKT-type interaction. As figure 11(c) indicates, the horizontal separation has already grown up to S r ≈ 2 by the time V r reaches its maximum (i.e. the TB has left the wake of the LB), and doubles during the next four time units at the end of which the horizontal velocity of the TB vanishes (t ≈ 23). From t ≈ 18 to t ≈ 23, V T B has dropped below V LB , which makes the vertical separation re-increase, from S ≈ 4.5 at t = 17.5 to S ≈ 5.3 at t = 23. The 25% drop of V T B from t ≈ 16 to t ≈ 19.5 emphasizes the fact that a substantial fraction of the kinetic energy of the liquid displaced by the TB has been spent in the meantime to move it laterally and balance the rate of work of the corresponding sideways force. In the next stage (23 t 28), the horizontal velocity of the TB changes sign and reaches a minimum of about -0.06V T B . This negative V r results in the reversed lateral drift already discussed in connection with figure 10. The horizontal velocity still describes some damped oscillations before the vertical and horizontal separations stabilize and reach values close to 5.6 and 4, respectively. Meanwhile, the inclination angle of the bubble pair stabilizes at θ ≈ 36 • . As already pointed out, this configuration is not entirely steady as the slight nonzero slopes of the curves in the right part of figure 11(c) confirm.

In other words, the vertical and transverse components of the interaction force driving the relative position of the two bubbles have not completely vanished yet. However, the remaining values of these components are very small, so that it takes an extremely long time for the system to reach a true steady state. This slow final evolution may be t=10 t=15 t=20 t=24.4

Figure 12. Pre-coalescence dynamics of a bubble pair with Ga = 20, Bo = 0.5; the streamlines in the cross-sectional plane are defined in the reference frame of the leading bubble.

connected to the findings of [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF] for spherical bubbles. Their figure 9(c) indicates that, provided Re > 25, the two bubbles repel each other whatever their separation distance when the inclination of their line of centres is less than a critical angle θ c ≈ 53 • . However, for 30 • < θ < θ c , this repulsive separation-dependent effect is very weak for S 5.0, which corresponds to the present situation. That the equilibrium angle is close to 37 • for Ga = 30, Bo = 0.3 instead of the above value for spherical bubbles is likely an effect of the significant oblateness of the bubbles considered here.

Head-on collision and coalescence

As figures 6 and 7 revealed, increasing the Bond number beyond the Ga-dependent threshold Bo c (Ga) leads to the numerical coalescence of the two bubbles. Most of the time, this coalescence is initiated by a head-on approach, but in some cases the bubbles may also approach each other in an asymmetric manner; e.g. the case (Ga, Bo) = (30, 0.5) in figure 7, which corresponds to near-threshold conditions. Actually, the approach configuration depends on the time elapsed since the bubble pair was released. As discussed in § 4, the attractive effect of the LB wake increases with the bubble deformation, owing to the direct relation between the interface curvature and the strength of the surface vorticity. This is why, for a given Ga, the time at which the two bubbles collide decreases as Bo increases. This leaves less time for non-axisymmetric disturbances to grow before the collision if Bo is significantly larger than Bo c (Ga), favouring the headon configuration. In the case of strongly deformed bubbles, the lift reversal mechanisms discussed in appendix C also act to increase the stability of the in-line configuration.

Here we select conditions (Ga, Bo) = (20, 0.5), i.e. M o = 7.8 × 10 -7 , as an archetype of the head-on coalescence scenario observed in the moderate-Reynolds-number range (Re ≈ 35). As expected from the previous discussion, figure 12 shows that the two bubbles rise almost in a straight line before they numerically coalesce (a tiny lateral deviation of the LB actually takes place in the late stage of the approach, see figure 13(b)). As far as the bubbles evolve independently, their deformation is quite large, characterized by aspect ratios χ ≈ 1.55 and χ ≈ 1.48 for the LB and TB, respectively. In both cases, the aspect ratio is smaller than the critical value χ c = 1.65 beyond which a standing eddy develops at the back of an isolated bubble [START_REF] Blanco | The structure of the axisymmetric high-reynolds number flow around an ellipsoidal bubble of fixed shape[END_REF]. This is why no such structure is present in figure 12(a). Beyond t ≈ 15, the shape of the two bubbles changes significantly. On the one hand, the front part of the LB flattens, due to the proximity of the TB which makes its rise speed increase (see figure 13(a)). On the other hand, the front part of the TB becomes more rounded, owing to the suction induced by the wake of the LB. Finally, the TB catches up with the LB, leaving only a thin liquid film in the gap and making the two bubbles behave essentially as a bluff compound body.

After coalescence takes place (see below), figure 13(a) indicates that the resulting bubble first undergoes a series of large-amplitude oscillations until it relaxes to an oblate shape with an aspect ratio close to 2.0 and a significant fore-aft asymmetry. Due to volume conservation, the radius of this bubble is 2 1/3 ≈ 1.26 times that of the initial bubbles, so that its characteristic parameters are Ga = 2 1/2 × 20 ≈ 28.3 and Bo = 2 2/3 × 0.5 ≈ 0.8, respectively. Under such conditions, the computations of Cano-Lozano et al. (2016a) predict that an isolated bubble still rises vertically, although it is close to the transition to a non-straight path. Present observations are in line with these earlier conclusions, since figure 13(a) indicates that the small horizontal velocity component present at the time of coalescence subsequently decreases over time.

Returning to the near-coalescence stage, figure 14 shows how the interface topology evolves in the stages that just precede and follow numerical coalescence. In its late stage, the film is almost flat. No dimple has formed at its periphery, unlike what is customarily observed with coalescing drops [START_REF] Hartland | The approach of a rigid sphere to a deformable liquid/liquid interface[END_REF][START_REF] Hartland | The profile of the draining film between a rigid sphere and a deformable fluid-liquid interface[END_REF][START_REF] Jones | The film drainage problem in droplet coalescence[END_REF]. This is because, for low-to-moderate Bo and fully mobile interfaces, a dimple forms only when the film has thinned down by several orders of magnitude [START_REF] Yiantsios | On the buoyancy-driven motion of a drop towards a rigid surface or a deformable interface[END_REF]). Indeed, with the same grid resolution, we observed a clear dimple for Bond numbers 1. Numerical coalescence takes place at t ≈ 24.44. Considering that the film starts to form when the gap is 0.5R-thick on the symmetry axis, one can estimate that the computation tracks the drainage process during approximately 1.15 time units. The arguments discussed in appendix B may then be used to estimate the time by which the film would actually rupture under real conditions. Figure 13(a) indicates that the dimensionless approach velocity V a = V T B -V LB of the two bubbles is approximately 0.5 during the early stage of the drainage. From this, the approach Weber number, W e a = V 2 a Bo, and the approach capillary number, Ca a = V a Bo/Ga = µV a /γ, are found to be 0.125 and 0.0125, respectively. Then, for nearly-spherical bubbles, the estimate (B 1) predicts a dimensionless inertial drainage time T di ≈ 0.27, while (B 2) (once properly transposed to the case of two identical bubbles) predicts a viscous drainage time T dv ≈ 2.5. Hence, the drainage is controlled by viscous effects and the limited resolution would shorten it by ≈ 1.35 time units, would the initial bubble oblateness be small. Given that the tandem rises with the average velocity V m = 1 2 (V T B + V LB ) ≈ 2.4, this implies that the vertical position at which coalescence actually happens would be underestimated by 3.25 bubble radii in this limit. The non-negligible bubble oblateness certainly increases the actual drainage time. For instance, the actual value of the parameter k i in (B 1) is 2.05 for χ = 1.5 [START_REF] Duineveld | Bouncing and coalescence of two bubbles in water[END_REF], so that the correct estimate for the inertial drainage time is T di ≈ 0.5. The quantitative influence of the bubble oblateness on the viscous drainage time is unknown but presumably increases also significantly T dv .

As the last two snapshots in figure 14(a) show, the neck radius of the resulting bubble grows very fast after coalescence happens. More specifically, the growth law (not shown) is found to be (t-t coal ) 0.36 over the very first stage following the coalescence time instant t coal . This is somewhat slower than the (t -t coal ) 1/2 self-similar behaviour observed in the experiments of [START_REF] Paulsen | Coalescence of bubbles and drops in an outer fluid[END_REF] and confirmed theoretically and numerically by [START_REF] Munro | Thin-sheet flow between coalescing bubbles[END_REF] and [START_REF] Anthony | Scaling laws and dynamics of bubble coalescence[END_REF], respectively. However [START_REF] Paulsen | Coalescence of bubbles and drops in an outer fluid[END_REF] noticed that the growth rate reduces gradually once the neck radius is larger than 0.3R. Here, due to the limitations inherent to minimum cell size δ min = R/272 and to the numerical procedure used to let the initial two interfaces merge, we barely observe the neck before it reaches this radius, which is certainly the reason for the above slower growth rate. M o 8 × 10 -7 . The DKT scenario is observed for the lower three values of Bo, while the pair with Bo = 0.2 eventually experiences head-on coalescence. For lower Bo, the interaction always leads to the side-by-side configuration and the lateral separation stabilizes at a value close to 2.5. While the Bond number does not have any noticeable influence on the final path in this regime, it affects the critical time by which the tumbling process sets in. Owing to the connection between the bubble oblateness and the amount of vorticity produced at its surface, the more oblate the LB is the stronger the attractive wake effect is. Therefore, at a given time, the larger Bo the shorter the separation between the two bubbles. From figures 15(b) -(d), it may be inferred that the axial symmetry of the flow breaks down when S = S c ≈ 2.6, a critical value reached in a shorter time as Bo increases.

Influence of deformation and initial conditions

Figure 16 shows the same path characteristics for Ga = 30 and Bond numbers ranging from 0.02 to 0.5, i.e. 9.9 × 10 -12 M o 1.5 × 10 -7 . Here, all cases with Bo 0.45 correspond to the ASE scenario, whereas the two bubbles coalesce in an asymmetric way for Bo = 0.5 (see figure 7). In all non-coalescing cases, it is found that the larger Bo the smaller the long-term lateral deviation of the TB (figure 16(a)). Hence, the long-term lateral separation decreases with Bo (figure 16(b)), since the LB hardly deviates from its initial vertical path except for Bo = 0.45. The vertical separation S is also seen to decrease with the Bond number. This is merely a geometrical consequence of the decrease of S r . Indeed, the shorter the time spent by the TB in its lateral motion, the shorter the slowing down of its vertical motion, hence the smaller the increase of the vertical separation during the lateral escape stage.

The long-term inclination of the bubble pair exhibits more complex variations. First, it increases with the Bond number up to Bo = 0.2. Then it slightly decreases for Bo = 0.3 until it experiences a sharp increase for Bo = 0.45. In the whole range Bo 0.3, the long-term inclination angle stands in the range 30 • < θ < 40 • , in agreement with the experimental findings of [START_REF] Kusuno | Experimental investigation of the motion of a pair of bubbles at intermediate Reynolds numbers[END_REF] obtained under similar conditions (Re < 150). The situation corresponding to Bo = 0.45 is specific. Indeed, the two bubbles are close to coalescing asymmetrically at some point (S ≈ 2 at t = 22.5). This is why the lateral motion of the TB significantly disturbs the subsequent motion of the LB which is seen to experience a series of damped oscillations before rising along a rectilinear, slightly inclined path. Moreover, since the vertical separation of the two bubbles is small just before the TB starts to drift laterally (S ≈ 3.8 at t = 18), the inclination of tandem after this drift has been completed is significantly larger (θ ≈ 60 • ) than those found with smaller Bo. We shall come back to the long-term evolution of this bubble pair later. Beyond these geometric features, the main information provided by figure 16 is that the smaller Bo is the earlier the side escape of the TB starts. This is in stark contrast with the observations made on figure 15 for the DKT regime, where bubble deformation is found to promote the tumbling process. Moreover, considering the critical time at which S r and θ depart from zero for the various Bo reveals that the corresponding vertical separation does not keep a constant value (figures 16(b) -(d)). Rather it decreases from S c ≈ 8 for Bo = 0.02 to S c ≈ 5.8 for Bo = 0.3. These values are significantly larger than the equilibrium separation predicted in the axisymmetric configuration which, according to figure 5, ranges from S e ≈ 6.1 for Bo = 0.02 to S e ≈ 4.3 for Bo = 0.3. Moreover, in all cases, this critical separation is much larger than that at which tumbling is found to start in the DKT regime. This finding indicates that, unlike the breakdown of the axial symmetry in the DKT scenario, the ASE mechanism is driven by a long-range interaction. This is in line with the conclusion of [START_REF] Yin | Lattice-Boltzmann simulation of finite Reynolds number buoyancydriven bubbly flows in periodic and wall-bounded domains[END_REF] who simulated the motion of buoyancy-driven suspensions of spherical non-coalescing bubbles at Re ≈ 10.

To get additional insight into the long-term behaviour of the bubble pair, figure 17 shows the evolution of the separation and inclination angle over longer times for the two sets (Ga, Bo) = (30, 0.3) and (Ga, Bo) = (30, 0.45) , i.e. M o = 3.3 × 10 -8 and M o = 1.1×10 -7 , respectively. In the former case, the two components of the separation are seen to slightly increase until the end of the computation but the inclination angle changes by less than 1 • over the last forty time units, reaching the 'asymptotic' value θ ≈ 36 • . Conversely, in the latter case, the two components of the separation exhibit significant variations throughout the computation. The vertical separation gradually tends to zero while S r goes on increasing even for t ≈ 100. At this final time, θ is close to 83 • and is still gently increasing, so that there is little doubt that the system eventually reaches a perfect side-by-side configuration. To understand why the two bubble pairs behave so differently on the long term, it must first be noticed that, at the end of the lateral drift of the TB (t ≈ 40), the distance S = (S 2 + S 2 r ) 1/2 between the two centroids is approximately 6.8 for Bo = 0.3 and 4.5 for Bo = 0.45, the Reynolds numbers being approximately 66 and 60, respectively. For spherical bubbles, the results of [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF] indicate that, for Re = 100, the torque acting on the tandem is negligibly small for θ 25 • when S 6 but keeps significant values whatever θ for S = 4.5. Although bubble deformation is large in the two sets considered in figure 17, this is a strong indication that, after the TB has drifted laterally, only the pair with Bo = 0.45 still experiences a noticeable torque. Consequently only this pair may reach the side-by-side configuration in a reasonable time. According to figure 16, S decreases with Bo. Hence, it may be concluded that in the ASE scenario, the distance between the two centroids after the TB has completed its drift is usually too large for the torque to remain significant. Therefore the inclination of the tandem remains almost unchanged in the subsequent stages. Only bubbles whose deformation is close to the coalescence threshold (here Bo ≈ 0.5) escape this rule, since the vertical separation of the corresponding pairs is very small when the side escape of the TB takes place.

Interestingly, for Bo = 0.45, S r goes on increasing at the end of the computation, suggesting that the interaction force is still nonzero and repels the two bubbles. At the same Re and S , the computational results of [START_REF] Legendre | Hydrodynamic interactions between two spherical bubbles rising side by side in a viscous liquid[END_REF] indicate that the interaction force between two spherical bubbles rising side-by-side is attractive. The difference between the two situations may be understood by noting once again that wake effects are much stronger for Bo = 0.45 than for Bo = 0. Since these effects are repulsive in the side-by-side configuration, the critical Reynolds number at which the interaction force switches from repulsive to attractive is significantly larger in the present case, which explains the observed behaviour.

Role of the TB shape in near-critical conditions for coalescence

The global geometric indicators reported in figure 16 reveal how the two bubbles get closer to coalescence as the Bond number increases. Nevertheless it is of interest to examine also how the bubble shapes vary with Bo close to the corresponding threshold, especially with respect to the potential influence of the bubble distortion on the magnitude and even the direction of the sideways force acting on the TB through the two mechanisms reviewed in appendix C. Figure 18 displays the evolution of the shapes and relative positions of the two bubbles for near-critical conditions at Ga = 30 and 15, respectively. It is striking that, for both values, the TB shape exhibits almost no left/right asymmetry until the late stage of the interaction. This is a strong indication that the A-mechanism discussed in appendix C cannot reduce significantly the sideways force. Similarly, the flow characteristics are such that the S-mechanism cannot take place. For instance, with Ga = 30 and Bo = 0.4, the aspect ratio of the TB at t = 14 is approximately 1.6 and its Reynolds number is close to 80. For these parameters, the path of an isolated bubble rising in a fluid at rest is stable and the numerical results of [START_REF] Adoua | Reversal of the lift force on an oblate bubble in a weakly viscous linear shear flow[END_REF] (their figure 5) indicate that finite-Reynolds-number effects decrease the lift force by only 15% with respect to the inviscid prediction. So, it can be concluded that within the parameter range considered in the present study, the lateral migration of the TB is merely driven by the standard shear-induced mechanism. Hence, close to the threshold, what makes the difference between non-coalescing and coalescing situations is essentially the enhancement of the attractive wake effect as the TB becomes increasingly oblate. For a sufficiently large Bo, this attractive effect becomes so strong that it delays the occurrence of the lateral instability of the TB to such an extent that, although this bubble subsequently migrates 'normally', its migration lasts for a too short time to prevent it from hitting the LB.

Influence of an initial angular deviation

Despite sophisticated bubble release systems, tiny initial lateral deviations can hardly be avoided in laboratory experiments, yielding small-but-nonzero angular deviations of the tandem with respect to the perfect in-line configuration. This calls for the investigation of the influence of an initial nonzero inclination on the subsequent dynamics. For this purpose, we systematically examined the impact of a nonzero initial inclination 0 • < θ 0 2 • throughout the (Ga, Re)-range considered in § § 5 and 6. Note that in all configurations considered in this subsection, the initial horizontal separation between the bubble centres is several times larger than the finest grid cells located on both sides of the bubbles surface (4 times larger for the smallest inclination). This initial configuration results in a 'macroscopic' asymmetry of the discretized solution, in which the tiny numerical asymmetries of the pressure field returned by the Poisson solver (see § 6.2 and appendix A) play no role.

For Ga = 10, we observed that under such conditions, the DKT scenario no longer takes place for Bo 0.1. Instead, the system follows a clear ASE evolution. For instance, selecting (Ga, Bo) = (10, 0.1), the TB starts drifting laterally when S ≈ 4 (resp. S ≈ 5) if θ 0 is set to 1 • (resp. 2 • ), leaving the path of the LB almost unaffected in both cases. Instead of ending up in the side-by-side configuration as it does when θ 0 = 0 • , the tandem then reaches an approximate final inclination of 50 • (resp. 40 • ). That θ 0 has such a dramatic influence on the existence of the DKT regime and the late geometry of the tandem is at variance with observations reported for rigid bodies falling in tandem. Indeed, experiments performed with short cylinders [START_REF] Brosse | Interaction of two axisymmetric bodies falling in tandem at moderate Reynolds numbers[END_REF] indicate that the DKT regime is still observed when a small initial lateral offset is imposed to the trailing body. The key difference with a bubble pair is the relative magnitude of attractive effects which is much stronger for rigid bodies, owing to the different vorticity generation mode at the body surface. Moreover, for a given body shape, shear-induced lift effects are significantly weaker for rigid bodies with Reynolds numbers of O(10) or larger, owing to the presence of large separated regions [START_REF] Kurose | Drag and lift forces on a rotating sphere in a linear shear flow[END_REF]. These two factors make the DKT configuration much more sensitive to small lateral deviations in the case of nearly-spherical bubbles.

Figure 19 provides the evolution of the tandem geometry in the case (Ga, Bo) = (30, 0.3) for initial inclinations θ 0 ranging from 0 • to 2 • . The system is found to follow the ASE scenario in all cases. However its final geometry strongly depends on θ 0 . In particular, as figure 19(d) shows, the final inclination of the line of centres decreases from 36 • for θ 0 = 0 • to 22 • for θ 0 = 2 • . This is because when θ 0 is nonzero, the TB starts drifting laterally soon after it is released from rest, which is not the case for θ 0 = 0 • . Indeed, with θ 0 = 0 • , the initial flow configuration is no longer axisymmetric and the TB faces an asymmetric wake as soon as vorticity generated at the LB surface has been advected downstream over an O(S)-distance. Consequently, the TB experiences a nonzero sideways force much earlier than in the perfect in-line configuration, where this force occurs only after the system has become unstable. A nonzero θ 0 shortens the initial time period after which S r starts to grow, as is made clear in figure 19(c). A direct consequence of this shortening is the fact that the vertical separation at which the lateral drift starts is larger for θ 0 = 0 • (S ≈ 8 for θ 0 = 2 • instead of S ≈ 6.8 for θ 0 = 0 • ). For this reason, weaker velocity gradients across the LB wake subsist for θ 0 = 0 • at the position of the TB when it starts drifting, resulting in a weaker sideways force. This makes the maximum of V r smaller (figure 19(b)), yielding a shorter final lateral position, hence a smaller final inclination of the tandem.

Figure 20 displays the influence of the Bond number on the evolution of bubble pairs with Ga = 30 when θ 0 is set to 2 • . This figure is the counterpart of figure 16 discussed above for θ 0 = 0 • . All pairs with Bo < 0.7, i.e. M o < 4.2×10 -7 , are seen to follow a clear ASE scenario. For the aforementioned reason, the final inclinations reported in figure 20 are significantly smaller than those observed in figure 16. In addition, they exhibit a marked and consistent decrease with the Bond number, a trend which is absent from figure 16. For pairs with Bo 0.8, θ is found to decrease over time until the two bubbles perfectly align vertically, which eventually forces them to coalesce. The near-critical case Bo = 0.7 is quite specific, and in many instances similar to the situation encountered for Bo = 0.45 in figure 16. In this case, the two bubbles are very close to coalescing at Y ≈ 58 in figure 21(a). The corresponding gap is so thin that the flow at the back of the LB is strongly disturbed, forcing the latter to deviate abruptly from its vertical path. This allows the lateral separation to increase beyond the critical value S rc 2 within a short time lapse (not visible in figure 21(b) which is limited to shorter times), allowing the tandem to avoid coalescence. Similar to the evolution displayed in figure 17(b), but through a sharper transition, the inclination of the tandem grows until the side-by-side configuration is reached. Given the small gap at which the transition takes place and the quite symmetric final lateral positions of the two bubbles, this specific evolution is closer to the DKT scenario than to a standard ASE evolution. The critical Bond number beyond which the two bubbles coalesce stands in the range 0.45 -0.5 in figure 16 (θ 0 = 0 • ), but is slightly larger than 0.7 for θ 0 = 2 • . In line with the discussion on figure 19, the reason for this marked increase stems directly from the longer vertical distance over which the shear-induced lift force acts on the TB when θ 0 = 0 • , and the slightly shorter lateral distance this bubble has to drift to avoid hitting the LB. These two factors imply that, compared to the reference case, a weaker positive lift force is sufficient to avoid coalescence when θ 0 is nonzero. The two lift reversal mechanisms discussed in appendix C being directly linked to the ability of the TB to deform, increasing the Bond number makes the lift force decrease and eventually change sign, other things being equal. Hence, in the presence of an initial angular deviation, coalescence can be avoided over a broader range of distortion of the TB, i.e. up to a larger Bond number.

The fate of all bubble pairs released with an initial inclination θ 0 = 2 • over the range 10 Ga 30, 0 < Bo 2.0 is summarized in figure 21. This figure is the counterpart of figure 6 obtained with θ 0 = 0 • . In line with the above discussion, the ASE regime observed when θ 0 = 2 • exists over a significantly broader range of Bo, hence for bubbles with a larger oblateness. For instance, bubbles with aspect ratios χ 1.3 (resp. 1.7) are found to coalesce for Ga = 15 (resp. 30) in the absence of any initial deviation. With θ 0 = 2 • , the corresponding critical aspect ratios raise beyond 1.5 (resp. 2.0) and the critical Reynolds numbers are close to 50 and 120, respectively.

Figure 22 shows how the bubble shapes and relative positions evolve close to the transition to coalescence for Ga = 30. Unlike those reported in figure 18, the TB shapes now exhibit a marked left/right asymmetry (although less pronounced, this trend subsists for lower Ga). The difference with the strict in-line configuration is that here the TB is fully immersed in an asymmetric flow as soon as the wake of the LB has developed within the gap. Consequently, its shape has to adapt to this asymmetric environment throughout the interaction process. While the TB is seen to migrate toward the right for Bo = 0.7, a tiny migration toward the left may be identified for Bo = 0.8, leading unavoidably to coalescence. The discussion in appendix C allows the origin of this reverse migration to be readily identified. The aspect ratio of the TB being only 1.85 for Bo = 0.8, the S-mechanism is not present, since it only takes place beyond a threshold χ cS ≈ 2.2 ( we observed a reversed migration due to this mechanism by increasing the Bond number beyond unity). In contrast, the egg-like shapes of the TB point to the A-mechanism.

The capillary number based on the rise velocity of this bubble at t = 14 only increases from 0.051 for Bo = 0.7 to 0.058 for Bo = 0.8. However this modest increase turns out to be sufficient for the deformation-induced force directed to the left to exceed the shear-induced lift directed to the right, preventing the lateral escape of the TB.

Influence of initial separation

All simulations considered up to now are based on an initial separation distance S 0 = 8. It is relevant to examine how this choice influences the upcoming evolution of the bubble pair. Since the in-line configuration is stable in the head-on coalescence regime, S 0 is not expected to have an influence in this regime. This is why we consider the influence of S 0 only in the DKT and ASE regimes.

Figure 23 displays the evolution of some characteristics of the bubble pair obtained by increasing S 0 from 8 to 12 in the case Ga = 10, Bo = 0.02. With S 0 = 8, this set of parameters yields the DKT configuration described in figures 8 and 9. Introducing an appropriate S 0 -dependent time shift t 0 and setting t * = t -t 0 (S 0 ) allows the t *evolutions of all quantities for the three initial separations to collapse on a single curve beyond t * ≈ 10. Hence a DKT scenario yielding the same final configuration is observed whatever the initial separation. This is no real surprise since, for this set of parameters, the in-line configuration becomes unstable only when the two bubbles get very close, which happens only for t * ≈ 40.

Figure 24 shows how the evolution of the same parameters vary with S 0 for Ga = 30, Bo = 0.3. In this case, the interaction process yields an ASE scenario whatever S 0 . However the final configuration now depends on the initial separation. In particular, the shorter S 0 the larger the final inclination of the tandem, with θ ≈ 47 • and θ ≈ 22 • for S 0 = 6 and S 0 = 12, respectively. Not unlikely, the larger S 0 is the longer it takes for the path of the TB to start bending. Examining the evolutions of the vertical separation and those of the rise speed of each bubble reveals that none of these three quantities exhibits a constant value by the time the TB starts drifting laterally. In contrast, it turns out that the difference V T B -V LB is close to 0.35 whatever S 0 when this lateral motion sets in. Therefore we conclude that the in-line configuration becomes unstable when the velocity difference between the two bubbles exceeds the above threshold, irrespective of the separation distance at which this threshold is reached.

Closely related to the influence of the initial separation is that of the sequential release of the two bubbles in actual laboratory experiments. Still with Ga = 30, Bo = 0.3, we ran a computation in which the TB was allowed to start rising only after the LB reached a dimensionless height S(t 0 ) = 8 above the point of release. Figure 25(a) displays the corresponding evolution of the vertical separation. Due to the time lapse the TB needs to reach its final rise speed, S reaches a maximum close to 10.2 before the attractive interaction sets in and the separation starts to decrease. This is why we compare the subsequent evolution with that obtained for S 0 = 10 when the two bubbles released simultaneously. For this purpose, we introduce an appropriate time shift t 0 and define a modified time t * = t + t 0 to make the two evolutions coincide in the early stage of the interaction. With this time shift, the two further evolutions of S almost superimpose, as do those of the lateral velocity of the TB (figure 25(b)) and the inclination angle (figure

25(c))

. Consequently it can be concluded that the time elapsed in between the release of two successive bubbles merely increases their 'initial' separation, defined as the distance that separates them before the wake of the LB starts to influence the rise of the TB.

Summary and concluding remarks

We carried out a series of three-dimensional simulations in order to dissect the physical mechanisms involved in the evolution of a pair of clean, deformable rising bubbles initially released in line. In this first part of the investigation, we focused on the parameter range 10 Ga 30, 0.02 Bo 2.0 which corresponds to inertia-dominated regimes in which the path of an isolated bubble remains straight and vertical. We made use of the Basilisk open source code and improved on the original version by implementing a specific AMR improvement [START_REF] Zhang | Vortex interactions between a pair of bubbles rising side by side in ordinary viscous liquids[END_REF]) allowing an automatic grid refinement in the gap left between the two bubbles when they come almost in contact. Previous theoretical and computational investigations of this configuration essentially considered spherical bubbles and steady or quasi-steady configurations. By allowing the bubbles to move deform freely, the present study provides a more realistic description of the hydrodynamic interactions governing the fate of the bubble pair.

To build on a reference case, we first ran axisymmetric time-dependent simulations. Similar to the predictions of [START_REF] Yuan | On the in-line motion of two spherical bubbles in a viscous fluid[END_REF] and [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF] for spherical bubbles, we found that the two bubbles stabilize a finite distance apart, provided their deformation remains moderate. Since the strength of the inertia-induced repulsive interaction force increases with Ga, so does the critical Bond number beyond which the bubbles come in contact and eventually coalesce. For this reason, the equilibrium distance depends on both Ga and Bo, or equivalently, on the Reynolds and Weber numbers. The empirical correlation (5.2) summarizes the corresponding findings.

We then turned to fully three-dimensional evolutions. The simulations revealed that, under conditions for which an equilibrium separation distance exists in the axisymmetric case, this configuration is never reached in the three-dimensional case because the in-line arrangement is unstable with respect to non-axisymmetric disturbances. As rationalized long ago [START_REF] Harper | On bubbles rising in line at large Reynolds numbers[END_REF], this is because, provided bubble deformation is moderate, any deviation of the TB from the axis of the LB wake is amplified, owing to the shearinduced lift force that tends to drive it out of the wake. This lateral drift takes two markedly different forms, depending on the strength of inertial and deformation effects. For Ga = O(10) and Bo 0.1, a DKT mechanism takes place. In this regime, the two bubbles deviate almost symmetrically from their initial paths and this deviation happens while they get very close to each other, the remaining gap being of the order of the bubble radius or even less. After the tumbling stage is completed, the two bubbles rise virtually side by side in a vertical plane, their horizontal centre-to-centre separation being approximately 3 bubble radii. In more inertial regimes, say 12 Ga 30, and low-to-moderate Bond numbers, the bubble pair evolves in such a way that the TB escapes laterally from the wake without significantly altering the path of the LB. In this Asymmetric Side Escape scenario, the lateral drift of the TB takes place while the two bubbles are still widely separated. In the subsequent stage, the tandem maintains a significant and almost constant inclination in the range 30 • θ 40 • , except for near-coalescence conditions under which the side-by-side configuration may eventually be reached. Whatever Ga, coalescence is avoided only up to a critical Bond number, the value of which increases with Ga, from Bo c ≈ 0.1 for Ga = 10 to Bo c ≈ 0.5 for Ga = 30. Actually, bubble deformation influences the evolution of the system in a number of ways. In particular, the fact that the vorticity generated on a bubble is directly proportional to the curvature of its surface implies that the strength of the wake-induced interaction acting on the TB increases sharply with the Bond number. This interaction being attractive, deformation effects are found to promote the suction of the TB along the centreline of the LB wake. This of course favours coalescence, but also tends to stabilize the in-line configuration when Bo < Bo c (Ga) by delaying the growth of nonaxisymmetric disturbances.

We performed a detailed examination of the influence of initial conditions, especially of a possible misalignment of the two bubbles. Even a minimal initial deviation (θ 0 2 • ) was found to have a dramatic influence on the evolution of the bubble pair. In particular, the DKT regime previously observed for Ga = O(10) and Bo 0.1 no longer takes place. Instead the system follows an ASE evolution. A nonzero initial inclination also promotes the ASE configuration toward larger bubble deformations, making the critical Bond number increase up to 0.8. The reason is that for θ 0 = 0 • the TB faces an asymmetric flow from the very beginning of its rise, which makes it able to drift laterally over a longer time than in the canonical θ 0 = 0 • case. This scenario is efficient to avoid coalescence as far as the sideways force is dominated by the classical shear-induced lift mechanism. However, mechanisms leading to a weakening or even a reversal of the overall transverse force exist for non-spherical bubbles. In particular, bubbles rising in a shear flow exhibit a non-axisymmetric shape, a feature known to produce a deformation-induced component of the transverse force with opposite sign compared to that of the shear-induced inertial lift. As the Bond number approaches its critical value, we observed that the asymmetry of the carrying flow in which the TB is immersed for θ 0 = 0 • yields a pronounced egg-like shape of the latter. This goes hand in hand with a reduction of the transverse force which eventually changes sign for Bo > Bo c , forcing the two bubbles to coalesce.

The present study clarifies the respective roles of inertial, viscous and capillary effects, as well as that of initial conditions, in the three-dimensional dynamics of the considered system. From the standpoint of the microstructure of bubbly suspensions, the main outcome is presumably the final geometry of the arrangement that emerges from the evolution of bubble pairs initially released in line, possibly with some small angular deviation. As we saw, only pairs made of nearly-spherical bubbles with O(10)-Galilei numbers end up in a side-by-side configuration and thereby tend to favour the formation of horizontal clusters. However, in most non-coalescing situations, the interaction process follows the ASE scenario, yielding final inclinations in the range 15 • -40 • and horizontal separations ranging from 2 to 5 radii, depending on Ga, Bo and θ 0 . Given the unavoidable variations of initial bubble positions in real bubbly flows, this significant range of nearequilibrium inclinations and separations guarantees a much more homogeneous spatial bubble distribution on the long term than what can be expected from simplified models assuming potential flow and/or spherical bubble shapes.

This study did not examine the influence of slight differences in the size of the two bubbles. Nevertheless the main consequences of such a difference may be inferred from present findings. Suppose that the LB is slightly larger than the TB. In the (Ga, Bo)range considered here, the corresponding increase in the buoyancy force makes the rise speed of the LB increase, lowering the positive velocity difference V T B -V LB during the axisymmetric stages of the interaction process. Therefore, compared to the reference case, the two bubbles maintain a larger separation at a given time, which favours the occurrence of an ASE-type scenario. Consequently, this size difference tends to broaden the subdomain corresponding to the ASE regime in the phase diagram of figure 6. Conversely, if the TB is slightly larger than the LB, the two bubbles are more prone to getting close to each other before the axial symmetry of the flow breaks down. Hence, if Ga and Bo are such that the system stands close to the DKT-ASE transition, this configuration favours the DKT scenario. Similarly, it lowers the critical Bond number Bo c (Ga) if the system is close to the coalescence threshold. The above predictions may be corroborated with the experimental observations of [START_REF] Kusuno | Lift force acting on a pair of clean bubbles rising in-line[END_REF] who considered two bubble pairs close to the DKT-ASE transition, with a LB corresponding to Ga = 13.5, Bo = 0.27 in both cases. With a TB 3.5% larger than the LB, they found the system to follow a DKT transition (their figures 4(c) and 5(c)). Conversely, they observed a clear ASE scenario when the TB was 3.5% smaller than the LB (figures 4(d) and 5(d)).

Another aspect that the present study leaves untouched is the fundamental question of the mathematical nature of the bifurcation that takes place when the axisymmetric configuration becomes unstable. Similarly, although the DKT and ASE regimes exhibit markedly different physical characteristics, whether or not they correspond to truly distinct unstable modes of the system remains unknown at this stage. Tackling these issues requires the development of an appropriate global linear stability approach. Numerical tools allowing the threshold and nature of bifurcations involved in the wake of fixed [START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF]Cano-Lozano et al. 2016b) or freely-moving [START_REF] Tchoufag | Linear instability of the path of a freely rising spheroidal bubble[END_REF] clean isolated bubbles with a prescribed shape have become available during the past decade. More recently, the same approach was extended to fully deformable isolated bubbles [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]. Further extending this approach to systems involving bubble pairs is certainly feasible and seems the natural next step capable of bringing new insight into these fundamental issues. that the initial configuration being axisymmetric, the transition to a non-axisymmetric state is governed by the asymmetry of numerical disturbances. Two possibilities exist to assess the influence of these disturbances. If one looks for a rigorous grid convergence study, imposing a well-defined small asymmetry to the initial conditions and examining how it changes the results as the grid is varied is the appropriate choice. However, a code user has no other control on the disturbances unavoidably present in the discretization procedure and the time-advancement algorithm than changing the user-defined numerical parameters, which in Basilisk are the grid refinement level and the tolerance on the Poisson solver (see below). Therefore, from a user point of view, it is more relevant to examine how the 'natural' disturbances that arise in the numerical solution influence the results, while the above two parameters are varied separately. This is the choice adopted in the tests reported below. The Poisson solver is the only part of the time-advancement algorithm used in Basilisk that does not preserve spatial symmetry up to machine accuracy [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. For this reason, we first examined how much the user-specified tolerance T applied to this solver influences the evolution of the bubble pair. This tolerance is defined as the maximum relative change during one time step of the fluid volume enclosed in a cell, i.e. the maximum of |∇ • u|∆t over the computational domain. The standard tolerance used throughout the computations is T = 1 × 10 -4 . In this test, we ran two additional simulations, with T = 1 × 10 -2 and 1 × 10 -6 , respectively. Figure 26 reveals that the smaller T is, the longer it takes for the axial symmetry of the system to be broken (see the evolution of S r in figure 26(b)). This was to be expected since reducing the tolerance reduces the asymmetry of the solution, delaying the transition to three-dimensionality. This of course has some impact on the final vertical separation of the two bubbles, hence on the inclination of their line of centres. In contrast, the bubble velocities before and after the lateral escape of the TB are left unchanged by the change in T , together with the final horizontal separation.

Then we carried out a grid convergence study, setting T = 1 × 10 -4 in all cases. Starting from the reference grid with ∆ min = R/68 used throughout the study, we ran the same case on a grid twice coarser (∆ min = R/34) and another grid twice finer (∆ min = R/136). Note that, since the time step decreases with ∆ min , the computation with ∆ min = R/136 was extremely time consuming. This is why it was stopped at t = 25, beyond which the dynamics of the flow is not expected to reveal any supplementary influence of the grid resolution. The results of these tests are presented in figure 27. Again, these results show that the only significant change in the velocities evolution is the time by which the TB starts to escape laterally, and consequently the final vertical separation of the two bubbles (although only evolutions for t 25 are displayed in the figure, the horizontal separation was found to relax to the same value with the two grids on which the runs were carried out over larger times). It is noticeable that the onset of three-dimensionality is reached earlier on the reference grid than on both the finer and coarser ones. To understand this surprising feature, it must be kept in mind that the time step limitation arises from capillary effects, implying ∆t ∝ ∆ 3/2 min . Therefore, the maximum error on |∇ • u| resulting from the Poisson solver is proportional to ∆ -3/2 min . Assuming that the cell size at the location where this maximum is reached is ∆, the corresponding error on u is proportional to ∆ -3/2 min ∆. Since the local cell volume is ∆ 3 , the contribution of the error to the local fluid momentum is proportional to ∆ -3/2 min ∆ 4 . Indeed, what determines the influence of the error on the onset of the lateral motion is the asymmetric component of the momentum over a cell rather than that of the local velocity u itself. The position of the maximum error within the computational domain is unknown a priori. It may vary from one grid to another, and this variation is responsible for the non-monotonic behaviour observed in figure 27. For instance, if the maximum is reached in the most refined subregion (i.e. very close to one interface) on the coarsest grid (∆ = ∆ min = R/34) while it is reached in the near wake on the other two grids (∆ = 4∆ min = R/17 and R/34, respectively), the asymmetric contribution to the momentum (normalized with (R/34) 5/2 ) is of O(1) on the coarsest grid, while it is of O(2 11/2 ) and O(2 3 ) on the intermediate and finest grids, respectively. In such a case, the onset of the TB lateral escape is expected to happen first on the intermediate grid, then on the finest one, and finally on the coarsest one, just as observed in the figure.

In summary, the above tests revealed that, as expected, the minimum cell size and tolerance on the Poisson solver influence the time by which the axial symmetry of the initial configuration breaks down. The smaller the tolerance T , the longer it takes for the bifurcation to take place on a given grid. The influence of the grid resolution is more complex, owing to the variability of the position at which the maximum asymmetry on the velocity field takes place. Because of this, the time by which the solution develops a significant three-dimensional component depends on both ∆ min and the local cell size ∆ at the location of this maximum. Obviously, a longer transition time being equivalent to as shorter separation between the two bubbles, the sensitivity of the system to these numerical parameters implies that the value of the critical Bond number Bo c (Ga) has a nonzero numerical 'error bar'. For instance, although the bubble pair with (Ga = 30, Bo = 0.45) is found to escape coalescence with ∆ min = R/64 and T = 1 × 10 -4 , it is very likely that it coalesces with a smaller T . The numerical uncertainty on Bo c (Ga) may be estimated by considering the time lag for the onset of the TB lateral escape resulting from the sensitivity of the system to T and ∆ min . Given the difference in the rise speed of the two bubbles by the time this escape starts, figures 26 and 27 indicate that the vertical separation between the two bubbles is approximately reduced by ∆S = 0.3 when T is reduced by two orders of magnitude or when the grid resolution is increased from ∆ min = R/68 to ∆ min = R/136. Then, figure 16(b) indicates that the minimum vertical separation reached during the lateral escape varies from 4.55 for Bo = 0.3 to 2.0 for Bo = 0.45, from which a rule of three suggests that a reduction of the vertical separation by ∆S = 0.3 implies a reduction of the critical Bond number by 1.8%.

Appendix B. Bubble coalescence in pure liquids

Coalescence of drops and bubbles in a suspending liquid has received considerable attention in the literature owing to its many applications. From the fluid mechanics viewpoint, most of the studies carried out during the second half of the past century attempted to determine the characteristics of the drainage of the film that forms in between two drops or bubbles when they approach each other (see [START_REF] Chesters | The modelling of coalescence processes in fluid-liquid dispersions: a review of current understanding[END_REF] and [START_REF] Chan | Film drainage and coalescence between deformable drops and bubbles[END_REF] for reviews). Assuming an axisymmetric geometry, asymptotic studies based on the lubrication approximation and numerical studies based on the boundary integral method (e.g. [START_REF] Chi | A theoretical study of the motion of a viscous drop toward a fluid interface at low Reynolds number[END_REF]) revealed the critical role of the drop/bubbleto-external fluid viscosity ratio in the drainage dynamics. Depending on how this ratio, λ say, compares with the film aspect ratio (the typical radius-to-thickness ratio of the near-contact region), four possible situations arise, corresponding to nearly immobile (λ 1/2 ), partially mobile (λ ∼ 1/2 ), mobile ( -1/2 λ 1/2 ) and fully mobile (λ -1/2 ) interfaces, respectively [START_REF] Davis | The lubrication force between two viscous drops[END_REF][START_REF] Chesters | The modelling of coalescence processes in fluid-liquid dispersions: a review of current understanding[END_REF]). In the first three cases, provided the drainage takes place under the action of a constant external force and the drops are nearly spherical except in the near-contact region, the minimum film thickness obeys a power law evolution [START_REF] Hartland | The approach of a rigid sphere to a deformable liquid/liquid interface[END_REF][START_REF] Jones | The film drainage problem in droplet coalescence[END_REF][START_REF] Yiantsios | On the buoyancy-driven motion of a drop towards a rigid surface or a deformable interface[END_REF][START_REF] Nemer | Coalescence of drops with mobile interfaces in a quiescent fluid[END_REF]. The corresponding exponent depends on the flow and boundary conditions, but stands in between 0 and -1 in all cases. A crucial consequence of this algebraic thinning rate is that the drainage requires an infinite time to be completed. Hence, non-hydrodynamic effects, in the first place the long-range London-van der Waals force, are required for coalescence to be achieved in a finite amount of time.

Things differ drastically with fully mobile interfaces, which is the relevant situation for bubbles in pure liquids (contamination by surfactants yields immobile or at least partially mobile interfaces; e.g. [START_REF] Vakarelski | Mobile-surface bubbles and droplets coalesce faster but bounce stronger[END_REF][START_REF] Vakarelski | Free-rising bubbles bounce more strongly from mobile than from immobile water-air interfaces[END_REF]). In this case, the flow in the gap is merely a plug flow [START_REF] Davis | The lubrication force between two viscous drops[END_REF], as interfaces offer no resistance to the squeezing of the film. Under such conditions, standard lubrication approximations do not hold unless the gap has become extremely thin [START_REF] Davis | The lubrication force between two viscous drops[END_REF][START_REF] Yiantsios | On the buoyancy-driven motion of a drop towards a rigid surface or a deformable interface[END_REF][START_REF] Nemer | Coalescence of drops with mobile interfaces in a quiescent fluid[END_REF]. The film thickness decreases exponentially over time during most of the drainage process, the decay rate depending on whether the Reynolds number Re a = ρV a R/µ based on the relative approach velocity V a of the bubbles is large or small [START_REF] Chesters | The modelling of coalescence processes in fluid-liquid dispersions: a review of current understanding[END_REF]. This exponential thinning law was confirmed experimentally [START_REF] Debrégeas | The life and death of "bare" viscous bubbles[END_REF]) and numerically [START_REF] Pigeonneau | Low-Reynolds-number gravity-driven migration and deformation of bubbles near a free surface[END_REF]) by considering buoyancy-driven bubbles reaching a free surface. This law holds as far as the film thickness can be considered uniform. However, similar to the case of immobile or partially mobile interfaces, the minimum film thickness initially located on the line of centres shifts gradually to the film periphery, giving rise to a 'dimple' corresponding to the transition region between the film and the outer flow. To take the influence of this dimple into account, [START_REF] Chesters | Bubble coalescence in pure liquids[END_REF] solved numerically the set of thinning equations with appropriate boundary conditions for initially spherical bubbles, assuming Re a 1, i.e. considering that inertial and capillary effects are in balance. They observed that, when the minimum film thickness has reduced sufficiently, the thinning velocity at the dimple position levels off at a value close to 0.1V a . Under constant-force conditions, the crucial consequence of this thinning evolution is that drainage is completed in a finite time without the need for a non-hydrodynamic force to intervene. Defining the approach Weber number W e a = ρV 2 a R/γ, the dimensionless inertial drainage time is then

T di = k i W e a /V a , ( B 1)
with k i ≈ 1.0 and V a = V a /(gR) 1/2 in the buoyancy-driven case of interest here. The value of k i was later slightly re-evaluated to k i ≈ 1.08 by [START_REF] Duineveld | Bouncing and coalescence of two bubbles in water[END_REF][START_REF] Duineveld | Bouncing and coalescence of bubble pairs rising at high Reynolds number in pure water or aqueous surfactant solutions[END_REF]. Within a liquid film bounded by two clean gas-liquid interfaces, the London-van der Waals force is attractive, thus shortening the coalescence time. However this force is known to be significant only when the distance between the two interfaces is less than 100 nm, which made Chesters & Hofman conclude that it barely shortens the coalescence process.

Obviously, coalescence occurs only if the relative approach velocity remains positive throughout the drainage. For this to be the case, only part of the kinetic energy resulting from the relative motion of the two bubbles (moving with velocities ±V a /2) must be converted into surface energy through the deformation of the bubble-fluid interface in the near-contact region. Still in the limit Re a 1, this criterion yields a critical Weber number W e ac beyond which the two bubbles bounce once or several times, until eventually coalescing when the approach Weber number has decreased sufficiently. Under potential flow assumptions, the kinetic energy associated with the bubble relative motion is proportional to the virtual mass (or added-mass) coefficient C M a in the corresponding direction, making W e ac vary linearly with C M a . For nearly-spherical bubbles, [START_REF] Chesters | Bubble coalescence in pure liquids[END_REF] established that the relative increase of the surface energy is at leading order (k i W e a /4) 2 , which yields W e ac = 2 3 k -2 i C M a . For two spheres in contact after having moved toward each other , C M a ≈ 0.80 [START_REF] Voinov | On the motion of two spheres in a perfect fluid[END_REF][START_REF] Miloh | Hydrodynamics of deformable contiguous spherical shapes in an incompressible inviscid fluid[END_REF], so that W e ac ≈ 0.45 [START_REF] Duineveld | Bouncing and coalescence of two bubbles in water[END_REF][START_REF] Duineveld | Bouncing and coalescence of bubble pairs rising at high Reynolds number in pure water or aqueous surfactant solutions[END_REF].

In addition, Chesters & Hofman pointed out that previous results also apply to a single bubble reaching a free surface, up to a simple geometric transformation. The outcome is that in this configuration (B 1) transforms into T di = 4k i W e a /V a , while the above criterion for the onset of bouncing becomes W e ac = 1 3 k -2 i C M a . Assuming the bubble to be spherical and the Froude number of the free surface, V 2 a , to be large, the relevant virtual mass coefficient is C M a ≈ 0.42 [START_REF] Miloh | Hydrodynamics of deformable contiguous spherical shapes in an incompressible inviscid fluid[END_REF], so that W e ac ≈ 0.12 [START_REF] Duineveld | Bouncing and coalescence of two bubbles in water[END_REF]. By tracking sub-millimeter size bubbles rising up to a free surface in ultrapure water, the same author determined the experimental threshold as W e ac ≈ 0.105 (corresponding to a rise Reynolds number Re ≈ 50). This agreement with the theoretical prediction is supported by more recent experiments [START_REF] Vakarelski | Free-rising bubbles bounce more strongly from mobile than from immobile water-air interfaces[END_REF] and provides an important support to the approach of [START_REF] Chesters | Bubble coalescence in pure liquids[END_REF], although some of its aspects have been questioned [START_REF] Chan | Film drainage and coalescence between deformable drops and bubbles[END_REF].

The above conclusions hold for nearly-spherical bubbles provided effects of the liquid viscosity are negligible. However, despite the uniform velocity profile in the film, viscous effects arise through normal stresses. For a given film thickness, these effects tend to increase the film radius, i.e. the area of the near-contact region. This results in a significant decrease of the film thinning rate as soon as Re a 10 ( [START_REF] Chesters | Bubble coalescence in pure liquids[END_REF]. In this Re a -range, this thinning rate decreases continually as the drainage proceeds, making the London-van der Waals force inescapable for coalescence to occur.

For low enough Re a , the film evolution is governed by a viscous-capillary balance, hence by the capillary number Ca a = W e a /Re a = µV a /γ. To the best of our knowledge, no theoretical model is available to predict the drainage time in this viscosity-dominated regime for fully mobile interfaces. However, recent experimental data may be used to obtain an empirical scaling law from which realistic estimates may be inferred. The case of air bubbles coalescing at a clean free surface after rising in a liquid 20 times more viscous than water was considered in detail by [START_REF] Vakarelski | Coalescence dynamics of mobile and immobile fluid interfaces[END_REF]. Small bubbles (R 0.1 mm) were observed to coalesce almost immediately on their arrival at the surface, while larger bubbles remained 'glued' there during some time without bouncing, until coalescence eventually occurred. A single tiny bounce was detected for bubbles larger than R ≈ 0.45 mm, corresponding to W e a 0.135. Overall, they found the residence (i.e. drainage) time T dv of the bubble at the surface, to be proportional to R 2 . Their results may be used as a basis to derive a generic empirical expression for the viscous drainage time, T dv . Considering that T dv = (g/R) 1/2 T dv is primarily driven by the capillary number in this O(1)-Reynolds number range and keeping in mind that the rise speed also grows like R 2 in this regime, it turns out that the simplest admissible scaling law for T dv is

T dv = k v Ca 3/2 a /V a , ( B 2)
with k v = 1.8 × 10 3 according to the above experimental data. It is then a simple matter to compare the inertial and viscous estimates for the drainage time in a given fluid and for a given bubble size. For instance, in ultrapure water (M o = 2.6 × 10 -11 ), the bounce/no bounce threshold for a bubble reaching a free surface is known to correspond to 0.335 mm-radius bubbles [START_REF] Duineveld | Bouncing and coalescence of two bubbles in water[END_REF]). Under such conditions, (B 1) (with the appropriate transformation) and (B 2) yield T di ≈ 0.32 and T dv ≈ 0.064, respectively. Hence the drainage is controlled by inertial effects and coalescence takes place very soon after the bubble collides with the free surface. Similarly, in the viscous liquid used by [START_REF] Vakarelski | Coalescence dynamics of mobile and immobile fluid interfaces[END_REF] (M o = 6.6 × 10 -5 ), the same two predictions for a bubble with R = 0.45 mm (i.e. just below the bounce/no bounce threshold) yield T di ≈ 0.68 and T dv ≈ 23.8, respectively. In this case, the drainage time estimated through the viscous law (B 2) is 35 times larger than that predicted by the inviscid approach, implying that the latter is irrelevant. This corresponds to a situation in which the bubble stays 'glued' to the free surface during a long time before coalescing. For a pair of bubbles rising in line under similar conditions, the two bubbles form a compound 'dumbbell' bubble during the drainage process, as observed by [START_REF] Sanada | Interactions and coalescence of bubbles in stagnant liquid[END_REF] and Watanabe & Sanada (2006). To estimate T dv in this case, (B 2) must be modified according to the geometric transformation of Chesters & Hofman, which yields

T dv = 1 2 k v Ca 3/2
a /V a . Most results reviewed above were obtained with nearly-spherical bubbles (Bo 1, in practice Bo 0.2). Influence of a significant distortion of the bubble shape on the coalescence process is complex because it involves antagonistic effects. On the one hand, the curvature of the near-pole region of an oblate bubble is smaller than that of a spherical bubble, making the area of the near-contact region larger. Thus the radial position of the dimple shifts outward, and a longer time is required to squeeze the film with a given approach velocity. [START_REF] Duineveld | Bouncing and coalescence of two bubbles in water[END_REF] solved the set of inviscid thinning equations for oblate bubbles with various aspect ratios. Compared to the reference case, his results show that the pre-factor k i involved in the prediction (B 1) for the drainage time is increased approximately by a factor of 2 (resp. 3) for χ = 1.5 (resp. χ = 2). On the other hand, an oblate body moving along its short axis displaces more fluid than a sphere, which translates into a larger virtual mass coefficient, hence a larger kinetic energy available for the drainage. Moreover, as explained in § 4, the bubble oblateness enhances wake effects, increasing the amount of fluid displaced by the bubble through the entrainment process in the wake. Therefore, for a given approach velocity, the kinetic energy of the fluid displaced by oblate bubbles may be significantly larger than the estimate based on the simple irrotational added-mass concept. Unfortunately, no drainage time prediction incorporating wake entrainment and/or viscous effects seems available for such bubbles.

Appendix C. Lift reversal mechanisms on distorted bubbles

As is well known, the shear-induced lift mechanism mentioned in § 4, hereinafter called L-mechanism, stems from the bending of the vorticity of the carrying shear flow past the bubble, a process that results in the formation of a pair of counter-rotating streamwise vortices in its wake [START_REF] Legendre | The lift force on a spherical bubble in a viscous linear shear flow[END_REF]. The vorticity produced at the bubble surface plays no role in this mechanism which is inviscid by nature (Lighthill 1956;[START_REF] Auton | The lift force on a spherical body in a rotational flow[END_REF]). However, for reasons discussed in § 4, the strength of the surface vorticity resulting from finite-Re effects increases sharply with the bubble oblateness. When an oblate bubble rises in a fluid at rest, the amount of vorticity produced at its surface becomes large enough for the axisymmetric wake to become unstable within a finite range of Reynolds number, Re S-(χ) Re Re S+ (χ), provided the bubble aspect ratio exceeds a critical value χ cS ≈ 2.2 [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF]. This mechanism is at the root of the path instability of millimeter-size bubbles rising in water [START_REF] Mougin | Path instability of a rising bubble[END_REF]. Similar to the above L-mechanism, it gives rise to a wake dominated by a pair of counter-rotating streamwise vortices, the sign of which is selected by some initial disturbance. This wake instability mechanism, hereinafter called S-mechanism, subsists when a strongly oblate bubble rises in a weak shear flow. The only difference is that the initial disturbance is now provided by the outer shear, so that the sign of the streamwise vorticity in each trailing vortex is no longer random. Rather, a detailed analysis reveals that the contributions of the L-and S-mechanisms to the vortex tilting term involved in the streamwise vorticity balance have opposite signs [START_REF] Adoua | Reversal of the lift force on an oblate bubble in a weakly viscous linear shear flow[END_REF]. For this reason, the sign of the overall sideways force depends on the relative strength of the two mechanisms. If the outer shear is strong enough, the L-mechanism dominates even for χ > χ cS and the lift force keeps the sign it would have for a spherical bubble. In contrast, if the shear is weak enough and Re and χ fall in the range where the Smechanism is active, the latter becomes dominant, yielding a reversed lift force. The lower Reynolds number beyond which lift reversal governed by the S-mechanism takes place in a weak shear flow is a sharply decreasing function of χ -χ cS , with Re S-≈ 155 for χ = χ cS and Re S-≈ 55 for χ = 2.5 for instance. Conversely, the upper Reynolds number beyond which the lift force recovers the sign predicted by the L-mechanism dramatically increases with χ -χ cS , from Re S+ = Re S-≈ 155 for χ = χ cS to Re S+ ≈ 680 for χ = 2.5 for instance [START_REF] Adoua | Reversal of the lift force on an oblate bubble in a weakly viscous linear shear flow[END_REF]. Because of this mechanism, one can suspect that bubbles experiencing a sufficient deformation may render the in-line configuration stable. Indeed, if the S-mechanism dominates, any deviation of the TB from the wake axis is expected to be counteracted by the reversed lift force.

A distinct mechanism may also lead to the same effect on a non-axisymmetric bubble (more generally a drop). This mechanism, hereinafter called A-mechanism, is a consequence of the asymmetric deformation experienced by a drop immersed in a shear flow. In the zero-Re limit, the drop is known to deform in such a way that its major and minor axes align with the eigen-directions of the associated strain rate [START_REF] Taylor | The viscosity of fluids containing small drops of another fluid[END_REF]. In the case the drop has an additional translation with respect to the fluid, this deformation induces a nonzero sideways force, even at Re = 0. When the Reynolds number is finite, this deformation-induced transverse force combines with the inertial shear-induced lift.

However, theoretical predictions indicate that the two mechanisms result in sideways forces of opposite signs, except for drops whose viscosity is close to that of the suspending fluid (see equation ( 49) in [START_REF] Magnaudet | Drag, deformation and lateral migration of a buoyant drop moving near a wall[END_REF]). The deformation-induced (resp. inertia-induced) force being proportional to the Weber (resp. Reynolds) number, the direction of the total lift force is governed by the capillary number Ca = W e/Re. This force changes sign for a critical value Ca = Ca c which, for a bubble with negligible internal viscosity, depends only on the relative shear rate βR/u T , β denoting the ambient shear rate. While the total lift force keeps the sign corresponding to the inertial shearinduced mechanism if Ca < Ca c , it acts in the opposite direction for larger Ca. In the case of a bubble, this mechanism results from the fact that the non-penetration condition and the normal stress balance have to be jointly satisfied at the gas-liquid interface. This requirement obviously subsists for Reynolds numbers larger than unity and so does the above mechanism, although inertial effects tending to make the bubble oblate combine with those of the outer shear to produce more complex asymmetric shapes compared to the low-Re configuration. Lift reversal due to the A-mechanism has been observed in computations [START_REF] Ervin | The rise of bubbles in a vertical shear flow[END_REF][START_REF] Sankaranarayanan | Lift force in bubbly suspensions[END_REF] and experiments performed with isolated bubbles rising in viscous liquids sheared in a Couette device [START_REF] Tomiyama | Transverse migration of single bubbles in simple shear flows[END_REF][START_REF] Aoyama | Lift force acting on single bubbles in linear shear flows[END_REF].

Figure 1 .

 1 Figure 1. Sketch of the problem. (a) initial configuration; (b) definition of some geometric parameters used to characterize the relative position of the two bubbles.

Figure 2 .

 2 Figure 2. Detail of a typical grid for a bubble pair with Ga = 30 and Bo = 0.3 at t = 8. (a): grid refinement in the wake and boundary layer regions; four refinement levels are shown, from dark blue (∆ = R/8.5) to yellow (∆ = ∆min = R/68). (b): grid detail within the leading bubble and in the neighbourhood of its surface. The gas-liquid interface is marked with a black (resp. red) line in (a) (resp. (b)).

Figure 3 .

 3 Figure 3. Topology-based AMR strategy employed to refine the grid from time t = t1 in (a) to t = t1 + ∆t in (b) and (c), as the liquid film separating the two bubbles gets squeezed.

Figure 4 .

 4 Figure 4. Final bubble shapes and separations observed in the axisymmetric configuration.Red bubbles maintain a finite separation, while blue bubbles are about to coalesce.

Figure 5 .

 5 Figure 5. Variation of the final equilibrium distance Se against: (a) the Reynolds number; and (b) the Weber number. Both Re and W e are based on the final rise velocity of the bubble pair.Each series identified with a given color corresponds to the same Ga and different Bo (increasing from top to bottom). In (a), the open circles correspond to Bo = 0.005, while the solid line and open triangles refer to the prediction (5.1) and the numerical results of[START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF] for spherical bubbles, respectively. The dotted lines in (a) and (b) correspond to the prediction (5.2), while the dash-dotted line in (a) represents this prediction evaluated for W e = 0.

  , namely S e (Re) = 4.40 log 10 Re -3.06 (5.1) As figure 5(a) reveals, present low-Bo predictions are in excellent agreement with (5.1) for Ga 20. For lower Ga, the numerical values of S e (Ga, Bo = 0.005) are found to lie slightly below the prediction (5.1). This is in line with the results of Hallez & Legendre (2011) which indicate for instance that S e (Re = 20) is approximately 5% less than predicted by (5.1). Figure 5(a) shows that present results for Bo = 0.005 (open circles) are in excellent agreement with those of the latter authors (open triangles) throughout the range of Reynolds number considered here. This establishes the accuracy of present computations in the in-line configuration.

Figure 6 .

 6 Figure 6. Set of configurations encountered in the three-dimensional simulations. (a) (Ga, Bo) phase diagram; (b) (Re, χ) phase diagram (for each (Ga, Bo) pair, Re and χ are the steady-statevalues determined with the corresponding isolated bubble); (c) typical trajectories illustrating each of the three configurations. Bullets: DKT scenario; solid squares: ASE scenario; solid triangles: collision followed by coalescence. In (c), the three panels from left to right correspond to (Ga, Bo) = (10, 0.1), (30, 0.3) and (20, 0.5), respectively. The thin solid lines in (a) are the iso-M o lines corresponding to different liquids, with water at the very bottom, then silicone oils T0-T11 of increasing viscosity from bottom to top (e.g.[START_REF] Zenit | Path instability of rising spheroidal air bubbles: a shapecontrolled process[END_REF] for the corresponding physical properties). The dashed line in (a) and (b) separates the subdomain of weakly deformed bubbles in which a finite equilibrium separation is reached in the axisymmetric case, from that in which the bubbles eventually coalesce.

Figure 8 .

 8 Figure 8. Path of a bubble pair following a DKT scenario (Ga = 10, Bo = 0.02). (a) Side and top views, with r the radial distance to the initial path and numbers referring to the dimensionless time at the corresponding position; (b) three-dimensional streamlines past the bubbles at different instants of time, in the reference frame of the leading bubble. A zoom of the flow field in the gap is provided for t = 30 and t = 38, to highlight the inception of the non-axisymmetric fluid motion.

Figure 9 .

 9 Figure 9. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (10, 0.02) undergoing a DKT interaction. (a): vertical velocity component of the LB (solid red line) and TB (dashed blue line); (b) same for the horizontal component; (c): vertical (red line, left axis) and horizontal (blue line, right axis) separations; (d) inclination of the line of centres with respect to the vertical. Open squares and circles in (a) and (c): axisymmetric prediction.

Figure 10 .

 10 Figure 10. Path of a bubble pair following an ASE scenario (Ga = 30, Bo = 0.3). (a) Side and top views of the path with numbers referring to the dimensionless time at the corresponding position of the TB; (b) three-dimensional streamlines past the bubbles at different instants of time, in the reference frame of the leading bubble.

Figure 11 .

 11 Figure 11. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (30, 0.3) undergoing an ASE interaction. (a): vertical velocity component of the LB (solid red line) and TB (dashed blue line); (b) same for the horizontal component; (c): vertical (red line, left axis) and horizontal (blue line, right axis) separations; (d) inclination of the line of centres with respect to the vertical. Open squares and circles in (a) and (c): axisymmetric prediction.

Figure 13 .Figure 14 .

 1314 Figure 13. Evolution of some characteristics of a bubble pair with Ga = 20, Bo = 0.5 undergoing coalescence. (a): vertical (solid lines, left axis) and horizontal (dash-dotted lines, right axis) components of the bubble velocity, with the red, blue and black lines corresponding to the LB, TB and final bubble, respectively; (b): vertical (red line, left axis) and horizontal (blue line, right axis) components of the separation. Open squares and circles: axisymmetric prediction. In (a), the shape of the resulting bubble is shown at several successive time instants during the transient 24.5 t 30 following coalescence. (a) (b)

7. 1 .

 1 Figure15shows how several characteristics of the path of bubble pairs with Ga = 10 vary with the Bond number in the range 0.02 Bo 0.2, i.e. 8 × 10 -10 M o 8 × 10 -7 . The DKT scenario is observed for the lower three values of Bo, while the pair with Bo = 0.2 eventually experiences head-on coalescence. For lower Bo, the interaction always leads to the side-by-side configuration and the lateral separation stabilizes at a value close to 2.5. While the Bond number does not have any noticeable influence on the final path in this regime, it affects the critical time by which the tumbling process sets in. Owing to the connection between the bubble oblateness and the amount of vorticity produced at its surface, the more oblate the LB is the stronger the attractive wake effect is. Therefore, at a given time, the larger Bo the shorter the separation between the two bubbles. From figures 15(b) -(d), it may be inferred that the axial symmetry of the

Figure 15 .

 15 Figure 15. Variations with the Bond number of some path characteristics for Ga = 10. (a) front view of the path (the solid and dashed lines correspond to the LB and TB, respectively); (b) vertical separation; (c) horizontal separation; (d) inclination of the line of centres. The square symbol denotes the position/time at which coalescence takes place for Bo = 0.2.

Figure 16 .

 16 Figure 16. Variations with Bo of some path characteristics for Ga = 30. (a) front view of the path (the solid and dashed lines correspond to the LB and TB, respectively); (b) vertical separation; (c) horizontal separation; (d) inclination of the line of centres. In (a), the square symbol indicates the position/time at which coalescence takes place for Bo = 0.5.

Figure 17 .

 17 Figure 17. Evolution of some characteristics of a bubble pair with (Ga, Bo) = (30, 0.45). (a): vertical (red line) and horizontal (blue line) separation distances; (b) inclination of the line of centres.

Figure 18 .

 18 Figure 18. Evolution of the bubble shapes and relative positions for near-critical conditions. (a): Ga = 30; (b): Ga = 15.

Figure 19 .Figure 20 .

 1920 Figure 19. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (30, 0.3) undergoing initial angular deviations θ0 from 0 • to 2 • . (a): vertical velocity of the LB (solid line) and TB (dashed line); (b) horizontal velocity of the TB; (c): vertical (solid line, left axis) and horizontal (dashed line, right axis) separations; (d) inclination of the line of centres. In (b), the bullets on the Sr-curve indicate the location where the maxima of Vr are reached.

Figure 21 .

 21 Figure 21. Phase diagram showing whether a bubble pair released with an angular deviation θ0 = 2 • evolves in the ASE regime (squares) or eventually coalesces (triangles). (a) (Ga, Bo) map; (b) (Re, χ) representation based on the steady-state values corresponding to the rise of the corresponding isolated bubble.

Figure 22 .Figure 23 .

 2223 Figure 22. Evolution of the bubble shapes and relative positions at Ga = 30 under near-critical Bo-conditions in the presence of an initial angular deviation θ0 = 2 • .

Figure 24 .

 24 Figure 24. Influence of the initial separation S0 on several characteristics of a bubble pair with (Ga, Bo) = (30, 0.3). (a): vertical velocity component of the LB (solid line) and TB (dashed line); (b) same for the horizontal component; (c): vertical (solid line) and horizontal (dashed line) separations; (d) inclination of the line of centres. Red, blue, green and orange lines correspond to S0 = 6, 8, 10 and 12, respectively.

Figure 25 .

 25 Figure 25. Influence of the sequential release of the two bubbles. (a) vertical separation; (b) lateral velocity; (c) inclination angle. Solid line: bubbles are released simultaneously with S0 = 8; blue dash-dotted line: same with S0 = 10 and a time shift t0 = 2.6; red dashed line: the TB is released from rest after the LB has travelled a distance S(t0) = 8.

Figure 26 .Figure 27 .

 2627 Figure 26. Influence of the tolerance imposed to the Poisson solver on the evolution of a bubble pair with Ga = 30, Bo = 0.3. (a): vertical velocity component of the LB (red line) and TB (blue line); (b) vertical (red line, left axis) and horizontal (blue line, right axis) components of the separation. Dotted, solid and dash-dotted lines refer to simulations performed with a tolerance of 1 × 10 -2 , 1 × 10 -4 and 1 × 10 -6 , respectively.

Table 1 .

 1 Dimensionless parameters characterizing the system. uT is the terminal rise speed of the bubble, and the subscript l refers to the properties of the liquid. may be chosen among those listed in table 1. The Reynolds (Re) and Weber (W e) numbers are generally preferred in theoretical studies. In contrast, in experiments and computations like those discussed here, the terminal rise speed u T is unknown a priori. This is why the Galilei (or Archimedes) number (Ga) and the Bond number (Bo) are usually selected and are used throughout the present study. The Morton number (M o) is frequently used in place of the Bond number, since M o = Bo 3 /Ga 4 . In what follows,

	Name	Abbreviation	Expression
	Reynolds number	Re	ρ l uT R/µ l
	Weber number	W e	ρ l u 2 T R/γ
	Galilei number	Ga	ρ l g 1/2 R 3/2 /µ l
	Bond number	Bo	ρ l gR 2 /γ
	Morton number	M o	gµ 4 l /ρ l γ 3

† Email address for correspondence: mjni@ucas.ac.cn ‡ Email address for correspondence: Jacques.Magnaudet@imft.fr

Acknowledgements.

The authors acknowledge the partial support of this research by NSFC ( Natural Science Foundation of China) under grants #51636009, #11872296, and #U 173227. J. Z. acknowledges the Young Elite Scientists Sponsorship Program under grant YESS No. 2018QNRC001.

Declaration of interests.

The authors report no conflict of interest.

Appendix A. Numerical tests

To assess the accuracy of the three-dimensional computations, two series of tests were performed, both with the parameters (Ga = 30, Bo = 0.3) corresponding to the ASE configuration discussed in § 6.3. To properly interpret these tests, it must be kept in mind