Jie Zhang

Ming-Jiu Ni

Jacques Magnaudet

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

The buoyancy-driven dynamics of a pair of gas bubbles released in line is investigated numerically, focusing on highly inertial conditions under which isolated bubbles follow non-straight paths. In an early stage, the second bubble always drifts out of the wake of the leading one. Then, depending on the ratios of the buoyancy, viscous and capillary forces which define the Galilei (Ga) and Bond (Bo) numbers of the system, five distinct regimes specific to such conditions are identified, in which the two bubbles may rise independently or continue to interact and possibly collide in the end. In the former case, they usually perform large-amplitude planar zigzags within the same plane or within two distinct planes, depending on the oblateness of the leading bubble. However, for large enough Ga and low enough Bo, they follow nearly vertical paths with small-amplitude erratic horizontal deviations. Increasing Bo makes the wake-induced attraction toward the leading bubble stronger, forcing the two bubbles to realign vertically one or more times along their ascent. During such sequences, wake vortices may hit the trailing bubble, deflecting its path and, depending on the case, promoting or hindering further possibilities of interaction. In some regimes, varying the initial distance separating the two bubbles modifies their lateral separation beyond the initial stage. Similarly, minute initial angular deviations favour the selection of a single vertical plane of rise common to both bubbles. These changes may dramatically affect the fate of the tandem as, depending on Bo, they promote or prevent future vertical realignments.

Introduction

In the first part of this investigation (Zhang, Ni & Magnaudet (2021), hereinafter referred to as ZNM), we analyzed the results of a series of simulations revealing the mechanisms governing the hydrodynamic interactions between two deforming gas bubbles released in line in a liquid at rest. The physical parameters were selected in such a way that the bubbles rose at moderate Reynolds number, and each bubble taken separately would ascend in a straight line. However, millimeter-size air bubbles rising in low-viscosity liquids, most notably in water, are subject to path instability. Consequently, they usually follow either zigzagging planar paths or more or less flattened spiralling paths, with in both cases large-amplitude horizontal excursions. These are the regimes on which this second part focuses.

As discussed in ZNM, interactions between two neighbouring bubbles (so-called pair interactions) play a key role in the microstructure of bubbly suspensions, as they govern to a large extent the bubble distribution and the agitation in the carrying liquid. Early computational studies of buoyancy-driven bubbly suspensions [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 the bubbles to keep a spherical shape and disregarded any possible influence of viscosity. Since the potential flow approximation predicts that two bubbles rising in line repel each other while they attract each other when rising side by side, such simulations inescapably concluded that the suspension dynamics lead to the formation of large horizontal bubble clusters. The next generation of simulations addressed more realistic conditions by considering the full Navier-Stokes equations, possibly including surface tension effects. With Reynolds numbers of O(10) to O(100) and spherical or weakly deformed bubbles, these simulations confirmed the tendency of bubble pairs to align horizontally, albeit less clearly than in the potential flow approximation [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows Part 2. Moderate 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] Esmaeeli | A direct numerical simulation study of the buoyant rise of bubbles at O(100) Reynolds number[END_REF][START_REF] Yin | Velocity fluctuations and hydrodynamic diffusion in finite-Reynolds-number sedimenting suspensions[END_REF]. With significantly deformed bubbles, the dynamics of bubbly suspensions was observed to depend crucially on the Reynolds number. More specifically, the simulations of [START_REF] Esmaeeli | A direct numerical simulation study of the buoyant rise of bubbles at O(100) Reynolds number[END_REF] with Reynolds numbers of O(100) revealed quite homogeneous bubble distributions, while for Reynolds numbers of O(10), marked vertical bubble alignments were noticed by [START_REF] Bunner | Effect of bubble deformation on the properties of bubbly flows[END_REF], suggesting that some 'chimney effect' is at work under such conditions. It must be pointed out that, for technical reasons, all of the above studies made use of some simplifying assumptions which depart from common physical conditions. In most of them, the liquid-to-gas density ratio was reduced by a factor of 20 to 50, making the bubble motion less sensitive to small fluctuations of the carrying liquid than under standard experimental conditions, while in others [START_REF] Yin | Velocity fluctuations and hydrodynamic diffusion in finite-Reynolds-number sedimenting suspensions[END_REF] bubbles were not allow to deform. These restrictions were removed in the most recent simulations [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF][START_REF] Innocenti | Direct numerical simulation of bubble-induced turbulence[END_REF]. However, the parameters were still chosen in such a way that bubbles are not subject to path instability. That is, the liquid viscosity and/or surface tension were sufficiently large to prevent an isolated bubble from rising along a zigzagging or spiralling path, making the predictions barely representative of most experiments performed in water. Even with the largest inertia-to-viscosity force ratio considered to date [START_REF] Innocenti | Direct numerical simulation of bubble-induced turbulence[END_REF], surface tension was selected so as to keep the bubbles only mildly deformed, preventing the occurrence of path instability. About one billion grid points were employed to track the flow generated by the rise of 256 bubbles in that study, but a ten times larger grid would be required to deal with regimes involving zigzagging or spiralling bubbles, owing to the very thin boundary layers and complicated wake structures involved.

Despite these demanding requirements, detailed numerical investigations of this highly inertial regime are needed, in view of its specificities and its ubiquity in natural and engineering air-water bubbly flows. For instance, it is known that the liquid velocity fluctuations generated by the rise of a dilute bubble swarm are highly anisotropic when the bubbles are only slightly distorted, the vertical fluctuations having a variance 4 -5 times larger than their horizontal counterpart [START_REF] Zenit | Measurements of the average properties of a suspension of bubbles rising in a vertical channel[END_REF]. In contrast, this ratio falls to values close to 2 for larger bubbles exhibiting path instability [START_REF] Risso | Velocity fluctuations in a homogeneous dilute dispersion of high-Reynolds-number rising bubbles[END_REF][START_REF] Riboux | Experimental characterization of the agitation generated by bubbles rising at high Reynolds number[END_REF]. Obviously, the reason for this sharp decrease stands in the large horizontal fluctuations of bubble positions along a zigzagging or spiralling path, which in turn induce large horizontal fluctuations in the liquid velocity. More importantly, the fraction of the liquid agitation resulting from the interacting bubble wakes is much larger in the case of zigzagging/spiralling bubbles. This contribution makes the whole carrying flow exhibit genuine turbulent properties over a significant range of scales. The mixing ability of this 'bubble-induced turbulence' is well highlighted by examining how a weakly diffusive species (e.g. dye) or a temperature gradient is mixed in a bubble swarm [START_REF] Alméras | Mixing by bubble-induced turbulence[END_REF][START_REF] Gvozdic | Experimental investigation of heat transport in homogeneous bubbly flow[END_REF]. In both cases, a clear increase of the relevant transfer coefficient with the gas volume fraction is observed in very dilute swarms (i.e. gas volume fractions of the order of 1%), together with a dramatic increase in the strength of the scalar fluctuations.

Although still far from such turbulent bubbly suspensions and the key statistical information gained from the aforementioned studies, the simplest arrangement capable of providing detailed insight into wake interactions in the relevant highly inertial regime is that involving a pair of bubbles. Focusing on such an elementary configuration offers a complementary point of view with respect to many-bubble configurations, since the hydrodynamic mechanisms involved in the interaction process can be identified and analyzed in a deterministic manner, and the few control parameters can be varied separately to check their influence. We refer the reader to the introduction of ZNM for a review of the available knowledge regarding the dynamics of pairs of spherical or weakly distorted bubbles in moderately inertial regimes as well as in the potential flow approximation. In the regime of interest here, some experimental studies considered the side-by-side arrangement [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]. In this configuration, the potential flow approximation predicts that the two bubbles are attracted toward each other. However, when they get close enough, their wakes interact directly and the sign of the transverse force acting on each bubble may reverse. This is why one objective of these studies was to check the predictions of the inviscid theory of [START_REF] Chesters | Bubble coalescence in pure liquids[END_REF] regarding the conditions under which the two bubbles bounce (possibly repeatedly) or rather coalesce during their ascent. It must be noticed that none of these studies addressed quantitatively the detail of wake interactions. Only one of them [START_REF] Sanada | Motion and coalescence of a pair of bubbles rising side by side[END_REF] visualized the wakes through a photochromic technique and could connect the reversal of the transverse bubble motion to the encounter of the two wakes.

Some experiments were also carried out in the in-line configuration on which the present investigation focuses, the most recent one being that of [START_REF] Kusuno | Lift force acting on a pair of clean bubbles rising in-line[END_REF]. However, all these studies considered Reynolds numbers of O(10 -100) which belong to the moderately inertial regime examined numerically by ZNM and, independently but with the same code, by [START_REF] Kusuno | Wake-induced lateral migration of approaching bubbles[END_REF]. A noticeable exception is the work of [START_REF] Filella | Interaction of two oscillating bubbles rising in a thin-gap cell: vertical entrainment and interaction with vortices[END_REF] in which wake interactions behind two bubbles rising in a thin-gap cell were scrutinized using time-resolved particle image velocimetry. However the corresponding two-dimensional dynamics makes the problem barely comparable to the three-dimensional configuration of interest here. The main outcome of the recent studies by [START_REF] Kusuno | Lift force acting on a pair of clean bubbles rising in-line[END_REF], [START_REF] Kusuno | Wake-induced lateral migration of approaching bubbles[END_REF] and ZNM is the existence of three regimes with markedly different dynamics, according to the values of the control parameters and to the detail of the initial conditions. In short, pairs of nearly-spherical bubbles released exactly in line and rising with a Reynolds number of O(10) follow the Drafting-Kissing-Tumbling (DKT) scenario widely observed with sedimenting spherical 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]. Conversely, pairs of bubbles exhibiting a sufficient oblateness always collide and eventually coalesce, most of the time in the headon configuration. Below this critical, Reynolds-number dependent, oblateness and beyond the narrow parameter range where the DKT mechanism takes place, the two bubbles follow an Asymmetric Side Escape (ASE) scenario in which the trailing bubble leaves the wake of the leading bubble through a vigorous lateral drift. Beyond this crucial stage, the trailing bubble rises along a new nearly vertical path, whereas the leading bubble is barely disturbed by the interaction and essentially goes on rising along its initial path, with however some marginal inclination. The tandem stabilizes in a configuration whose final geometry, i.e. lateral separation and inclination of the line of centres with respect to the vertical, depends significantly on the control parameters and is very sensitive to initial conditions. Since new bubble pairs form continuously in a swarm, the DKT and ASE scenarios are self-repeating in a real suspension. Given that the latter is by far the most frequent in the moderately inertial regime and yields significantly different bubble pair geometries according to tiny changes in the initial conditions, it is no surprise that bubbly suspensions with the corresponding characteristics are far more homogeneous than predicted by the simplistic potential flow approximation.

Still considering the in-line configuration, our aim in this second part is to explore the regime in which bubbles rise with Reynolds numbers of O(100 -1000). In most cases, an isolated bubble would then either perform large-amplitude planar zigzags or follow a spiralling path. Similar to ZNM, we vary the control parameters and the initial conditions so as to cover a significant range of physical situations and analyze different evolution scenarios, with a specific attention on those involving the direct interaction between the trailing bubble and the three-dimensional time-dependent wake of the leading bubble. We present the problem and summarize the numerical approach in § 2 (a specific test aimed at confirming that the grid resolution is adequate even for the highest Reynolds numbers considered is discussed in Appendix A). Then the results obtained under nominal initial conditions are discussed in § 3. Influence of the initial conditions, i.e. a slight angular deviation or a change in the vertical separation, is discussed in § 4. A summary of the main findings and some prospects are presented in § 5.

Problem statement and outline of the numerical approach

Since the problem was introduced in detail in ZNM, only a brief account is presented here. A pair of initially spherical gas bubbles with radius R is released in line with a vertical centre-to-centre separation S 0 (figure 1(a)). Starting from rest, the two bubbles rise freely under the effect of buoyancy. The evolution of the tandem geometry is tracked by considering the vertical and horizontal dimensionless separations S(t) = S(t)/R and S r (t) = S r (t)/R, respectively, and the angular deviation of the line of centres with respect to the vertical, θ(t) (figure 1(b)). In addition to the initial separation S 0 = S 0 /R (hereinafter set to S 0 = 8 except specified otherwise), the flow and bubble dynamics are characterized by the Galilei and Bond numbers respectively defined as

Ga = ρ l g 1/2 R 3/2 /µ l , Bo = ρ l gR 2 /γ , (2.1)
where ρ l and µ l are the density and viscosity of the carrying liquid, γ is the surface tension, and g denotes gravity. Once the terminal velocity u T of each bubble is known, the terminal Reynolds and Weber numbers, Re = ρ l u T R/µ l and W e = ρ l u 2 T R/γ, may be determined. Moreover, the Morton number M o = Bo 3 /Ga 4 = gµ 4 l /ρ l γ 3 is useful to identify the dynamics of bubble pairs in a given fluid, irrespective of the bubble size. In ZNM, Ga and Bo were varied in the range 10 Ga 30 and 0.01 Bo 1.0, respectively. The corresponding terminal Reynolds numbers were such that 10 Re 120, so that an isolated bubble followed a rectilinear path. Here we consider the parameter range 30 < Ga 90 and 0.02 Bo 1.0. In that range, the shape of an isolated bubble is close to an oblate spheroid. For instance, a 1.8 mm-diameter air bubble rising in water has a Bond number Bo ≈ 0.11 and a Galilei number Ga ≈ 85. Its aspect ratio χ, defined as the ratio between the length of its major axis and its minor axis, is χ ≈ 2.0 and its terminal Reynolds number is Re ≈ 350. Similarly, a 3 mm air bubble rising in silicone oil T2, whose kinematic viscosity is twice that of water but surface tension is nearly four times smaller [START_REF] Zenit | Path instability of rising spheroidal air bubbles: a shapecontrolled process[END_REF], corresponds to Ga = 90 and Bo ≈ 1.0. Although the value of Ga is close for both bubbles, the terminal Reynolds number of the second one is only 150, owing to its larger oblateness (χ ≈ 3.1). Under such high-Ga low-Bo conditions, an isolated bubble rises in a straight line only if is shape is close enough to a sphere. In contrast, bubbles such that χ 2.2 follow a zigzagging or a (possibly flattened) spiralling path [START_REF] Cano-Lozano | Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability[END_REF]. Then the bubble wake is no longer axisymmetric. Rather, this wake is dominated by a pair of counter-rotating trailing vortices in which the streamwise vorticity is concentrated. A significant part of the paper examines on the role of these trailing vortices on the dynamics of the bubble pair, especially on the direct interaction in some regimes of the trailing bubble with the vortex pair emanating form the leading bubble. The above two examples indicate that the parameter range on which this study focuses corresponds to millimeter-size bubbles rising in low-viscosity liquids, especially water. It must be kept in mind that the zigzagging/spiralling style of rise is specific to such conditions under which bubbles with a low-to moderate Bond number may have a large Galilei number. Bubbles with similar Ga rising in high-viscosity liquids, or centimeter-size bubbles rising in low-viscosity liquids, have large Bond numbers. Such bubbles exhibit a skirted or spherical cap shape and rise in a straight line, even if the corresponding Reynolds number is large [START_REF] Bhaga | Bubbles in viscous liquids: shapes, wakes and velocities[END_REF]. This makes the interaction scenarios for such bubbles distinct from those considered here.

The results to be discussed below are obtained by solving the three-dimensional timedependent Navier-Stokes equations valid throughout the flow (i.e. including capillary effects) with the open source flow solver Basilisk (http://basilisk.fr) described in [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]. The numerical schemes employed in this solver are detailed in the above two references and are summarized in ZNM. In particular, a geometrical volume of fluid approach is employed to track and advance the liquid-gas interface. A sophisticated adaptive mesh refinement technique makes it possible to locally refine the grid close to the interface and within high vorticity regions, a feature that greatly enhances the computational efficiency. In ZNM, an additional grid refinement strategy was activated to capture the flow within the thin gap separating bubbles under near-contact conditions. Such conditions are met here for Bond numbers of O(1). However, since the focus of the present paper is on the role of the leading bubble wake in the dynamics of the tandem, collisions are not examined in detail. That is, the aforementioned thin-film grid refinement strategy is not activated, making us unable to distinguish between the collisions that under real conditions are followed by a bounce of the two bubbles and those leading to their coalescence.

Similar to ZNM, we make use of a cubic numerical domain, with size (240R) 3 . This large size guarantees that artificial confinement effects are negligibly small throughout the considered range of parameters. A free-slip condition is imposed on all four lateral boundaries, while a periodic condition is prescribed on the top and bottom boundaries. Given the importance of air-water bubbly flows in applications, the entire set of computations is carried out with liquid-to-gas density and viscosity ratios set to 10 3 and 10 2 , respectively. The spatial resolution is refined down to ∆ min = R/68 close to the bubble interface and to ∆ = R/17 in the wake. With this discretization, ZNM established that the evolutions of the rising speed of each bubble are grid-independent up to Ga = 30. However, the highest Reynolds numbers considered here are typically four times larger than in ZNM. Therefore the boundary layers may be twice as thin, making it necessary to check whether or not the above resolution remains sufficient under such conditions. To this aim, a grid convergence test was carried out with the set of parameters Ga = 90, Bo = 0.05 for which the Reynolds number of an isolated bubble is close to 470. The results of this test are discussed in Appendix A, from which it can be concluded that the above resolution properly captures the details of the flow throughout the Reynolds number range considered here.

In the course of the paper, we shall often refer to the study of Cano-Lozano et al. (2016) (hereinafter referred to as CL16) who computed the path of isolated bubbles close to the onset of path instability. These computations were carried out with the Gerris open solver. This code is the direct predecessor of Basilisk and makes use of the same algorithms. The grid resolution employed in CL16 was similar to the present one. These remarks are important because they imply that the results of this previous study can be safely used as a reference to compare the evolution of a bubble pair computed here with that of the corresponding isolated bubble reported in CL16. Throughout the paper, all results are normalized using the characteristic length R and time R/g.

Results in the reference case

Overview

Figure 2 summarizes the various interaction scenarios observed in the parameter range covered by the present investigation, 30 < Ga 90, Bo 1.0. The results from ZNM at Ga = 30 are also reported in figure 2(a) to serve as a reference. At this moderate Ga, the Asymmetric Side Escape (ASE) scenario is observed up to Bo ≈ 0.45, beyond which the two bubbles collide and coalesce. At Ga = 40 and 50, the ASE scenario is still present, but only for Bond numbers less than a critical value Bo = Bo o ≈ 0.2. Beyond this threshold, once the trailing bubble (hereinafter abbreviated as TB) has drifted out of the wake of the leading bubble (abbreviated as LB), both bubbles exhibit oscillatory paths. For Bond numbers just beyond Bo o (Bo = 0.3 -0.4 in figure 2(a)), a first oscillatory regime is observed, with the two bubbles zigzagging virtually in the same plane and independently from each other (left panel in figure 2(c)). We call it the Coplanar Independent Zigzagging regime (hereinafter abbreviated as CIZ). At Ga = 40 this regime subsists until a second critical Bond number Bo = Bo c ≈ 0.5 beyond which the bubbles collide. Following the remark in § 3, we stopped the corresponding simulations when the two bubbles were very close to each other and did not attempt to capture the next steps of their dynamics. Increasing the Galilei number to Ga = 50, two new regimes are identified in between the CIZ and collision regimes. First, for Bo = 0.4, the two bubbles still rise independently but their paths stand in two distinct preferential planes, with only tiny excursions of each bubble out of each of them. For this reason, we qualify the corresponding regime as : asymmetric side escape (ASE) regime; •: small-amplitude independent erratic (SIE) regime; : collision regime; , , : CIZ, NIZ and IFS oscillating regimes, respectively. (c): front and bottom views of the paths corresponding to the last three regimes for (Ga, Bo) = (50, 0.3), (70, 0.3) and (90, 0.5), respectively (the red solid and blue dashed lines refer to the LB and TB, respectively). The dash-dotted line in (a) is the neutral curve corresponding to the onset of path instability for an isolated bubble (CL16); the thin solid lines are the iso-M o lines corresponding to different liquids, with the GaInSn liquid metal (Galinstan) at the very bottom (M o = 1.37 × 10 -13 ), then water (M o = 2.54 × 10 -11 at 20 • C), then silicone oils T0 (M o = 1.8 × 10 -10 ), T2 (M o = 1.6 × 10 -8 ) and T5 (M o = 6.2 × 10 -7 ) from bottom to top (see e.g. [START_REF] Zenit | Path instability of rising spheroidal air bubbles: a shapecontrolled process[END_REF] for the characteristics of these oils). Non-coplanar Independent Zigzagging (hereinafter abbreviated as NIZ). Clearly the CIZ-NIZ transition corresponds to a loss of the planar symmetry. Then, increasing the Bond number to 0.5 reveals another type of evolution, which we identify as the Interacting Flattened Spiralling (IFS) regime. Here again, the two bubbles do not rise within the same plane after some time but, as will be made clear later, the motion of the TB continues to be deeply influenced by that of the LB during a significant part of the rise, if not throughout it. The two paths look like very flattened spirals, with the midplane of the TB path slowly rotating in such a way that in several cases the two midplanes are almost orthogonal at the end of the simulation (right panel in figure 2(c)). Last, the two bubbles collide beyond Bo c ≈ 0.7. For higher Galilei numbers (Ga 70), the succession of the NIZ, IFS and collision regimes is observed beyond Bo = 0.1. In contrast, no CIZ regime is encountered for smaller Bond numbers. In that range, instead of performing large-amplitude planar zigzags after the TB has moved out from the LB wake, the two bubbles rise independently from each other with only small-amplitude erratic horizontal excursions. This defines the Small-amplitude Independent Erratic (SIE) regime.

To properly interpret figure 2, it is important to keep in mind what the evolution of an isolated bubble with the same parameters would be. The critical curve corresponding to the onset of path instability for an isolated bubble as determined by CL16 is reported in figure 2(a) (red dash-dotted line). This curve is found to coincide with the transition from the ASE scenario to the CIZ and SIE regimes up to Ga = 70 (the curve is uncertain for Bo 0.05 and Ga 70, as only Bond numbers of O(0.1) or larger were considered in that Ga-range in CL16). This is a clear indication that the CIZ and SIE regimes (and a fortiori the NIZ and IFS regimes), correspond to conditions under which the rectilinear path of each bubble taken alone would be unstable. At this point it may be useful to remind that the path instability of an isolated bubble starting from rest arises through a supercritical Hopf bifurcation, which is the reason why the path subsequently exhibits periodic lateral oscillations (Tchoufag et al. 2014b). Global linear stability analysis allows the spatial structure of the most unstable mode and the associated growth rate and frequency to be predicted at the threshold (the corresponding reduced frequency or Strouhal number, i.e. the frequency normalized by the rise speed, stands typically in the range 0.025 -0.06, increasing with Bo). The style of path (zigzagging vs. spiralling), together with the amplitude of the lateral excursions, may also be predicted close to the threshold through a weakly nonlinear analysis [START_REF] Tchoufag | Weakly Nonlinear Model with Exact Coefficients for the Fluttering and Spiraling Motion of Buoyancy-Driven Bodies[END_REF].

For bubble pairs released in line and such that Bo < Bo c (Ga), the lateral drift of the TB triggers the path instability, giving rise to one of the above five scenarios. The central issue to be examined in the rest of the paper is then to determine whether the subsequent evolution of the two bubbles is similar to that they would follow if rising alone, or whether they instead continue to interact and the observed non-straight paths are still influenced by this interaction. Before discussing the various regimes in detail, it is worth pointing out that not all of them have the same status. Some arise through a true bifurcation, i.e. a change in the symmetries of the system. This is for instance the case of the CIZ and SIE regimes which correspond to the loss of stability of the rectilinear path of an isolated bubble. This is also the case of the NIZ regime which, increasing Bo while keeping Ga fixed, corresponds to the loss of coplanarity of the two paths. In contrast, the definition of the IFS regime is essentially qualitative, this regime being encountered when 'significant' interactions between the two paths subsist in the course of their successive oscillations. Last, the collision regime corresponds merely to a limit case in which the attractive interaction between the two bubbles is strong enough for the gap that separates them to virtually vanish at some point of their ascent. Obviously, the transitions from one regime to another depend in a complex manner on the whole set of control parameters, including the initial conditions of the system. As will be made apparent in § 4, small changes in these conditions may have a deep influence on the development of the interaction process. This is why in certain regions of the (Ga, Bo)-plane, the regime map of figure 2 may change dramatically according to the initial conditions, even going in some cases as far as the disappearance of one of the above regimes.

CIZ regime

Although this regime is only encountered within a narrow range of parameters, it is the one with the simplest characteristics. This is why we discuss it first, based on the parameter set (Ga, Bo) = (50, 0.3). An isolated bubble with these parameters follows a planar zigzagging path and its terminal Reynolds number and aspect ratio are Re ≈ 125 and χ ≈ 2.1, respectively (CL16). The evolution of the corresponding paths and streamlines for the bubble pair are shown in figure 3. Not surprisingly, the front view of the paths (panel (a)) reveals that the initial in-line configuration breaks down in the very asymmetric manner typical of the ASE scenario (12 t 20). Then the TB immediately starts to follow a large-amplitude planar zigzagging path. In contrast, the LB slowly transitions to a similar path and reaches a nearly-saturated stage only near the end of the computation. In this 'final' stage, the two oscillatory paths have almost identical frequencies and amplitudes, and these quantities are similar to those found with an isolated bubble (green line). These features suggest that the two bubbles rise independently of each other. The bottom view of the paths (panel (b)) shows that, up to tiny deviations, the entire sequence takes place within a single vertical plane. This plane is that defined by the initial ASE interaction, i.e. the vertical path of the LB and the inclined path of the TB during its lateral escape at early times (t 16). This lateral drift provides a finite-amplitude disturbance that triggers the path instability of both bubbles. However the corresponding asymmetry is much weaker in the neighbouring of the LB than around the TB, which results in a much slower development of the oscillatory motion of the former. This difference in the magnitude of the flow asymmetry about the two bubbles may be appreciated by examining the distortion of the streamlines and the inclination of the bubble minor axis at t ≈ 18-19 (panel (c)). Remarkably, under present conditions, the path instability mechanism conserves the planar symmetry dictated by this early evolution of the system. As we shall see later, this is an exception rather than a general rule. The evolution of some characteristics of the motion is displayed in figure 4. The early evolution of the rising speed (panel (a), left axis), is similar to that observed in the ASE scenario detailed in ZNM. That is, once the TB enters the wake of the LB, it is strongly accelerated by the corresponding 'sheltering' effect, while the rising speed of the LB only experiences a slight increase resulting from the reduction of velocity gradients in its wake caused by the flow at the front of the TB (4 t 13). As a result, the vertical separation decreases sharply (panel (c)) and the aspect ratio of the TB (panel (a), right axis) drops dramatically because of the suction induced by the low pressure at the back of the LB. Then, for the reasons discussed in ZNM, the in-line arrangement becomes unstable and the TB starts to drift laterally, as the records of the transverse separation S r (panel (c)) and transverse velocity (panel (b)) indicate. Owing to this lateral drift, part of the potential energy of the TB is used to 'feed' its horizontal motion, making its rising speed, V T B , drop dramatically and become even lower than that of the LB, V LB , which in turn results in a re-increase of the vertical separation beyond t ≈ 15.

At t ≈ 19, the TB has left the wake of the LB and starts zigzagging in a vertical plane. The first period of this oscillating motion is still influenced by the large-amplitude shows that the oscillations of its lateral velocity, V r LB , still have a slightly smaller amplitude than those of V r T B at t ≈ 160. Because of this small difference, a slightly larger fraction of the potential energy of the LB is converted into the kinetic energy associated with its rise, making V LB be still slightly larger than V T B , as figure 4(a) confirms. This is the origin of the gradual increase of the vertical separation still present at the end of the sequence in figure 4(c). The final geometry of the tandem may be anticipated from the late stages displayed in figures 4(c) -(d): the average transverse separation will be close to 5.5, a distance at which the interaction between the two bubbles is extremely weak at the current Reynolds number according to the predictions of [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF] for spherical bubbles, and the average inclination of the line of centres will be approximately 24 • . Due to the slight increase in the frequency of the LB path oscillations as they grow, the phase shift between the two paths has decreased continuously during the transient (figure 4(b)). However, there is a priori no reason for this phase shift to vanish eventually, given the quasi-independence of the two paths. Consequently, both S r and θ will continue to oscillate slightly about their mean value.

Another perspective into the dynamics of the tandem is provided by the evolution of the streamwise vorticity distribution in the wake of the two bubbles. Indeed, a spheroidal bubble rising in a straight vertical line having an axisymmetric wake, vorticity is purely azimuthal in this case and has no streamwise, i.e. vertical, component. Successive snapshots of a selected iso-contour of the vertical vorticity, ω y = ±0.5, are plotted in figure 5. Two vortex threads in which ω y periodically changes sign emerge in the TB wake right after its initial lateral drift. In contrast, no structure corresponding to the selected |ω y |-level is detected in the LB wake until t ≈ 35. This finding confirms that only the TB performs large-amplitude lateral oscillations at earlier times, since such oscillations directly result from the lift force originating in the double-threaded streamwise vortices.

Then the size of the vortex threads past the LB increases until t ≈ 140, beyond which the two pairs of threads keep a similar size, as could be anticipated from figure 4(b). It may be noted that, throughout their nearly-parallel rise, the separation between the two bubbles is large enough to avoid each of them to hit the pair of threads emitted by the other. This confirms that in the present regime only minimal interactions subsist between the two bubbles after the initial ASE stage.

NIZ regime

We now increase the ratio of inertial to viscous effects by examining the results obtained with (Ga, Bo) = (70, 0.3) which, according to the phase diagram of figure 2, belongs to the NIZ regime. Under such conditions, an isolated bubble follows a flattened spiralling motion, and its terminal rising speed and aspect ratios are Re ≈ 168 and χ ≈ 2.2, respectively (CL16). The major two differences in the evolution of this bubble pair compared to the previous one may be appreciated from figures 6(a) -(b). First, both bubbles now follow large-amplitude oscillatory paths soon after the initial ASE interaction is completed. Second, these oscillatory motions take place in two distinct vertical planes. These two figures also show that both paths have frequencies and amplitudes very similar to those of the corresponding isolated bubble, which again indicates that the two bubbles rise almost independently beyond the initial ASE stage.

Figure 6(c) displays the evolution of the bubbles shape and orientation, together with that of the streamlines pattern. The comparison with figure 3(c) in the early stage, for instance at t = 12, is enlightening. Clearly, the front part of the LB is significantly flatter in figure 6(c), and the rear part is slightly more rounded. Also, the standing eddy at the back of the same bubble is much larger. That an increase in Ga (or Re) increases the fore- aft asymmetry of an isolated rising bubble in the above way is a well-documented effect [START_REF] Ryskin | Numerical solution of free-boundary problems in fluid mechanics. Part 2. Buoyancy-driven motion of a gas bubble through a quiescent liquid[END_REF][START_REF] Zenit | Path instability of rising spheroidal air bubbles: a shapecontrolled process[END_REF]. It is directly related to the influence of the shear-free condition on the pressure distribution in the liquid at the interface, i.e. to the presence of a boundary layer around the bubble. Since the flattening of the front is stronger than the rounding of the rear, increasing Ga makes the bubble aspect ratio increase. This in turn has a direct influence on the magnitude of the azimuthal vorticity generated at the gas-liquid interface, since the maximum of this surface vorticity (reached near the buble equator) is known to vary as χ 8/3 for large χ [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF]). This increase in the amount of vorticity produced at the surface of the LB is the direct cause of the larger size of the standing eddy noticed above. It is also the cause of the vigorous path instability revealed by the evolution of the LB path. Indeed, it has been established that wake instability past a perfectly spheroidal bubble sets in when the aspect ratio exceeds a threshold χ c ≈ 2.2 [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF]. As the record in figure 7(a) indicates, with (Ga, Bo) = (70, 0.3) the aspect ratio of the LB stays beyond 2.25 from t ≈ 5 to t ≈ 30. Therefore, the conditions required for the wake (hence the path) of this bubble to become unstable are fulfilled, which explains why it has already performed the first half of a large-amplitude zigzag at t = 30. In comparison, with (Ga, Bo) = (50, 0.3), figure 4(a) indicates that the aspect ratio of the LB exceeds 2.2 only during a short transient near t = 15. This is why the lateral excursions of this bubble grow much more slowly. 7 confirm that the dynamics of the tandem have reached a fullydeveloped state soon after t = 30. Indeed, the vertical separation stabilizes at S ≈ 9.7 and no longer varies. Owing to the phase shift between the two paths revealed by figure 6(a), the horizontal separation S r experiences periodic large-amplitude oscillations with a reduced frequency St ≈ 0.031. These oscillations make S r vary from 2.25 to 10, which in turn results in an inclination of the tandem varying from 13 • to 45 • . As figure 7(a) shows, the rising speed and aspect ratio of each bubble experience significant oscillations during this fully-developed stage. The frequency of these oscillations is twice that of the other quantities, the rising speed reaching its maximum twice during a complete zigzag. The oscillations of the aspect ratio are enslaved to those of the rising speed, being driven by the instantaneous value of the Weber number. A noticeable feature of the transition from the initial transient to the fully-developed state is the significant decrease in the rising speed of the LB (similar to that of the isolated bubble), from V LB ≈ 2.65 for t 28 to an average value close to 2.45 at later times. As discussed in CL16, this reduction results on the one hand from the increased dissipation in the wake associated with the double-threaded trailing vortices, and on the other hand from the part of the bubble potential energy spent to 'feed' its lateral excursions rather than its rise.

As pointed out before, the bottom view of the path (figure 6(b)) reveals that the LB oscillates in a plane which does not coincide with that selected by the TB during its initial lateral drift, plane A say. An early indication of this angular splitting is provided by the two perpendicular views at t = 17 in figure 6(c). Indeed, the flow at the back of the LB is symmetrical neither in plane A nor in the perpendicular vertical plane, plane B say. Consequently, transverse forces move the LB out of these two planes and force it to oscillate in a vertical plane having an intermediate orientation. Given the above discussion, the reason for this behaviour is obvious: the oblateness of the LB being beyond the wake instability threshold χ c , non-axisymmetric disturbances round this bubble develop independently from those generated by the TB drift. Because of this, the two bubbles do not select the same plane of rise and the system does not preserve a planar symmetry.

SIE regime

For Ga 70 and Bo 0.1, a specific regime takes place. Here, once the TB has drifted out of the LB wake, the two bubbles do not follow large-amplitude zigzagging or spiralling paths, nor do they rise strictly in a straight line. Figure 8 shows a typical example of this regime, corresponding to Ga = 90 and Bo = 0.05. As panel (a) reveals, the LB almost rises vertically throughout its ascent, with however some tiny meandering. On the other hand, after having started its lateral drift very early through the usual ASE mechanism, the TB continues to drift gradually well after it has moved out of the LB wake. In figure 8(b) it is noticed that the TB keeps on rising within a vertical plane during some time (t 17), before its path becomes three-dimensional. Then the horizontal excursions of this bubble exhibit erratic orientations, as do those of the LB from the very beginning of its ascent. These horizontal motions remain of small amplitude, typically some tenths of the bubble radius, to be compared to 2 -4 radii in figures 3 and 6. Clearly, the two bubbles no longer interact beyond the initial ASE stage. The above characteristics are shared by all bubble pairs standing in the SIE regime.

A remarkable feature in the streamline patterns displayed in figure 8(c) is that no standing eddy exists past the two bubbles. Their final aspect ratio is close to 1.75 (figures 2(b) and 9(a)), an oblateness for which the numerical results of [START_REF] Blanco | The structure of the axisymmetric high-Reynolds number flow around an ellipsoidal bubble of fixed shape[END_REF] for an isolated bubble indicate that a standing eddy exists only if the Reynolds number is less than 160. Here the final Reynolds number of the two bubbles is almost three times larger (Re ≈ 470, see figure 2(b)), implying a significantly lower accumulation of vorticity at the back of the bubble. Consequently, the absence of standing eddies in figure 8(c) is no surprise. What may sound more surprising is that the system exhibits path instability (albeit with weak manifestations) while the wake instability past a fixed spheroidal bubble only takes place if the aspect ratio is beyond χ c ≈ 2.2 [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF]. Similar small-amplitude chaotic paths were identified in CL16 in the same (Ga, Bo) range as the present SIE regime. As discussed in that reference, the reason why the path of a bubble may be unstable while its wake is still intrinsically stable is that freely-moving bubbles must satisfy overall constraints related to the zero-torque and constant-force conditions. These constraints are responsible for the existence of specific instability modes, which may under certain circumstances become unstable before (i.e. at lower χ) than those associated with the wake instability. Possible shape oscillations play no role in this phenomenon, as it also exists (and is better documented) for rigid bodies. For instance, numerical simulation [START_REF] Auguste | Falling styles of disks[END_REF]) and linear stability analysis (Tchoufag et al. 2014a) reveal that, for certain solid-to-fluid density ratios, path instability of freely-falling disks occurs at a critical Reynolds number more than three times smaller than the threshold of wake instability past a disk held fixed at normal incidence in a uniform stream. Figure 9 helps understand the reason and consequences of the unexpectedly long drift of the TB. In particular, panel (a) reveals that, while the rising speed of the LB has reached its final value at t ≈ 13, its aspect ratio still increases significantly until t ≈ 20. This is an indication that the bubbles have not yet reached their final shape during the stage when they interact. The reason for this delay is that the rate at which the bubble shape adapts to the dynamical conditions is limited by viscous diffusion. As the corresponding characteristic time is very long under such highly inertial conditions (the viscous time scale being proportional to Re), the aspect ratio χ(t) at a given time is lower than what the quasi-static approximation of [START_REF] Moore | The velocity of rise of distorted gas bubbles in a liquid of small viscosity[END_REF] based on the Weber number W e(t) built with the rising speed V (t) would predict. A similar argument applies to the evolution of the TB wake. During the time the TB drifts across the flow region disturbed by the presence of the LB (say t 10 according to figure 8(c)), its wake exhibits the classical double-threaded structure resulting from the ambient shear created by this disturbance (see the snapshots at t = 8 and 10). Once the TB has reached the flow region left undisturbed by the LB, its wake tends to recover its axial symmetry since it is intrinsically stable. However the return to this symmetry is governed by the disappearance of the trailing vortices and the reorientation of the bubble whose minor axis must realign with the vertical. Viscous mechanisms playing a central role in both processes, the return to an axisymmetric flow structure past the TB is far from instantaneous, as the slow decay of the transverse velocity in figure 9(b) and the nonzero inclination of the bubble equator . In anticipation to the regime discussed below with the same Ga but a ten times larger Bo (see figure 12(c)), one may also notice in figure 9(c) that the vertical separation does not drop during the initial period when the TB stands in the LB wake. Indeed, such low-Bo bubbles being only moderately deformed and their rise Reynolds number being very large even at this early stage, their wakes remain thin, making the attractive influence of the LB too weak to reduce significantly S below its initial value.

IFS regime

We now focus on the set of parameters (Ga, Bo) = (90, 0.5). With these parameters, an isolated bubble follows a flattened spiralling path (CL16), with final characteristics Re ≈ 175 and χ ≈ 2.6. Figure 10 shows some streamlines past the two bubbles during their initial in-line rise. Owing to the large curvature of the LB in the equatorial region, a large standing eddy develops behind it. Asymmetries in the streamlines pattern become visible at t ≈ 13. At this moment, the standing eddy at the back of the LB and the disturbance flow past the TB are directly connected, since the closure zone of the former almost coincides with the upper end of the closed streamlines ahead of the TB. Therefore, no intermediate region with a nearly-parallel flow subsists in between the two bubbles, unlike the situation prevailing at earlier times. For this reason, the two bubbles are tightly coupled and any tilting of one of them tilts the other in the same direction, as the last two snapshots in the series confirm. This is why they both start to drift in the same direction. Moreover, the flow asymmetry associated with this lateral motion makes both of them take a pronounced egg-like shape, with the pointed end directed toward the direction of the drift. That the two bubbles initially drift in the same direction is strikingly different from what was observed in the previous three regimes, which all involve an initial ASEtype stage.

The bottom view of the paths at larger times is displayed in figure 11, the initial (red) stage corresponding to the time interval during which the two bubbles drift in the same direction. The path of the corresponding isolated bubble is also shown to serve as a reference (green trace in the bottom part of the figure). It reveals that the bubble follows an almost perfect planar zigzagging motion. In the case of the bubble pair, after a transitional (blue) stage, say beyond t 35, the LB first follows a nearly planar zigzagging path (orange stage), while the TB starts to describe a flattening spiral with a 5 : 1 aspect ratio. Beyond this first complete oscillation, the two bubbles follow significantly different evolutions. The initial planar zigzagging path of the LB turns gradually into a flattened spiralling motion whose aspect ratio decreases over time and is close to 4.5 : 1 at the end of the simulation. Meanwhile, the corresponding path experiences a weak clockwise precession, the mid-plane of the oscillations having rotated by approximately 10 • over the five complete oscillations covered by the simulation. Conversely, the path oscillations of the TB become more and more two-dimensional: while the first (orange) flattened spiral beyond t = 35 has a 6.5 : 1 aspect ratio, the last (black) oscillation is virtually a planar zigzag. In the meantime, the mid-plane of the oscillations has rotated clockwise by approximately 60 • , revealing a precession six times faster than that of the LB. As a result, the two horizontal traces of the paths which were initially nearly aligned make an angle close to 80 • at the end of the simulation.

Figure 12(a) shows how the rising speed and aspect ratio of the two bubbles vary along their rise. Both quantities experience strong and rapid variations during the stage 15 t 25 which corresponds to the setting up of the zigzagging/spiralling motions. Beyond this point, they reach a 'developed' state characterized by well-defined mean values on which quasi-periodic oscillations with a dominant Strouhal number St ≈ 0.042 superimpose. These oscillations are far from sinusoidal and a series of higher harmonics with a significant amplitude is present. This is especially visible in the evolution of the aspect ratio which exhibits numerous spikes, reflecting the fact that the strong and asymmetric deformation of both bubbles involves a combination of modes. It is worth noting that the average rising speed of each bubble is close to 1.9, to be compared with values close to 2.5 in the CIZ and NIZ regimes (see figures 4 and 7), or even 5.3 in the case of the SIE regime (figure 9). This is a direct consequence of the larger bubble deformation which implies stronger velocity gradients in the near-bubble flow, from which a larger drag results. Quite remarkably, the transverse velocities seen in figure 12(b) exhibit very similar evolutions, both in amplitude and frequency, despite the significant differences noticed in the record of the lateral motions. Differences are however visible during short specific stages, most notably near the extrema.

The records of the horizontal and vertical separations in figure 12(c) reveal interesting features. In particular, apart from a small bump during the stage 23 t 33, the vertical separation continuously decreases over time, finally stabilizing around S ≈ 2.6. As a result, beyond t ≈ 70, the inclination of the line of centres (not shown) stabilizes around θ ≈ 65 • with ±5 • oscillations. This small final vertical separation is to be compared with the final values S ≈ 10 or even larger found in figures 4, 7 and 9. Most of the decrease of S takes place before the horizontal oscillations fully develop. The reason for this is easily identified from figures 10 and 11. Since the two bubbles drift laterally in the same direction during a significant time, the 'shielding' effect of the LB lasts for a much longer time than in the previous examples, allowing the TB to get much closer to the LB before the zigzagging/spiralling oscillations fully develop. This makes the vertical separation decrease continuously up to t ≈ 23, when it reaches a first minimum close to 3.2. The evolution of the horizontal separation is qualitatively similar to those observed in the CIZ and NIZ regimes. In particular the average separation is again close to 6 at the end of the sequence. In contrast, the amplitude of the S r -oscillations is nearly three times less than in figure 7(c), a result of the gradual reduction of the lateral excursions of the LB. Both the average value and the oscillations of S r are still growing at the end of the sequence, confirming that the two bubbles keep on following different horizontal dynamics.

Figure 13 details the path and wake structure of the two bubbles during three selected stages of their rise. In each of them, the tandem first drifts toward the right with the TB ahead of the LB, as shown in the inset of panel (a). Then the direction of the drift reverses and the LB takes the lead. In the first half of each stage, say 21.5 < t < 28 in (a), the LB wake hardly alters the motion of the TB which drifts ahead of it. In contrast, in the second half, say 28 < t < 34.5 in (a), the double-threaded vortex structure shed by the LB hits the TB, especially through the negative (green) thread, and thereby influences its path. In (b) this interaction is responsible for the smaller radius of curvature of the TB path observed during the time interval 40.4 < t < 46.7, a feature that initiates the clockwise precession characterizing the subsequent evolution of this path. This scenario repeats itself and the clockwise precession of the TB path increases every time the negative vortex thread hits the TB. However, as figure 11 shows, this precession gradually makes the paths of the two bubble develop in orthogonal planes. This tends to keep the transverse distance between the two bubbles larger during their ascent, thus reducing the opportunities for the LB wake to directly influence the TB. The reason why the interaction process described above results in a precession of the TB is easily understood using figure 14 and fundamental laws of vortex dynamics. In this figure, both bubbles are drifting toward the left and the streamwise vorticity is negative (resp. positive) in the upper (resp. lower) thread of both pairs (panel (a)). This corresponds to an antagonistic interaction configuration in which the induced velocity field tends to make the two pairs repel each other. For this reason, the vortex pair 1 (corresponding to the TB) is deviated in the anticlockwise direction, making the TB rotate clockwise since the force acting on the bubble is opposite to that acting on the fluid (panels (b) -(d)). The same reasoning indicates that the LB rotates anticlockwise. Another example of interacting paths is displayed in figure 15, the control parameters being now (Ga, Bo) = (50, 0.5). Similar to the previous case, the two bubbles start to drift in the same direction, with the TB well ahead of the LB (panel (b)), implying little influence of the wake of the latter on the TB in this early stage. Then at t = 23, the horizontal motion of the TB changes sign, which somewhat later (t = 25) brings the two bubbles back to the in-line configuration, implying S r ≈ 0 and θ ≈ 0 (panels (d) and (e)). During this second stage, the dynamics of the TB is strongly influenced by the LB wake, the streamwise vortices shed downstream hitting directly the TB as the iso-contours of the vertical vorticity in figure 15(b) indicate. This intense interaction first manifests itself in the sharp decrease (resp. increase) of the TB aspect ratio (resp. rising speed), as the dashed curves in the same panel confirm. Then a violent deflection of the TB path toward the X < 0 direction takes place at t = 27, initiating a spiralling anticlockwise motion. This is also a consequence of the above interaction: as a slight misalignment of the two vertical planes in which the paths take place has developed gradually (see the successive positions of the blue and red circles beyond t = 23 in panel (a)), the impulse transferred to the TB by the streamwise vortices has a nonzero transverse component with respect to its ongoing path, from which the observed lateral deflection follows. The reaction of the LB takes place somewhat later (t = 30), resulting in a similar and even more abrupt deflection of its path in the opposite direction (X > 0). This intense sequence is succeeded by a long transitional stage (30 t 45) during which the two bubbles follow somewhat erratic paths and the intensity of their interaction decays, owing to the gradual increase of their horizontal separation. At t = 45, the tandem is close to the side-by-side configuration (panel (e)). Then, each path evolves almost independently in a much more conventional way. On the one hand, the TB follows a largeamplitude planar zigzagging path with a crest-to-crest amplitude of 3 -4 bubble radii.

On the other hand, the LB describes a flattened spiral with a major axis 2 -2.5 bubble radii long. Interestingly, the bubble rotation switches from anti-clockwise to clockwise at some point. The lateral excursions of the TB being larger than those of the LB, the horizontal velocities follow the same trend. This leaves a smaller fraction of the potential energy available for the rise of the TB, making the rising speed of the LB slightly larger, which in turn tends to increase the vertical separation as panel (d) confirms. In summary, the dynamics examined here represents an intermediate case. The two bubbles strongly interact for t 40, a first 'life' during which the system stands in the IFS regime. Then, they evolve almost independently with the average direction of their oscillations lying in two different vertical planes. This makes their behaviour similar to that observed in the NIZ regime, except that here the LB follows a very flattened spiralling path instead of a strictly planar zigzagging path.

Influence of some parameters

Influence of the Bond number

The influence of the Bond number on the fate of the bubble pair was discussed in detail in ZNM for Ga 30. It was shown that the direct connection between the bubble oblateness and the magnitude of the vorticity generated at the gas-liquid interface results in a strong increase of the wake 'sheltering' effect with the Bond number. For this reason, the larger Bo is, the stronger the attraction of the TB toward the LB is and the later its lateral escape starts. This is what makes the two bubbles eventually collide and coalesce beyond the critical Bond number Bo c introduced earlier. A side effect is that, for Bo < Bo c , the final horizontal separation of the two bubbles decreases as Bo increases, given the shorter time during which the TB is able to drift. The above mechanism still holds in the parameter range considered here. However, in this higher-Ga range, the fact that vortex shedding takes place makes the horizontal separation reached after the initial drift of the TB a crucial parameter for the possibility of further interactions. In other words, increasing Bo (still with Bo < Bo c ) makes the system shift from the CIZ or NIZ Figure 17. Variations of the path geometry (bottom view) with the initial angular deviation θ0 for three typical bubble pairs which, in the reference case θ0 = 0 • , evolve in the CIZ, NIZ and IFS regimes (from top to bottom). Red, blue and green lines correspond to θ0 = 0 • , 1 • and 2 • , respectively. A shift in the initial X-position has been applied to each pair for readability; the numbers along the green dashed line in the central row indicate the corresponding time instant. Initially, the two bubbles stand a dimensionless horizontal distance Sr = S0 tan θ0 apart, with the TB shifted in the X > 0 direction. regimes to the IFS regime.

Figure 16, which displays the horizontal trace of the paths, confirms the expected tendency. For Ga = 50, the gradual reduction of the horizontal separation as Bo increases is clearly seen in the ASE regime (Bo 0.2). Then, the dramatic change of the trace from Bo = 0.3 to 0.5 helps appreciate how the mechanisms involved in the interaction of the two bubbles modify the geometry of their paths during the succession of the CIZ, NIZ and IFS regimes. The change of the horizontal trace through the NIZ-IFS transition that takes place in the range 0.4 < Bo < 0.5 is also well visible for Ga = 90. At this high Ga, bubble pairs with Bo 0.1 rise in the SIE regime. This is why in figure 16 the two bubbles of the corresponding pairs are seen to perform only tiny horizontal displacements once the initial ASE drift is completed.

Influence of an initial angular deviation

In ZNM it was shown that for Ga 30 a small nonzero initial inclination θ 0 between the two bubbles has a dramatic influence on their fate. The initial configuration being already non-axisymmetric in this case, the lateral drift of the TB starts much earlier, i.e. at significantly larger separations, than in the θ 0 = 0 • case. As a consequence, the DKT interaction which characterizes the dynamics for Ga = O(10) and Bo 0.2 when θ 0 = 0 • no longer exists with an initial inclination as weak as 2 • . The critical Bond number beyond which the two bubbles collide and eventually coalesce also increases significantly with θ 0 , owing to the longer time offered to the TB to leave the LB wake. Last, when the two bubbles do not collide, their final horizontal separation is smaller than in the θ 0 = 0 • case, since the shear that drives the lateral drift of the TB is weaker, owing to the larger S at which this drift takes place. Here we question the influence of θ 0 on the dynamics of the bubble pair at higher Ga. As an introduction, figure 17 shows how the horizontal trace of the paths varies when the initial inclination is increased up to 2 • . Clearly this parameter has a dramatic influence on the three-dimensionality of the path, tending to constrain it to develop in (or close to) the vertical plane containing initially the two bubbles, here the (X, Y ) plane. The reason for this is clear. For an isolated bubble, path instability arises through a bifurcation that preserves a symmetry plane in the wake [START_REF] Mougin | Path instability of a rising bubble[END_REF]Tchoufag et al. 2014b). The orientation of this plane is arbitrary and is usually dictated by some initial disturbance. Here, due to the nonzero θ 0 , the flow past each bubble does not strictly preserve an axial symmetry about the vertical direction, even at short time. For this reason, the bifurcation leading to non-straight paths is imperfect, the (X, Y ) plane providing a preferential orientation to the symmetry plane. Therefore, the zigzagging motion of the TB resulting from this imperfect bifurcation takes place in this preferential plane, until the LB wake possibly provides a subsequent disturbance having a component out of that plane. Figure 17 confirms that this scenario holds in all cases. Because of the presence of a preferential vertical plane, configurations such as that found for Ga = 90, Bo = 0.5 with θ 0 = 0 • , for which the two paths evolve in a complicated way and eventually stabilize in nearly perpendicular planes, no longer exist. Instead, for the same Ga and Bo but θ 0 = 2 • , the influence of the initial angular deviation is seen to be sufficient to make both bubbles perform large-amplitude nearly planar zigzags within vertical planes whose orientation almost coincides with that of the (X, Y ) plane. A major consequence of this modified geometry is that, beyond the initial ASE stage, the two bubbles maintain a sufficient separation for their interaction to be very weak. In other terms, while the bubble pair with Ga = 90, Bo = 0.5 belongs to the IFS regime when θ 0 = 0 • , it rather stands in the CIZ regime (or the weakly three-dimensional NIZ regime) when θ 0 = 2 • . We observed the same qualitative change for all bubble pairs standing in the IFS regime in the strictly vertical initial configuration. In other words, in the (Ga, Bo) range we explored, we found that the IFS regime does not subsist when θ 0 = 2 • . Conversely, as we shall see below, introducing an initial inclination may instead promote interactions throughout the ascent of the bubbles in the other regimes, even though the two planes of rise coincide or make a small angle.

To illustrate the influence of the initial deviation in more detail, we consider the bubble pair with (Ga, Bo) = (50, 0.3). As figure 17 revealed, any nonzero initial deviation makes the corresponding tandem rise within the (X, Y ) plane. The earlier start of the TB drift triggered by the slightly asymmetric initial configuration is clearly seen in figures 18(a) and (c). This earlier start has a direct consequence on the future rise of the LB: since S is larger than in the reference case during the TB drift (S ≈ 8 instead of S ≈ 5 -6 for θ 0 = 0 • according to figure 18(d)), the disturbance induced near the LB is smaller, delaying the triggering of its path instability and reducing its growth rate (compare the blue and red solid curves in panels (a) and (c)). The weaker shear rates encountered at this larger S by the TB during its drift across the LB wake drastically reduce the lateral extension of the paths in the X-direction, i.e. the maxima of S r , confirming the indication of figure 17. Indeed, the mean value of S r which is close to 4.5 when θ 0 = 0 • is reduced to 2.1 for both θ 0 = 1 • and θ 0 = 2 • (figure 18(d)). The large-amplitude zigzags performed by the TB and their phase difference with the slowly growing oscillations of the LB make S r experience large variations, with crest-to-crest amplitudes as large as 2.5 in the late stage of the simulations. The fact that the minima of the horizontal separation periodically reach values in the range 1 S r 2 while the vertical separation is large (S 9) implies that the TB repeatedly comes back into the LB wake. Because of this, the dynamics of the TB remain influenced by its interaction with that wake throughout its rise. Actually, this interaction even tends to strengthen as time proceeds, the vertical separation increasing by only 35% from t = 20 to t = 80 while the minimal inclination of the tandem reduces by nearly 60% and becomes less than 10 • for t > 50 (figure 18(e)). This continuous interaction, combined with the large amplitude of the zigzags performed by the TB, makes V T B experience growing periodic variations with a relative magnitude up to 5% in the late stage; in contrast the rising speed of the LB stays remarkably constant (and independent of θ 0 ) beyond t = 20 (figure 18(b)). In summary, the evolution of this bubble pair for θ 0 = 0 • reveals the existence of a regime that was not encountered in the reference case: the system still evolves within a single vertical plane, both bubbles perform large-amplitude zigzagging motions, but the TB continues to experience the influence of the LB wake throughout its rise, far from the 'independent' evolution discussed in § 3.2.

The influence of θ 0 on the evolution of the bubble pair discussed in § 3.3 (Ga = 70, Bo = 0.3) is also worthy of some comments. The characteristics of the corresponding dynamics are detailed in figure 19, the bottom view of the paths being displayed in the central row of figure 17. The general trends of this dynamics are similar to those already identified in the strictly in-line configuration. In short, the two bubbles soon develop large-amplitude planar zigzagging motions and, beyond a certain stage, evolve independently from each other in two distinct vertical planes. This makes this pair belong to the NIZ regime whatever θ 0 in the range 0 • -2 • . What is of interest here is the way the two bubbles interact before this 'independent' state, especially the processes that lead to the selection of the vertical plane in which each of them eventually rises. In the original in-line configuration, the vertical plane in which the LB oscillates is selected by the instability that develops in its near wake (as in the case of an isolated bubble with similar characteristics), while the TB oscillates in a plane whose orientation is dictated by its initial drift. Since the two selection mechanisms are independent, the two planes have different orientations from the very beginning of the LB oscillations. This is no longer the case for θ 0 = 1 • and 2 • . Here, not only is the flow asymmetry introduced by the initial inclination of the bubble pair sufficient to select the direction of the initial TB drift, but it also almost controls the orientation of the symmetry plane resulting from the (now imperfect) bifurcation initiating the zigzagging motion of the LB. Therefore the TB oscillates within the (X, Y ) plane and the LB does the same within a plane making only a small angle with the former. Obviously, this configuration makes the TB more prone to experience disturbances that later emanate from the LB wake. This could already be the case at t ≈ 20: at this instant, the TB has already completed its first zigzag while the LB only starts its own, so that the large phase difference makes the two bubbles almost realign vertically (figures 19(a), (d) -(e)). However, as the first snapshot inserted in panels (b) and (d) reveals, the streamwise vortices in the LB wake are still weak and the vertical separation is large (S ≈ 10). Therefore the TB is barely disturbed by the LB wake and starts a new zigzag, still in the (X, Y ) plane. Then, the oscillations of the LB path saturate quickly (t ≈ 25), so that the streamwise vortices soon develop downstream and reach the TB at t ≈ 31 (second snapshot in panels (b) and (d)), although the two bubbles are widely separated at this stage (S ≈ 12, S r ≈ 5.5). Since the symmetry plane of the LB wake makes a small but nonzero angle with the (X, Y ) plane, the momentum transferred to the TB during the collision with the streamwise vortices has a nonzero component perpendicular to the latter plane. This deflects the path of the TB toward another plane distinct from the previous two. From there on, the possibility for another interaction sequence of the same nature is much reduced, and the system now evolves in a NIZ configuration very similar to that observed for θ 0 = 0 • . 4.3. Influence of the initial separation in the zigzagging/spiralling regimes So far the initial separation between the two bubbles was maintained at S 0 = 8. We now examine how varying S 0 may affect the interaction process in the parameter range where the CIZ and NIZ scenarios were previously encountered. Figure 20 displays several characteristics of the dynamics obtained by decreasing S 0 from 8 to 4 for (Ga, Bo) = (50, 0.3). As the bottom view of the paths in panel (a) indicates, the system keeps its planar symmetry and the average horizontal separation does not change much as S 0 is varied (see panels (b) and (e) for a quantitative confirmation). Therefore, the bubble pair remains in the CIZ regime at least down to S 0 = 4. The minimum reached by the vertical separation during the initial axisymmetric stage decreases significantly as S 0 is reduced, from S ≈ 5 for S 0 = 8 to S ≈ 2.3 for S 0 = 4. Moreover, the smaller S 0 the shorter the time required to reach this minimum. A direct consequence of this shorter initial stage, already observed at lower Ga in ZNM, is that the smaller S 0 is, the earlier the initial lateral drift of the TB starts, as panels (b) and (d) confirm. At the same time, the reduction of the S-minimum implies that the magnitude of the asymmetric disturbance induced by this drift in the vicinity of the LB increases as S 0 is reduced, triggering the lateral motion of the latter. This is why the zigzagging motion of the LB reaches a saturated state much earlier for S 0 5 (t ≈ 20) than for larger initial separations (t 120 for S 0 = 8, see figure 4). The behaviour of the system for S 0 > 8 may be extrapolated by combining the information provided by figure 20(e). First, since the S-minimum reached during the initial stage increases with S 0 , the two bubbles cannot collide during that stage, and the initial drift of the TB out of the LB wake is likely to take place whatever S 0 . Beyond this stage, the long-term value of S increases monotonically with S 0 while the average lateral separation stays almost S 0 -independent. From this monotonic behaviour in the range 4 S 0 8, it may reasonably be inferred that the system remains in the CIZ regime for S > 8, with possible significant changes only in the final vertical separation, hence in the inclination of the tandem. A qualitatively similar conclusion is reached by considering the pair with (Ga, Bo) = (70, 0.3) and varying the initial separation in the range 5 S 0 20 (not shown): the system stays in the NIZ regime whatever S 0 , although its final geometry, especially the vertical separation of the two bubbles and the orientation of their line of centres, varies significantly with the initial separation.

We then examine the changes in the dynamics of the pair with (Ga, Bo) = (50, 0.5) when the initial separation is varied from S 0 = 8 to S 0 = 12. Figure 21(a) indicates that increasing the initial separation makes the horizontal projection of the paths flatter. Qualitatively, the system first transitions from the mixed IFS-NIZ evolution described in § 3.5 for the reference case S 0 = 8 to a clear IFS evolution for S 0 = 10. Then, for S 0 = 12, the two paths are almost coplanar, so that the observed dynamics look like a CIZ evolution, with the difference that the LB follows a very flattened spiralling path rather than a strictly planar zigzagging motion. For both S 0 = 10 and S 0 = 12, the tandem first exhibits an evolution qualitatively similar to that described in § 3.5, with essentially some time shift. For instance, the system realigns at t ≈ 31 -32, instead of t ≈ 25 in the reference case. This realignment leads potentially to a direct interaction between the TB and the trailing vortices shed by the LB. However, this interaction weakens dramatically as S 0 increases as we shall see. This is because the LB stands right at the threshold of path instability, its aspect ratio being very close to 2.2. Under such nearly critical conditions, the disturbance provided by the initial drift of the TB is crucial for triggering the path instability of the LB (see the discussion in § 3.3). Therefore, the larger S 0 is, the longer it takes for the trailing vortices in the LB wake to grow.

For S 0 = 10, the LB wake is significantly but not completely developed by the time the two bubbles realign. Hence, similar to the reference case (figure 15), the TB undergoes a violent lateral deviation after the tandem has realigned vertically at t ≈ 31. However, here the two bubbles start to spiral right away, while for S 0 = 8 they drift erratically during a long transient before starting to spiral. The absence of this transient is a consequence of the weaker lateral momentum transferred by the trailing vortices to the TB in the present case, owing to the aforementioned incomplete development of the LB wake and the larger vertical separation at the corresponding time (S ≈ 4.5 instead of S ≈ 2, according to figure 21(b)). As a result, the horizontal separation between the two bubbles stays much smaller than in the reference case, even at large time (see figure 21(c), where the S rminima are seen to be now about S r ≈ 0.7, instead of S r ≈ 3 for S 0 = 8). Hence, as the last two snapshots in figure 21(d) confirm, the TB keeps on being under the influence of the LB wake at regular time intervals throughout its rise, providing a new example of an IFS evolution.

Beyond t ≈ 32, the evolution observed when S 0 = 12 differs dramatically from the above one. Indeed, as the orientation of the black and purple segments in the left panel of figure 21(e) indicate, the TB deviates only by a small angle from the path it followed in the earlier stage. For the aforementioned reason, the strength of the trailing vortices past the LB is still weak in this case by the time the two bubbles realign (left snapshot in figure 21(e)). Moreover, the vertical separation is significantly larger than with S 0 = 10 (S ≈ 7 instead of S ≈ 4.5). Both features cooperate to limit the influence of the LB wake on the TB during this crucial stage, so that the initial planar symmetry of the system is only marginally altered. At longer times, the TB essentially performs large planar zigzags. However it may occasionally interact with the LB wake which is now fully developed. During such interactions, some lateral momentum is transferred to the TB, inducing some slow precession of its plane of oscillation (see the right snapshot in figure 21(e) and the change in the horizontal trace at t = 56.5).

Influence of the initial separation in the collision regime

Here we consider the parameter set (Ga, Bo) = (50, 1.0) and vary the initial separation to analyze its influence in the collision regime identified for S 0 = 8 in figure 2. Some characteristics of the dynamics observed by increasing the initial separation up to S 0 = 28 are displayed in figure 22. The two bubbles are found to collide for S 0 20, the collision time increasing significantly with S 0 (figures 22(b) -(d)). They stay significantly apart from each other until the end of the simulation for the largest two separations. Nevertheless, as figure 22(d) reveals, the vertical separation decreases consistently over time whatever S 0 , suggesting that the two bubbles would also eventually collide later for these large initial separations. Present observations are reminiscent of those of [START_REF] Stewart | Bubble interaction in low-viscosity liquids[END_REF] and [START_REF] Brücker | Structure and dynamics of the wake of bubbles and its relevance for bubble interaction[END_REF] with large rising bubbles having Bond numbers in the range 1 -10 and Galilei numbers in the range 100 -500. In these experiments, the zigzagging/spiralling bubbles were always found to collide despite the lateral excursions of the TB. Coalescence never happened, presumably because the interfaces were contaminated in these experiments. Instead, the bubbles repelled each other after collision, leading to a self-repeating scenario.

According to figure 22(a), an isolated bubble with the same Ga and Bo (black line) first follows a flattened spiralling path before describing planar zigzags, an evolution consistent with the findings of CL16. Nevertheless, it takes some time for the initial straight vertical path to become unstable, which explains why the two bubbles collide in the head-on configuration for S 0 = 8. Path instability of the isolated bubble sets in at t ≈ 25 and saturates at t ≈ 70 according to figure 22(c). Consistently, all pairs with S 0 > 8 perform oscillatory paths. For S 0 > 16, the saturated path of the LB takes the form of a planar zigzag or (for S 0 = 24) a flattened spiral. Simultaneously, the TB follows a slowly precessing zigzagging path. A noticeable feature is that the amplitude of the horizontal excursions decreases with S 0 , approximately from 2.5 for S 0 = 20 to 1.5 for S 0 = 28, and reaches minima close to zero at regular time intervals (figure 22(d)). These characteristics, together with the aforementioned decrease of S over time, reinforce the view that a collision happens in all cases. In figure 22(b), the records of V LB are seen to collapse on a single curve, indicating that the LB is barely disturbed by the presence of the TB throughout its rise. The only exception is the very late stage before a collision happens in the head-on (S 0 = 8) or oblique (S 0 = 16) configuration, in which case V LB sharply increases, owing to the vertical impulse transferred to the LB by the TB which is about to touch it. A close look at the plots in figures 22(b) -(c) for S 0 = 28 is of interest. Beyond t ≈ 70, the vertical and transverse velocities of both bubbles appear to have reached a saturated state. However the rising speed of the TB and the maxima of its transverse velocity are slightly larger than those of the LB, indicating that some interaction is still present. Indeed, since the horizontal distance separating the two bubbles remains small, the TB keeps on experiencing the slight velocity defect present in the far wake of the LB.

The remaining issue is to understand why, as figure 22(d) shows, the minimum horizontal separation between the two bubbles returns repeatedly to near-zero values, which leads unavoidably to a collision, and presumably to coalescence in a good number of cases as far as the interfaces are uncontaminated. For this it is necessary to consider the characteristics of the corresponding isolated bubble, whose aspect ratio and terminal Reynolds number are χ ≈ 2.55 and Re ≈ 77, respectively. [START_REF] Adoua | Reversal of the lift force on an oblate bubble in a weakly viscous linear shear flow[END_REF] computed the shear-induced lift force acting on an oblate spheroidal bubble held fixed in a linear shear flow over a wide range of conditions. With the above characteristics, their results indicate that, if the relative shear rate (defined as the shear-induced velocity difference at the bubble scale, normalized by the local relative fluid velocity) stands below a Redependent threshold, the lift force is oriented in the direction opposite to that it would have on a spherical or moderately oblate bubble, other things being equal. The reasons for this change of sign were summarized in appendix C of ZNM. In short, beyond a critical Re-dependent oblateness, the sign (but not the strength) of the streamwise vorticity in each of the trailing vortices of the bubble wake is dictated by the vorticity produced at its surface, and no longer by the upstream vorticity. A detailed analysis of the streamwise vorticity balance reveals that, for a given upstream vorticity, this change results in a switch of the sign of the streamwise vorticity in the trailing vortices, hence in a reversal of the shear-induced lift force acting on the bubble. Quantitatively, this lift reversal may be appreciated by writing the dimensional lift force in the classical form F = C L MU × ω ∞ , with M the mass of fluid that can be contained within the bubble volume, ω ∞ the upstream vorticity and U the relative fluid velocity at the position of the bubble centroid. In the inviscid limit (i.e. with no vorticity generated at the bubble surface by the shearfree condition), the lift coefficient on a spherical bubble is C L = 0.5 [START_REF] Auton | The lift force on a spherical body in a rotational flow[END_REF]. Still in this limit, C L is an increasing function of the aspect ratio, with C L ≈ 1.45 for χ = 2.55. These predictions are only marginally altered by finite-Re effects as far the relative shear rate Sr = R||ω ∞ ||/||U|| is such that Sr 0.1. In contrast, for smaller shear rates, the lift force changes sign beyond Re ≈ 37 for this specific oblateness, and an interpolation of the results of [START_REF] Adoua | Reversal of the lift force on an oblate bubble in a weakly viscous linear shear flow[END_REF] predicts C L ≈ -0.6 for Re = 77 and Sr = 0.02. The manner the above mechanism influences the evolution of the bubble pair in the situation of interest here is as follows. When the TB oscillates along a zigzagging path, the streamwise vorticity in each of its trailing vortices changes sign periodically, as it would for an isolated zigzagging bubble [START_REF] Brücker | Structure and dynamics of the wake of bubbles and its relevance for bubble interaction[END_REF][START_REF] Mougin | Wake-induced forces and torques on a zigzagging/spiralling bubble[END_REF]. The upstream vorticity ω ∞ 'felt' by the TB results from the axisymmetric velocity defect in the LB wake, and this upstream (≈ horizontal) vorticity is responsible for an additional amount of streamwise (≈ vertical) vorticity in the TB wake, compared to the isolated configuration. When the vertical separation between the two bubbles is moderate, say typically S 10, the TB stands in the near wake of the LB. There, the relative shear rate is large enough (Sr 0.1) for the above two sources of streamwise vorticity to cooperate and provide a lift force acting to move the TB away from the wake centreline. Nevertheless, at such short separations, the sheltering effect acting along the wake axis is so strong that the TB does not succeed in drifting laterally by a sufficient distance to avoid colliding with the LB. In contrast, for larger S, the above two mechanisms combine in an antagonistic manner, making the streamwise vorticity slightly more intense when the TB is close to the inflection points of the zigzag (i.e. to the centreline of the LB wake where ω ∞ vanishes) than when it is close to its extremities. This may be confirmed from the curves corresponding to S = 20 and S = 28 in figure 22(c): as the transverse velocity V r reaches its maxima at the inflection points of the path, the fact that these maxima are slightly larger for the TB than for the LB implies that the lift force responsible for the generation of V r is somewhat stronger for the former. This decrease of the lift force acting on the TB as it moves away from the centreline of the LB wake limits the amplitude of its lateral excursions. The weaker Sr, i.e. the larger S, the stronger this effect as figure 22(a) confirms. Hence, in the course of the oscillations of the LB path, the number of opportunities for the lateral separation between the two bubbles to fall to near-zero values increases with S (as seen in figure 22(d)), making a collision unavoidable at some point.

Summary and concluding remarks

In this second part of our investigation, we carried out three-dimensional numerical simulations of the flow past a pair of identical bubbles initially released in line in the parameter range 40 Ga 90, 0.02 Bo 1.0. In that range, provided the Bond number exceeds a Ga-dependent threshold Bo o (Ga) decreasing from Bo o ≈ 0.3 for Ga = 40 to Bo o ≈ 0.1 for Ga = 90, an isolated bubble follows a zigzagging or spiralling path. Considering a fixed initial separation between the two bubbles and assuming their line of centres to be initially exactly vertical, the simulations allowed us to build a phase map gathering the encountered flow regimes in the (Ga, Bo) plane. For Bond numbers less than the above threshold and Ga 70, i.e. weakly deformed bubbles and still moderate Reynolds numbers, the system evolves according to the ASE scenario discussed in detail for smaller Ga in ZNM. Beyond a second threshold ranging from Bo c ≈ 0.5 for Ga 40 to Bo c ≈ 0.7 for Ga 50, the two bubbles always collide. We did not attempt to compute the collision process in detail, but analyzed the mechanisms that make the collision unavoidable under such conditions, despite the large lateral excursions of both bubbles. For Bond numbers in between the above two thresholds, the simulations revealed the existence of three markedly different regimes. The first of them arises for Bo Bo o and Ga 50. There, the initial lateral drift of the TB is succeeded by a 'fully developed' stage in which the two bubbles rise independently within the same vertical plane along large-amplitude zigzagging paths, which defines the Coplanar Independent Zigzagging (CIZ) regime. Increasing the Bond number beyond the upper limit of this regime, one first encounters a Non-coplanar Independent Zigzagging (NIZ) regime. Here again, after the initial ASE stage, the two bubbles rise independently along zigzagging or flattened spiralling paths but the two paths stand in distinct vertical planes. The selection of the CIZ or NIZ regime depends on the ability of the LB to develop an unstable path on its own. If the aspect ratio of the LB is below 2.2, its path can become unstable only due to the external disturbance provided by the lateral drift of the TB. This initial condition forces the two paths to take place within the same plane, leading to the CIZ regime. In contrast, when the oblateness of the LB is large enough, its path instability develops independently from the lateral motion of the TB. This makes the two bubbles rise in distinct planes, yielding the NIZ regime. Increasing the Bond number beyond the upper limit of that regime, but still for Bo < Bo c , the system enters the Interacting Flattened Spiralling (IFS) regime. Here again, both bubbles follow zigzagging or flattened spiralling paths. However the attraction provided by the LB wake is strong enough for the TB to reenter this wake one or more times in the course of its ascent. During such sequences, the trailing vortices shed by the LB directly hit the TB, which results in marked successive lateral deflections of its path. In the most inertial regimes considered here (Ga 70), the CIZ regimes was not observed. Instead, for low Bond numbers (Bo 0.1, which at such Ga is typically encountered in water and liquid metals), the bubble pair first follows an ASE scenario, beyond which the two paths perform small-amplitude erratic horizontal motions independently from each other. This defines the Small-amplitude Independent Erratic (SIE) regime, the existence of which results from the fact that the wake of an isolated bubble rising under such high-Ga low-Bo conditions is intrinsically stable but its path is not, the instability being then driven by the overall force and torque constraints imposed by Newton's second law. Increasing the Bond number beyond the upper limit of the SIE regime, the NIZ-IFS-collision sequence observed for Ga = O(50) beyond the CIZ regime is recovered.

We then showed that, even slight, a nonzero misalignment of the two bubbles from the reference in-line configuration favours oscillations within the same vertical plane. Then, depending on the lateral separation the two bubbles may reach, this constraint may promote or reduce the possibility for further interactions between them. Typically, these interactions weaken when the bubble oblateness is large (Bo ≈ 0.5), i.e. the system stands in the IFS regime for θ 0 = 0 • . Indeed, the initial lateral drift of the TB takes place quite close to the LB in this case, so that the shear 'felt' by the TB is strong and allows the minima reached by the lateral separation to be sufficiently large to avoid the TB re-entering the LB wake. For this reason, the IFS regime no longer subsists for θ 0 = 2 • within the (Ga, Bo) range we explored. The reverse happens when the bubble oblateness is moderate (Bo ≈ 0.3), i.e. the system stands in the CIZ regime for θ 0 = 0 • . There, as the initial drift of the TB takes place at quite large distances from the LB, the shear encountered by the TB during this stage is small and yields small minima of the lateral separation in later stages. Under such conditions, the TB returns repeatedly within the LB wake in the course of its oscillations and is directly hit by the trailing vortices every time the tandem is nearly vertical, maintaining the system in close interaction even in the long term. Intermediate scenarios may happen in the NIZ regime, with the geometrical constraint induced by the initial inclination promoting interactions between the two bubbles during some time, until a lateral deviation of the TB under the action of the LB wake makes the two planes of rise distinct and allows the two bubbles to continue their ascent independently.

The role of the distance separating initially the two bubbles was found to depend dramatically on the considered regime. Variations in S 0 only change the geometry of the final arrangement in the CIZ and NIZ regimes, i.e. for Bo 0.4, reducing (resp.

increasing) the inclination of the bubble pair if the initial separation is increased (resp. reduced). Things are more complex in the IFS regime, i.e. for Bo ≈ 0.5 and Ga 50, because the evolution of the system then crucially depends on both the intensity and the orientation of the interaction taking place between the TB and the LB wake when the two bubbles realign vertically for the first time. If this interaction is strong, the lateral deviation undergone by the TB is large enough for the two paths to remain sufficiently far apart that no direct interaction further happens. Increasing S 0 weakens the above lateral deviation, thus maintaining the two paths closer to each other. This of course favours the possibility of new realignments, hence of further interactions. However these interactions weaken since the vertical separation between the two bubbles increases. So, while some increase beyond the reference value S 0 = 8 tends to maintain permanently the system in the IFS regime, a larger increase makes it transition to the NIZ or CIZ regime after some time. For Bond numbers of O(1), the two bubbles collide no matter how large S 0 . With S 0 10, the approach to collision is driven by the usual attraction mechanism resulting from the longitudinal pressure gradient along the LB wake. Although the lift force acting on the TB tends classically to extract it from that wake, this effect is too weak compared to the wake-induced suction to produce a sufficiently large deviation of the TB avoiding collision. For larger S 0 , the TB stands further downstream in the LB wake, so that the longitudinal pressure gradient (hence the attractive effect) is significantly weaker than in the previous case. In this regime, the Reynolds number and the bubble oblateness being large (Re = O(100), χ 2.5) while the transverse shear perceived by the TB is small, the shear-induced lift force tends to push it toward the wake centreline, according to the reversal mechanism identified by [START_REF] Adoua | Reversal of the lift force on an oblate bubble in a weakly viscous linear shear flow[END_REF]. As a consequence, the TB stays within the LB wake even for large S 0 . The residual longitudinal attraction is then sufficient to make collision inevitable again.

The results discussed in this paper reveal a rich and frequently nonintuitive phenomenology. Beyond the role of the control parameters Ga, Bo, θ 0 and S 0 identified in ZNM for Ga 30, the large-amplitude lateral oscillations of zigzagging or spiralling bubbles rising at higher Ga add an extra level of complexity. Retrospectively, it turns out that this is partly due to fact that the growth rate of these oscillations differs generally for the two bubbles, owing to the differences in the initial disturbances that trigger them. Because of this, the phase difference between these oscillations may under certain circumstances bring the orientation of the tandem back to the vertical or close to it at various stages of the ascent, allowing the trailing vortices present in the LB wake to hit the TB, thereby modifying its dynamics.

Present results are in principle of direct use to understand the behaviour of high-Reynolds-number bubbly dilute plumes released from a single injection point, or from an array of injectors provided the distance between two of them is much larger than the bubble radius. However, the findings of § 4 suggest that their direct applicability is hampered by the sensitivity of the system to initial conditions, especially to the initial angular deviation, which was shown to change the evolution of the bubble pair dramatically in some regimes. In particular, the quantitative relevance of the results obtained for θ 0 = 0 • is limited by the experimental uncertainly on θ 0 , which can hardly be reduced to less than 1 • . Regarding nearly homogeneous high-Reynolds-number suspensions, these results clearly represent only a first step toward a satisfactory understanding of the dynamics at stake. Indeed, there is an infinity of possible initial orientations for each bubble pair in such suspensions, and the in-line configuration considered here is only one of them. In particular, bubble pairs released side by side may exhibit deeply different dynamics, as the 'risk' of a head-on collision and the probability that the two bubbles follow zigzagging paths located within the same plane are much smaller in that case. Therefore, the side-by-side configuration requires specific investigations similar to the present one. The role of size variations within the bubble population also needs to be specifically examined. A slight size difference between the two bubbles in a pair induces a difference between their rise speeds, thus modifying the evolution of their vertical separation when they interact. Hence, depending on whether the larger bubble stands behind the second bubble or ahead of it, this difference reinforces or weakens the interaction. How this simple effect modifies the dynamics of a group of bubbles, depending on the relative position of the larger/smaller test bubble with respect to its neighbours, is far from trivial. Last but not least, collisions are frequent in bubbly suspensions, as soon as the gas volume fraction exceeds a few percent. This makes it necessary to determine the fate of colliding bubbles, i.e. whether they bounce or coalesce depending on the flow conditions. This aspect, which was deliberately disregarded here for computational reasons, represents a major technical challenge.

The above 'warnings' suggest that the route toward a detailed understanding of the dynamics of high-Reynolds-number bubbly suspensions is still long and will require massive computational resources, even if simulations are restricted to a few hundreds of bubbles. Nevertheless, the present study identified basic mechanisms that will help to decipher the complex dynamics of such suspensions, provided they are 'reasonably' dilute. Among these mechanisms, one may especially mention the ability or inability of the path instability of the LB to develop on its own depending on the bubble oblateness, the momentum transfer resulting from the collision between the TB and the trailing vortices of the LB most of the times the former re-enters in the wake of the latter, and the consequences of this transfer on the changes of direction of the TB path. These mechanisms, together with those already known at lower Reynolds number (e.g. the attraction of the TB in the LB wake, or the instability of the in-line configuration), survive whatever the initial conditions and the number of bubbles, and the present findings provide useful information to clarify their role under realistic conditions, especially in air-water flows involving millimeter-size bubbles. The inset in (a) provides a precise view of the residual differences between the rising speeds predicted on the two grids at t ≈ 20. The small high-frequency oscillations visible in this inset result from the dynamic changes of the grid near the bubble surface in this low-Bo high-Ga regime.

(a) (b) R/136), approximately three (resp. six) grid cells stand within the boundary layer. If grid convergence is reached, the rising speed of each bubble should be grid-independent, except during the early stage corresponding to the onset of the lateral drift of the TB across the LB wake. Indeed, the initial axial symmetry of the system being broken by an instability, the transition to a non-axisymmetric state is governed by the asymmetry of 'natural' numerical disturbances present in the code. Since these disturbances depend on the grid resolution (among other things), the onset of the lateral motion varies with the spatial resolution and no grid independence may be expected during that stage. A direct consequence of this grid-dependent onset is that the vertical and horizontal separations between the two bubbles cannot be strictly grid-independent on the long term, as they result from the time integration of the bubble velocities (see Appendix A in ZNM for a detailed discussion). The test was carried out on the reference grid with ∆ min = R/68 and on a twice as fine grid with ∆ min = R/136. The computational time required on this refined grid is very large, owing to the larger number of cells and smaller time step. This is why the corresponding run was stopped at t ≈ 23, beyond which the two bubbles rise independently. On both grids, the tolerance applied to the Poisson solver for the pressure field, defined as the maximum relative change of the fluid volume enclosed in a cell over one time step, was fixed to 1 × 10 -4 . The results of this test are summarized in figure 23. The inset in panel (a) reveals that the rising speeds reached at t = 23 on the two grids differ by less than 0.4%, and the difference is still reducing at that time. The bubble aspect ratio being close to 1.75, the contribution of the boundary layer and wake to the drag, hence to the rising speed, is approximately 3.5% at this Reynolds number [START_REF] Moore | The velocity of rise of distorted gas bubbles in a liquid of small viscosity[END_REF]. Therefore it can be concluded that this contribution, while small, is properly captured with ∆ min = R/68, despite the limited number of cells standing within the boundary layer. This establishes the grid convergence of present high-Ga results in the sense defined above. The differences observed on the rising speeds during the lateral drift of the TB and their consequences on the two components of the separation may easily be rationalized. The onset of the TB drift is seen to happen somewhat later on the finer grid. This delay implies that, on that grid, the vertical separation between the two bubbles is smaller when the TB starts drifting, resulting in a smaller vertical separation during some time (see the red curves in figure 23(b) in the time interval 7 t 10). However this shorter S implies that the TB experiences a larger ambient shear during some time, hence a stronger lift force. This is why the growth rate of the horizontal separation is larger on the finer grid. Owing to this stronger transverse force, a larger fraction of the potential energy of the TB is converted into kinetic energy associated with the transverse motion, at the expense of the vertical motion. This is confirmed in figure 23(a), where the minimum rising speed of the TB during its lateral drift is seen to be slightly smaller on the finer grid. Due to this more severe transient reduction of the TB speed, the vertical separation observed on that grid becomes larger beyond t ≈ 11. The difference stops increasing at the end of the TB drift but persists in the long term, making the final S slightly larger on the finer grid.

Figure 1 .

 1 Figure 1. Sketch of the problem. (a) Initial configuration; (b) definition of the geometric parameters. (a) (b)

Figure 2 .

 2 Figure 2. Configuration map based on simulation results: (a) (Ga, Bo)-plane; (b) (Re, χ)-plane.: asymmetric side escape (ASE) regime; •: small-amplitude independent erratic (SIE) regime; : collision regime; , , : CIZ, NIZ and IFS oscillating regimes, respectively. (c): front and bottom views of the paths corresponding to the last three regimes for (Ga, Bo) = (50, 0.3), (70, 0.3) and (90, 0.5), respectively (the red solid and blue dashed lines refer to the LB and TB, respectively). The dash-dotted line in (a) is the neutral curve corresponding to the onset of path instability for an isolated bubble (CL16); the thin solid lines are the iso-M o lines corresponding to different liquids, with the GaInSn liquid metal (Galinstan) at the very bottom (M o = 1.37 × 10 -13 ), then water (M o = 2.54 × 10 -11 at 20 • C), then silicone oils T0 (M o = 1.8 × 10 -10 ), T2 (M o = 1.6 × 10 -8 ) and T5 (M o = 6.2 × 10 -7 ) from bottom to top (see e.g.[START_REF] Zenit | Path instability of rising spheroidal air bubbles: a shapecontrolled process[END_REF] for the characteristics of these oils).

Figure 3 .

 3 Figure 3. Paths of a bubble pair with (Ga, Bo) = (50, 0.3) following a CIZ evolution. (a) Side view, with r the transverse distance to the initial path and numbers indicating the time instant at the corresponding position; (b) bottom view; (c) snapshots of the streamlines past the tandem in the reference frame of the LB. In (a) and (b), the red solid and blue dashed lines refer to the LB and TB, respectively, while the green dash-dotted line corresponds to the path of the corresponding isolated bubble.

Figure 4 .

 4 Figure 4. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (50, 0.3). (a) Rising speed (colored lines, left axis) and aspect ratio (black lines, right axis) of the LB (solid lines) and TB (dashed lines); (b) same for the horizontal velocity component (numbers indicate the corresponding time instant); (c) vertical (red line, left axis) and transverse (blue line, right axis) separations; (d) inclination of the line of centres with respect to the vertical. In (a) -(b) the color code is identical to that in figure 3.

Figure 5 .

 5 Figure 5. Successive snapshots of the iso-surface ωy = -0.5 (green) and ωy = 0.5 (red) during the rise of a bubble pair with (Ga, Bo) = (50, 0.3).

Figure 6 .

 6 Figure 6. Paths of a bubble pair with (Ga, Bo) = (70, 0.3) following a NIZ evolution. (a) Side view, with r the transverse distance to the initial path and numbers indicating the time instant at the corresponding position; (b) bottom view; (c) snapshots of the three-dimensional streamlines past the two bubbles in the reference frame of the LB. At t = 17, streamlines are shown in two perpendicular vertical planes, plane A containing both the initial vertical path and that followed by the TB during its first lateral drift. In (a) -(b), the color code is similar to that in figure 3; the side view in (a) was obtained by considering the path of each bubble in its own plane and juxtaposing the three paths in the same plane.

Figure 7 .

 7 Figure 7. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (70, 0.3). (a) Rising speed (colored lines, left axis) and aspect ratio (black lines, right axis) of the LB (solid lines), TB (dashed lines) and isolated bubble (dash-dotted lines); (b) same for the horizontal velocity component; (c) vertical (red line, left axis) and transverse (blue line, right axis) separations; (d) inclination of the line of centres with respect to the vertical.

Figure 8 .

 8 Figure 8. Paths of a bubble pair in the SIE regime (Ga = 90, Bo = 0.05). (a) Side view, with numbers indicating the time instant at the corresponding position; (b) bottom view; (c) snapshots of the three-dimensional streamlines past the two bubbles in the reference frame of the LB. The wake structure at t = 8 and 10 is also shown using the iso-surfaces ωy = ±3. In (a) -(b), the color and line codes are similar to those of figure 3.

Figure 9 .

 9 Figure 9. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (90, 0.05). (a) Rising speed (colored lines, left axis) and aspect ratio (black lines, right axis) of the LB (solid lines) and TB (dashed lines); (b) same for the horizontal velocity component; (c) vertical (red line, left axis) and transverse (blue line, right axis) separations; (d) inclination of the line of centres with respect to the vertical.

Figure 10 .

 10 Figure10. Three-dimensional streamlines past a bubble pair with (Ga, Bo) = (90, 0.5) during the early stage of their interaction (in the LB reference frame).

Figure 11 .

 11 Figure11. Evolution of the horizontal trace of the paths for a pair of bubbles with (Ga, Bo) = (90, 0.5); the solid and dash-dotted lines refer to the LB and TB, respectively. Successive periods of the path are coloured in red/blue/orange/green/purple/pink/cyan/black; numbers along the TB path indicate the corresponding time instant. The green trace in the bottom part is that of the corresponding isolated bubble, the initial (X, Z) position of which has been shifted for readability.

Figure 12 .

 12 Figure 12. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (90, 0.5). (a): rising speed (colored lines, left axis) and aspect ratio (black lines, right axis) of the LB (solid lines) and TB (dashed lines); (b) same for the horizontal velocity component; (c) vertical (red line, left axis) and transverse (blue line, right axis) separations. In (a) -(b), the color code is similar to that in figure 3.

Figure 13 .Figure 14 .

 1314 Figure 13. Bubble-wake interactions during three successive stages of the spiralling/zigzagging motion of a bubble pair with (Ga, Bo) = (90, 0.5). (a) -(c): bottom view of the path for 21.5 < t < 34.5, 34.5 < t < 46.7 and 72.6 < t < 85.6, respectively (the inset in (a) shows the two successive arrangements of the bubble pair in the horizontal direction). (d) -(f ) iso-surfaces ωy = +3 (red) and ωy = -3 (green) of the vertical vorticity at t = 32, 45 and 83, respectively (the top row displays two orthogonal side views, while the bottom row shows the bottom view).

Figure 15 .

 15 Figure 15. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (50, 0.5). (a): bottom view of the path of the LB (solid line) and TB (dashed line); (b): rising speed (colored lines, left axis) and aspect ratio (black lines, right axis); (c): Z-component of the horizontal velocity; (d): vertical (red, left axis) and horizontal (blue, right axis) separations; (e): inclination of the line of centres. The embedded snapshots in (b) show the iso-surfaces ωy = ±3 past the two bubbles during the time interval 23 t 26.

Figure 16 .

 16 Figure16. Bottom views of the paths for various Bo at Ga = 50 and Ga = 90. From left to right: Bo = 0.02, 0.05, 0.1, 0.2, 0.3, 0.5 and 1.0 (the latter with a yellow square corresponding to a head-on collision). Solid and dashed lines refer to the LB and TB, respectively. A shift has been applied to the initial X-position of each pair for readability.

Figure 18 .

 18 Figure 18. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (50, 0.3) undergoing initial angular deviations θ0 from 0 • to 2 • . (a): side view of the LB (solid line) and TB (dashed line) paths; (b): rising speed; (c): transverse velocity; (d): vertical (solid lines, left axis) and transverse (dashed lines, right axis) separations; (e): inclination of the line of centres.

Figure 19 .

 19 Figure 19. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (70, 0.3) undergoing initial angular deviations θ0 from 0 • to 2 • . For caption, see figure 18. The two snapshots straddling panels (b) and (d) show the iso-surfaces ωy = ±0.5 at two different instants of time for θ0 = 2 • .

Figure 20 .

 20 Figure 20. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (50, 0.3) for initial separations in the range 4 S0 8. (a) and (b): bottom and side views of the path of the LB (solid line) and TB (dashed line), respectively; (c) and (d): rising speed and transverse velocity of the two bubbles (same convention), respectively; (e) vertical (solid line, left axis) and transverse (dashed line, right axis) separations.

Figure 21 .

 21 Figure 21. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (50, 0.5) for initial separations in the range 8 S0 12. (a): bottom view of the paths; (b) -(c): vertical and transverse separations, respectively; (d) -(e): enlarged horizontal trace of the paths (left panel) and vertical vorticity distribution (identified with the iso-surfaces ωy = ±1) at different instants of time in various vertical or horizontal planes for S0 = 10 and S0 = 12, respectively.In (a) -(c) the red, blue and green lines correspond to S0 = 8, 10 and 12, respectively. In the insets of (d) -(e), the black (resp. purple) part of the paths corresponds to the rectilinear (resp. spiralling/zigzagging) stage that precedes (resp. succeeds) the realignment of the two bubbles at t ≈ 32.

Figure 22 .

 22 Figure 22. Evolution of several characteristics of a bubble pair with (Ga, Bo) = (50, 1.0), with the initial separation varied in the range 8 S0 28. (a): bottom view of the paths; (b) -(c): vertical and transverse velocities of the LB (solid line) and TB (dashed line), respectively; (d): vertical (solid line, left axis) and transverse (dashed line, right axis) separations. The squares in (b)-(d) indicate the point at which the two bubbles collide. In (a) a shift in the initial X-position has been applied to each pair for readability, while in (b) -(c) the results corresponding to S0 = 24 are not shown for the same reason.

Figure 23 .

 23 Figure 23. Influence of grid resolution on the evolution of a bubble pair with Ga = 90, Bo = 0.05. (a): Rising speed 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. Solid and dotted lines refer to simulations performed with ∆min = R/68 and ∆min = R/136, respectively.The inset in (a) provides a precise view of the residual differences between the rising speeds predicted on the two grids at t ≈ 20. The small high-frequency oscillations visible in this inset result from the dynamic changes of the grid near the bubble surface in this low-Bo high-Ga regime.

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

Acknowledgements.

The authors acknowledge the supports from the NSFC (Natural Science Foundation of China) under grants 11872296, 51636009 and 51806173, and from the Fundamental Research Funds for the Central Universities (xhj032021011-02). J.Z. also acknowledges the Young Talent Support Plan of Xi'an Jiaotong University.

Declaration of interests.

The authors report no conflict of interest.

Appendix A. Grid convergence in a high-Reynolds-number case

Several tests were reported in ZNM to assess the quality of the computational results. In particular, a grid convergence study with minimum grid spacings ∆ min = R/34, R/68 and R/136 showed that grid convergence is obtained with ∆ min = R/68 at least up to Ga = 30, which corresponds to Reynolds numbers of O(100). However in the present investigation, low-Bo bubble pairs with Ga = 90 reach Reynolds numbers of O(500). The corresponding boundary layer being typically twice as thin as in the previous case, an extra grid convergence test is in order. For this purpose, we select the parameters Ga = 90, Bo = 0.05, a configuration belonging to the SIE regime. In that case, the final Reynolds number of each bubble is about 470, so that with ∆ min = R/68 (resp.