1. Introduction
In the first part of this investigation (Zhang, Ni & Magnaudet Reference Zhang, Ni and Magnaudet2021), hereinafter referred to as Part 1), we analysed 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, millimetre-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 Part 1, 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 (Sangani & Didwania Reference Sangani and Didwania1993; Smereka Reference Smereka1993) assumed the bubbles kept 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 leads 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 (Esmaeeli & Tryggvason Reference Esmaeeli and Tryggvason1999; Bunner & Tryggvason Reference Bunner and Tryggvason2002; Esmaeeli & Tryggvason Reference Esmaeeli and Tryggvason2005; Yin & Koch Reference Yin and Koch2008). With significantly deformed bubbles, the dynamics of bubbly suspensions was observed to depend crucially on the Reynolds number. More specifically, the simulations of Esmaeeli & Tryggvason (Reference Esmaeeli and Tryggvason2005) 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 Bunner & Tryggvason (Reference Bunner and Tryggvason2003), 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 (Yin & Koch Reference Yin and Koch2008) bubbles were not allowed to deform. These restrictions were removed in the most recent simulations (Loisy, Naso & Spelt Reference Loisy, Naso and Spelt2017; Innocenti et al. Reference Innocenti, Jaccod, Popinet and Chibbaro2021). 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 (Innocenti et al. Reference Innocenti, Jaccod, Popinet and Chibbaro2021), surface tension was selected so as to keep the bubbles only mildly deformed, preventing the occurrence of path instability. Approximately 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 (Zenit, Koch & Sangani Reference Zenit, Koch and Sangani2001). In contrast, this ratio falls to values close to 
$2$ for larger bubbles exhibiting path instability (Risso & Ellingsen Reference Risso and Ellingsen2002; Riboux, Risso & Legendre Reference Riboux, Risso and Legendre2010). 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 (Alméras et al. Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015; Gvozdic et al. Reference Gvozdic, Alméras, Mathai, Zhu, van Gils, Verzicco, Huisman, Sun and Lohse2018). 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 analysed 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 Part 1 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 (Duineveld Reference Duineveld1998; Sanada et al. Reference Sanada, Sato, Shirota and Watanabe2009; Kong et al. Reference Kong, Mirsandi, Buist, Peters, Baltussen and Kuipers2019). 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 Chesters & Hofman (Reference Chesters and Hofman1982) 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 (Sanada et al. Reference Sanada, Sato, Shirota and Watanabe2009) 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 Kusuno, Yamamoto & Sanada (Reference Kusuno, Yamamoto and Sanada2019). However, all these studies considered Reynolds numbers of 
${O}(10\unicode{x2013}100)$ which belong to the moderately inertial regime examined numerically by Part 1 and, independently but with the same code, by Kusuno & Sanada (Reference Kusuno and Sanada2021). A noticeable exception is the work of Filella, Ern & Roig (Reference Filella, Ern and Roig2020) 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 Kusuno et al. (Reference Kusuno, Yamamoto and Sanada2019), Kusuno & Sanada (Reference Kusuno and Sanada2021) and Part 1 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 (Joseph et al. Reference Joseph, Fortes, Lundgren and Singh1986; Fortes, Joseph & Lundgren Reference Fortes, Joseph and Lundgren1987). Conversely, pairs of bubbles exhibiting a sufficient oblateness always collide and eventually coalesce, most of the time in the head-on 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\unicode{x2013}1000)$. In most cases, an isolated bubble would then either perform large-amplitude planar zigzags or follow a spiralling path. Similar to Part 1, we vary the control parameters and the initial conditions so as to cover a significant range of physical situations and analyse different evolution scenarios, with specific attention to 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 the Appendix). 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.
2. Problem statement and outline of the numerical approach
Since the problem was introduced in detail in Part 1, 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 1a). 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 
$\bar {S}(t)=S(t)/R$ and 
$\bar {S}_r(t) = S_r(t)/R$, respectively, and the angular deviation of the line of centres with respect to the vertical, 
$\theta (t)$ (figure 1b). In addition to the initial separation 
$\bar {S}_0=S_0/R$ (hereinafter set to 
$\bar {S}_0 = 8$ except specified otherwise), the flow and bubble dynamics are characterized by the Galilei and Bond numbers respectively defined as
\begin{equation} Ga = \rho_l g^{1/2}R^{3/2}/\mu_l, \quad Bo = \rho_l gR^2/\gamma, \end{equation}
where 
$\rho _l$ and 
$\mu _l$ are the density and viscosity of the carrying liquid, 
$\gamma$ 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 = \rho _lu_TR/\mu _l$ and 
$We = \rho _lu_T^2R/\gamma$, may be determined. Moreover, the Morton number 
$Mo=Bo^3/Ga^4= g\mu _l^4/\rho _l\gamma ^3$ is useful to identify the dynamics of bubble pairs in a given fluid, irrespective of the bubble size. In Part 1, 
$Ga$ and 
$Bo$ were varied in the range 
$10 \leq Ga \leq 30$ and 
$0.01 \leq Bo \leq 1.0$, respectively. The corresponding terminal Reynolds numbers were such that 
$10 \lesssim Re \lesssim 120$, so that an isolated bubble followed a rectilinear path. Here, we consider the parameter range 
$30 < Ga \leq 90$ and 
$0.02 \leq Bo \leq 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\approx 0.11$ and a Galilei number 
$Ga\approx 85$. Its aspect ratio 
$\chi$, defined as the ratio between the length of its major axis and its minor axis, is 
$\chi \approx 2.0$ and its terminal Reynolds number is 
$Re\approx 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 (Zenit & Magnaudet Reference Zenit and Magnaudet2008), corresponds to 
$Ga=90$ and 
$Bo\approx 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 (
$\chi \approx 3.1$). Under such high-
$Ga$ low-
$Bo$ conditions, an isolated bubble rises in a straight line only if its shape is close enough to a sphere. In contrast, bubbles such that 
$\chi \gtrsim 2.2$ follow a zigzagging or a (possibly flattened) spiralling path (Cano-Lozano et al. Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016). 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 the role of these trailing vortices in the dynamics of the bubble pair, especially in the possible direct interaction of the trailing bubble with the vortex pair emanating from the leading bubble. The above two examples indicate that the parameter range on which this study focuses corresponds to millimetre-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 centimetre-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 (Bhaga & Weber Reference Bhaga and Weber1981). This makes the interaction scenarios for such bubbles distinct from those considered here.
Sketch of the problem. (a) Initial configuration; (b) definition of the geometric parameters.
The results to be discussed below are obtained by solving the three-dimensional time-dependent Navier–Stokes equations valid throughout the flow (i.e. including capillary effects) with the open source flow solver Basilisk (http://basilisk.fr) described in Popinet (Reference Popinet2009, Reference Popinet2015). The numerical schemes employed in this solver are detailed in the above two references and are summarized in Part 1. 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 Part 1, 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 Part 1, 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 
$\varDelta _{min} = R/68$ close to the bubble interface and to 
$\varDelta = R/17$ in the wake. With this discretization, Part 1 established that the evolutions of the rise 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 Part 1. 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. With 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 the Appendix, 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. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016) (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 
$\sqrt {R/g}$.
3. Results in the reference case
3.1. Overview
Figure 2 summarizes the various interaction scenarios observed in the parameter range covered by the present investigation, 
$30< Ga\leq 90$, 
$Bo\leq 1.0$. The results from Part 1 at 
$Ga=30$ are also reported in figure 2(a) to serve as a reference. At this moderate 
$Ga$, the ASE scenario is observed up to 
$Bo\approx 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} \approx 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\unicode{x2013}0.4$ in figure 2a), a first oscillatory regime is observed, with the two bubbles zigzagging virtually in the same plane and independently from each other (left panel of figure 2c). 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}\approx 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 non-coplanar independent zigzagging (hereinafter abbreviated as NIZ). Clearly the CIZ–NIZ transition corresponds to a loss of the planar symmetry of the system. 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 2c). Last, the two bubbles collide beyond 
$Bo_{c} \approx 0.7$. For higher Galilei numbers (
$Ga\geq 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 of each other with only small-amplitude erratic horizontal excursions. This defines the small-amplitude independent erratic (SIE) regime.
Configuration map based on simulation results: (a) 
$(Ga, Bo)$-plane; (b) 
$(Re, \chi )$-plane. violet square: ASE regime; magenta bullet: small-amplitude independent erratic (SIE) regime; red triangle: collision regime; orange diamond, blue diamond, green diamond: coplanar independent zigzagging, non-coplanar independent zigzagging and interacting flattened spiralling 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 leading bubble and trailing bubble, 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-
$Mo$ lines corresponding to different liquids, with the GaInSn liquid metal (Galinstan) at the very bottom (
$Mo=1.37\times 10^{-13}$), then water (
$Mo=2.54\times 10^{-11}$ at 
$20\,^\circ {\rm C}$), then silicone oils T0 (
$Mo=1.8\times 10^{-10}$), T2 (
$Mo=1.6\times 10^{-8}$) and T5 (
$Mo=6.2\times 10^{-7}$) from bottom to top (see e.g. Zenit & Magnaudet Reference Zenit and Magnaudet2008 for the characteristics of these oils).
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\lesssim 0.05$ and 
$Ga\gtrsim 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, Magnaudet & Fabre Reference Tchoufag, Magnaudet and Fabre2014b). 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 (Tchoufag, Fabre & Magnaudet Reference Tchoufag, Fabre and Magnaudet2015).
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.
3.2. 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\approx 125$ and 
$\chi \approx 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\lesssim t\lesssim 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\lesssim 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 neighbourhood 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\approx 18\unicode{x2013}19$ (panel c). Remarkably, under the 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.
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,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.
The evolution of some characteristics of the motion is displayed in figure 4. The early evolution of the rise speed (panel a, left axis), is similar to that observed in the ASE scenario detailed in Part 1. That is, once the TB enters the wake of the LB, it is strongly accelerated by the corresponding ‘sheltering’ effect, while the rise 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 \lesssim t \lesssim 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 Part 1, the in-line arrangement becomes unstable and the TB starts to drift laterally, as the records of the transverse separation 
$\bar {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 rise speed, 
$V_{TB}$, 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\approx 15$.
Evolution of several characteristics of a bubble pair with 
$(Ga,Bo) = (50,0.3)$. (a) Rise speed (coloured 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 colour code is identical to that in figure 3.
At 
$t \approx 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 disturbance provided by the initial lateral drift. However, this influence quickly fades away and the lateral excursions reach a saturated periodic state at 
$t\approx 35$. Beyond this point, the TB kinematics are characterized by a rise speed 
$V_{TB}\approx 2.5$, a Strouhal number 
$St=f/V_{TB}\approx 0.029$ and a crest-to-crest amplitude 
$\bar {a}\approx 1.75$ (figure 3a), with maximum horizontal velocities (figure 4b) of approximately 
$0.15V_{TB}$. These characteristics superimpose onto those of the corresponding isolated bubble, confirming that the LB no longer influences the TB. The oscillating component of the LB motion follows a strikingly different evolution. Indeed, it develops much more gradually and figure 4(b) shows that the oscillations of its lateral velocity, 
$V_{r_{LB}}$, still have a slightly smaller amplitude than those of 
$V_{r_{TB}}$ at 
$t\approx 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_{TB}$, 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 figure 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 Hallez & Legendre (Reference Hallez and Legendre2011) for spherical bubbles, and the average inclination of the line of centres will be approximately 
$24^\circ$. 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 4b). However, there is a priori no reason for this phase shift to vanish eventually, given the quasi-independence of the two paths. Consequently, both 
$\bar {S}_r$ and 
$\theta$ will continue to oscillate slightly about their mean value.
Another perspective on 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, 
$\omega _y=\pm 0.5$, are plotted in figure 5. Two vortex threads in which 
$\omega _y$ periodically changes sign emerge in the TB wake right after its initial lateral drift. In contrast, no structure corresponding to the selected 
$|\omega _y|$-level is detected in the LB wake until 
$t\approx 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\approx 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.
Successive snapshots of the iso-surfaces 
$\omega _y = -0.5$ (green) and 
$\omega _y = 0.5$ (red) during the rise of a bubble pair with 
$(Ga,Bo) = (50,0.3).$
3.3. 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 rise speed and aspect ratios are 
$Re\approx 168$ and 
$\chi \approx 2.2$, respectively (CL16). The major two differences in the evolution of this bubble pair compared with the previous one may be appreciated from figure 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.
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 colour 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 6(c) displays the evolution of the bubbles shape and orientation, together with that of the streamline 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 (Ryskin & Leal Reference Ryskin and Leal1984; Zenit & Magnaudet Reference Zenit and Magnaudet2008). 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 bubble equator) is known to vary as 
$\chi ^{8/3}$ for large 
$\chi$ (Magnaudet & Mougin Reference Magnaudet and Mougin2007). 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 
$\chi _c\approx 2.2$ (Magnaudet & Mougin Reference Magnaudet and Mougin2007). 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\approx 5$ to 
$t\approx 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.
Evolution of several characteristics of a bubble pair with 
$(Ga,Bo) = (70,0.3)$. (a) Rise speed (coloured 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.
Panels (b-d) in figure 7 confirm that the dynamics of the tandem has reached a fully developed state soon after 
$t=30$. Indeed, the vertical separation stabilizes at 
$\bar {S}\approx 9.7$ and no longer varies. Owing to the phase shift between the two paths revealed by figure 6(a), the horizontal separation 
$\bar {S}_r$ experiences periodic large-amplitude oscillations with a reduced frequency 
$St\approx 0.031$. These oscillations make 
$\bar {S}_r$ vary from 
$2.25$ to 
$10$, which in turn results in an inclination of the tandem varying from 
$13^\circ$ to 
$45^\circ$. As figure 7(a) shows, the rise 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 rise speed reaching its maximum twice during a complete zigzag. The oscillations of the aspect ratio are enslaved to those of the rise 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 rise speed of the LB (similar to that of the isolated bubble), from 
$V_{LB}\approx 2.65$ for 
$t\lesssim 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 6b) 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 
$\chi _c$, non-axisymmetric disturbances around 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.
3.4. SIE regime
For 
$Ga\geq 70$ and 
$Bo\leq 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\lesssim 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 with 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.
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 
$\omega _y=\pm 3$. In (a,b), the colour and line codes are similar to those of figure 3.
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 2b and 9a), an oblateness for which the numerical results of Blanco & Magnaudet (Reference Blanco and Magnaudet1995) 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\approx 470$, see figure 2b), 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 
$\chi _c\approx 2.2$ (Magnaudet & Mougin Reference Magnaudet and Mougin2007). 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 
$\chi$) 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 (Auguste, Magnaudet & Fabre Reference Auguste, Magnaudet and Fabre2013) and linear stability analysis (Tchoufag, Fabre & Magnaudet Reference Tchoufag, Fabre and Magnaudet2014a) revealed 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.
Evolution of several characteristics of a bubble pair with 
$(Ga,Bo) = (90,0.05)$. (a) Rise speed (coloured 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 9 helps understand the reason and consequences of the unexpectedly long drift of the TB. In particular, panel (a) reveals that, while the rise speed of the LB has reached its final value at 
$t\approx 13$, its aspect ratio still increases significantly until 
$t\approx 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 
$\chi (t)$ at a given time is lower than what the quasi-static approximation of Moore (Reference Moore1965) based on the Weber number 
$We(t)$ built with the rise 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\lesssim 10$ according to figure 8c), 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 non-zero inclination of the bubble equator at 
$t=15$ in figure 8(c) confirm. This is why the TB continues to drift laterally until the uppermost position reached in the computation. A direct consequence of this sustained drift is that the long-term lateral separation reaches much larger values than in the previous regimes (
$\bar {S}_r\approx 15$ according to figure 9(c), to be compared with average values in the range 5–6 in the CIZ and NIZ regimes). This in turn translates into a significantly larger final inclination of the bubble pair, close to 
$55^\circ$ (figure 9d). In anticipation to the regime discussed below with the same 
$Ga$ but a ten times larger 
$Bo$ (see figure 12c), 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 
$\bar {S}$ below its initial value.
3.5. 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 \approx 175$ and 
$\chi \approx 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\approx 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 ASE-type stage.
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).
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\gtrsim 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^\circ$ 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^\circ$, 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^\circ$ at the end of the simulation.
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 position of which has been shifted for readability.
Figure 12(a) shows how the rise speed and aspect ratio of the two bubbles vary during their rise. Both quantities experience strong and rapid variations during the stage 
$15\lesssim t\lesssim 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\approx 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 rise 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.
Evolution of several characteristics of a bubble pair with 
$(Ga,Bo)=(90,0.5)$. (a) rise speed (coloured 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 colour code is similar to that in figure 3.
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\lesssim t\lesssim 33$, the vertical separation continuously decreases over time, finally stabilizing around 
$\bar {S}\approx 2.6$. As a result, beyond 
$t\approx 70$, the inclination of the line of centres (not shown) stabilizes around 
$\theta \approx 65^\circ$ with 
$\pm 5^\circ$ oscillations. This small final vertical separation is to be compared with the final values 
$\bar {S}\approx 10$ or even larger found in figures 4, 7 and 9. Most of the decrease of 
$\bar {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 for 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\approx 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 
$\bar {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 
$\bar {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 bubbles 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.
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 
$\omega _y = +3$ (red) and 
$\omega _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).
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 (respectively positive) in the upper (respectively 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.
Sketch illustrating the mechanism driving the horizontal precession of the TB path observed for 
$(Ga,Bo) =(90,0.5)$ (bottom view). (a) Deviation of a vortex pair in the presence of another vortex pair (the red and green dashed lines schematize the positive and negative vortex threads); (b–d) successive stages of the clockwise deviation of the TB as it approaches the LB and undergoes the influence of its wake. The blue arrows indicate the direction of the horizontal drift of the corresponding bubble/vortex pair.
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 
$\bar {S}_r \approx 0$ and 
$\theta \approx 0^{\circ}$ (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 (respectively increase) of the TB aspect ratio (respectively rise 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 non-zero 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$).
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) rise speed (coloured 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 
$\omega _y = \pm 3$ past the two bubbles during the time interval 
$23 \leq t \leq 26$.
This intense sequence is succeeded by a long transitional stage (
$30\lesssim t \lesssim 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 large-amplitude 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 rise 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\lesssim 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.
4. Influence of some parameters
4.1. Influence of the Bond number
The influence of the Bond number on the fate of the bubble pair was discussed in detail in Part 1 for 
$Ga\leq 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 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\leq 0.2$). Then, the dramatic change of the trace from 
$Bo=0.3$ to 
$0.5$ helps us 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 clearly visible for 
$Ga=90$. At this high 
$Ga$, bubble pairs with 
$Bo\lesssim 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.
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.
4.2. Influence of an initial angular deviation
In Part 1 it was shown that, for 
$Ga\leq 30$, a small non-zero initial inclination 
$\theta _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 
$\theta _0=0^\circ$ case. As a consequence, the DKT interaction which characterizes the dynamics for 
$Ga={O}(10)$ and 
$Bo\lesssim 0.2$ when 
$\theta _0=0^\circ$ no longer exists with an initial inclination as weak as 
$2^\circ$. The critical Bond number beyond which the two bubbles collide and eventually coalesce also increases significantly with 
$\theta _0$, owing to the longer time offered to the TB to leave the LB wake. Finally, when the two bubbles do not collide, their final horizontal separation is smaller than in the 
$\theta _0=0^\circ$ case, since the shear that drives the lateral drift of the TB is weaker, owing to the larger 
$\bar {S}$ at which this drift takes place. Here, we question the influence of 
$\theta _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^\circ$. 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 initially containing 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 (Mougin & Magnaudet Reference Mougin and Magnaudet2002; Tchoufag et al. Reference Tchoufag, Magnaudet and Fabre2014b). The orientation of this plane is arbitrary and is usually dictated by some initial disturbance. Here, due to the non-zero 
$\theta _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 
$\theta _0=0^\circ$, 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 
$\theta _0=2^\circ$, 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 
$\theta _0=0^\circ$, it rather stands in the CIZ regime (or the weakly three-dimensional NIZ regime) when 
$\theta _0=2^\circ$. 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 
$\theta _0=2^\circ$. 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.
Variations of the path geometry (bottom view) with the initial angular deviation 
$\theta _0$ for three typical bubble pairs which, in the reference case 
$\theta _0=0^\circ$, evolve in the CIZ, NIZ and IFS regimes (from top to bottom). Red, blue and green lines correspond to 
$\theta _0=0^{\circ }$, 
$1^{\circ }$ and 
$2^{\circ }$, 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 
$\bar {S}_r=\bar {S}_0\tan \theta _0$ apart, with the TB shifted in the 
$x>0$ direction.
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 non-zero 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 18(c). This earlier start has a direct consequence for the future rise of the LB: since 
$\bar {S}$ is larger than in the reference case during the TB drift (
$\bar {S}\approx 8$ instead of 
$\bar {S}\approx 5\unicode{x2013}6$ for 
$\theta _0=0^\circ$ according to figure 18d), 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 
$\bar {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 
$\bar {S}_r$, confirming the indication of figure 17. Indeed, the mean value of 
$\bar {S}_r$ which is close to 
$4.5$ when 
$\theta _0=0^\circ$ is reduced to 
$2.1$ for both 
$\theta _0=1^\circ$ and 
$\theta _0=2^\circ$ (figure 18d). The large-amplitude zigzags performed by the TB and their phase difference with the slowly growing oscillations of the LB make 
$\bar {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\lesssim \bar {S}_r\lesssim 2$ while the vertical separation is large (
$\bar {S}\gtrsim 9$) implies that the TB repeatedly comes back into the LB wake. Because of this, the dynamics of the TB remains 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^\circ$ for 
$t>50$ (figure 18e). This continuous interaction, combined with the large amplitude of the zigzags performed by the TB, makes 
$V_{TB}$ experience growing periodic variations with a relative magnitude of up to 
$5\,\%$ in the late stage; in contrast, the rise speed of the LB stays remarkably constant (and independent of 
$\theta _0$) beyond 
$t=20$ (figure 18b). In summary, the evolution of this bubble pair for 
$\theta _0\neq 0^\circ$ 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.
Evolution of several characteristics of a bubble pair with 
$(Ga,Bo) = (50,0.3)$ undergoing initial angular deviations 
$\theta _0$ from 
$0^{\circ }$ to 
$2^{\circ }$. (a) Side view of the LB (solid line) and TB (dashed line) paths; (b) rise speed; (c) transverse velocity; (d) vertical (solid lines, left axis) and transverse (dashed lines, right axis) separations; (e) inclination of the line of centres.
The influence of 
$\theta _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 
$\theta _0$ in the range 
$0^\circ$–
$2^\circ$. 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 
$\theta _0=1^\circ$ and 
$2^\circ$. 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\approx 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 (figure 19a,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 (
$\bar {S}\approx 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\approx 25$), so that the streamwise vortices soon develop downstream and reach the TB at 
$t\approx 31$ (second snapshot in panels b and d), although the two bubbles are widely separated at this stage (
$\bar {S}\approx 12,\bar {S}_r\approx 5.5$). Since the symmetry plane of the LB wake makes a small but non-zero angle with the 
$(x,y)$ plane, the momentum transferred to the TB during the collision with the streamwise vortices has a non-zero 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 
$\theta _0=0^\circ$.
Evolution of several characteristics of a bubble pair with 
$(Ga,Bo) = (70,0.3)$ undergoing initial angular deviations 
$\theta _0$ from 
$0^{\circ }$ to 
$2^{\circ }$. For caption, see figure 18. The two snapshots straddling panels (b) and (d) show the iso-surfaces 
$\omega _y = \pm 0.5$ at two different instants of time for 
$\theta _0=2^\circ$.
4.3. Influence of the initial separation in the zigzagging/spiralling regimes
So far, the initial separation between the two bubbles has been maintained at 
$\bar {S}_0 = 8$. We now examine how varying 
$\bar {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 
$\bar {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 
$\bar {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 
$\bar {S}_0=4$. The minimum reached by the vertical separation during the initial axisymmetric stage decreases significantly as 
$\bar {S}_0$ is reduced, from 
$\bar {S}\approx 5$ for 
$\bar {S}_0 = 8$ to 
$\bar {S}\approx 2.3$ for 
$\bar {S}_0 = 4$. Moreover, the smaller 
$\bar {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 Part 1, is that the smaller 
$\bar {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 
$\bar {S}$-minimum implies that the magnitude of the asymmetric disturbance induced by this drift in the vicinity of the LB increases as 
$\bar {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 
$\bar {S}_0\leq 5$ (
$t\approx 20$) than for larger initial separations (
$t \gtrsim 120$ for 
$\bar {S}_0 = 8$, see figure 4). The behaviour of the system for 
$\bar {S}_0>8$ may be extrapolated by combining the information provided by figure 20(e). First, since the 
$\bar {S}$-minimum reached during the initial stage increases with 
$\bar {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 
$\bar {S}_0$. Beyond this stage, the long-term value of 
$\bar {S}$ increases monotonically with 
$\bar {S}_0$ while the average lateral separation stays almost 
$\bar {S}_0$-independent. From this monotonic behaviour in the range 
$4\leq \bar {S}_0\leq 8$, it may reasonably be inferred that the system remains in the CIZ regime for 
$\bar {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\leq \bar {S}_0\leq 20$ (not shown): the system stays in the NIZ regime whatever 
$\bar {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.
Evolution of several characteristics of a bubble pair with 
$(Ga,Bo) = (50,0.3)$ for initial separations in the range 
$4\leq \bar {S}_0\leq 8$. (a,b) Bottom and side views of the path of the LB (solid line) and TB (dashed line), respectively; (c,d) rise speed and transverse velocity of the two bubbles (same convention), respectively; (e) vertical (solid line, left axis) and transverse (dashed line, right axis) separations.
We then examine the changes in the dynamics of the pair with 
$(Ga, Bo) = (50, 0.5)$ when the initial separation is varied from 
$\bar {S}_0 = 8$ to 
$\bar {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 
$\bar {S}_0 = 8$ to a clear IFS evolution for 
$\bar {S}_0 = 10$. Then, for 
$\bar {S}_0 = 12$, the two paths are almost coplanar, so that the observed dynamics looks 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 
$\bar {S}_0 = 10$ and 
$\bar {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\approx 31\unicode{x2013}32$, instead of 
$t\approx 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 
$\bar {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 
$\bar {S}_0$ is, the longer it takes for the trailing vortices in the LB wake to grow.
Evolution of several characteristics of a bubble pair with 
$(Ga,Bo) = (50,0.5)$ for initial separations in the range 
$8\leq \bar {S}_0 \leq 12$. (a) Bottom view of the paths; (b,c) vertical and transverse separations, respectively; (d,e) enlarged horizontal trace of the paths (left image) and vertical vorticity distribution (identified with the iso-surfaces 
$\omega _y=\pm 1$) at different instants of time in various vertical or horizontal planes for 
$\bar {S}_0 = 10$ and 
$\bar {S}_0 = 12$, respectively. In (a–c) the red, blue and green lines correspond to 
$\bar {S}_0 =8$, 
$10$ and 
$12$, respectively. In the insets of (d,e), the black (respectively purple) part of the paths corresponds to the rectilinear (respectively spiralling/zigzagging) stage that precedes (respectively succeeds) the realignment of the two bubbles at 
$t\approx 32$.
For 
$\bar {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\approx 31$. However, here, the two bubbles start to spiral right away, while for 
$\bar {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 (
$\bar {S}\approx 4.5$ instead of 
$\bar {S}\approx 2$, according to figure 21b). As a result, the horizontal separation between the two bubbles stays much smaller than in the reference case, even at large time (see figure 21c, where the 
$\bar {S}_r$-minima are seen to be now approximately 
$\bar {S}_r\approx 0.7$, instead of 
$\bar {S}_r\approx 3$ for 
$\bar {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\approx 32$, the evolution observed when 
$\bar {S}_0 = 12$ differs dramatically from that above. Indeed, as the orientation of the black and purple segments in the left panel of figure 21(e) indicates, 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 21e). Moreover, the vertical separation is significantly larger than with 
$\bar {S}_0 = 10$ (
$\bar {S}\approx 7$ instead of 
$\bar {S}\approx 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$).
4.4. 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 analyse its influence in the collision regime identified for 
$\bar {S}_0=8$ in figure 2. Some characteristics of the dynamics observed by increasing the initial separation up to 
$\bar {S}_0=28$ are displayed in figure 22. The two bubbles are found to collide for 
$\bar {S}_0 \leq 20$, the collision time increasing significantly with 
$\bar {S}_0$ (figure 22b–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 
$\bar {S}_0$, suggesting that the two bubbles would also eventually collide later for these large initial separations. Present observations are reminiscent of those of Stewart (Reference Stewart1995) and Brücker (Reference Brücker1999) 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.
Evolution of several characteristics of a bubble pair with 
$(Ga,Bo) = (50,1.0)$, with the initial separation varied in the range 
$8\leq \bar {S}_0\leq 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 
$\bar {S}_0=24$ are not shown for the same reason.
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 
$\bar {S}_0 = 8$. Path instability of the isolated bubble sets in at 
$t\approx 25$ and saturates at 
$t\approx 70$ according to figure 22(c). Consistently, all pairs with 
$\bar {S}_0 > 8$ perform oscillatory paths. For 
$\bar {S}_0 > 16$, the saturated path of the LB takes the form of a planar zigzag or (for 
$\bar {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 
$\bar {S}_0$, approximately from 
$2.5$ for 
$\bar {S}_0 =20$ to 
$1.5$ for 
$\bar {S}_0 =28$, and reaches minima close to zero at regular time intervals (figure 22d). These characteristics, together with the aforementioned decrease of 
$\bar {S}$ over time, reinforce the view that a collision happens in all cases. In figure 22(b), the records of 
$V_\mathrm {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 (
$\bar {S}_0=8$) or oblique (
$\bar {S}_0=16$) configuration, in which case 
$V_\mathrm {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 figure 22(b,c) for 
$\bar {S}_0 =28$ is of interest. Beyond 
$t\approx 70$, the vertical and transverse velocities of both bubbles appear to have reached a saturated state. However, the rise 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 
$\chi \approx 2.55$ and 
$Re\approx 77$, respectively. Adoua, Legendre & Magnaudet (Reference Adoua, Legendre and Magnaudet2009) 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 
$Re$-dependent 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 Part 1. 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 
$\boldsymbol {F}=\mathcal {C}_L\mathcal {M}\boldsymbol {U}\times \boldsymbol {\omega }_\infty$, with 
$\mathcal {M}$ the mass of fluid that can be contained within the bubble volume, 
$\boldsymbol {\omega }_\infty$ the upstream vorticity and 
$\boldsymbol {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 shear-free condition), the lift coefficient on a spherical bubble is 
$\mathcal {C}_L=0.5$ (Auton Reference Auton1987). Still in this limit, 
$\mathcal {C}_L$ is an increasing function of the aspect ratio, with 
$\mathcal {C}_L\approx 1.45$ for 
$\chi =2.55$. These predictions are only marginally altered by finite-
$Re$ effects as far the relative shear rate 
$Sr=R\|\boldsymbol {\omega }_\infty \|/\|\boldsymbol {U}\|$ is such that 
$Sr\gtrsim 0.1$. In contrast, for smaller shear rates, the lift force changes sign beyond 
$Re\approx 37$ for this specific oblateness, and an interpolation of the results of Adoua et al. (Reference Adoua, Legendre and Magnaudet2009) predicts 
$\mathcal {C}_L\approx -0.6$ for 
$Re=77$ and 
$Sr=0.02$.
The manner in which 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 (Brücker Reference Brücker1999; Mougin & Magnaudet Reference Mougin and Magnaudet2006). The upstream vorticity 
$\boldsymbol {\omega }_\infty$ ‘felt’ by the TB results from the axisymmetric velocity defect in the LB wake, and this upstream (
$\approx$ horizontal) vorticity is responsible for an additional amount of streamwise (
$\approx$ vertical) vorticity in the TB wake, compared with the isolated configuration. When the vertical separation between the two bubbles is moderate, say typically 
$\bar {S}\lesssim 10$, the TB stands in the near wake of the LB. There, the relative shear rate is large enough (
$Sr\gtrsim 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 
$\bar {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 
$\boldsymbol {\omega }_\infty$ vanishes) than when it is close to its extremities. This may be confirmed from the curves corresponding to 
$\bar {S}=20$ and 
$\bar {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 
$\bar {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 
$\bar {S}$ (as seen in figure 22(d)), making a collision unavoidable at some point.
5. 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 \leq Ga \leq 90$, 
$0.02 \leq Bo \leq 1.0$. In that range, provided the Bond number exceeds a 
$Ga$-dependent threshold 
$Bo_o(Ga)$ decreasing from 
$Bo_o\approx 0.3$ for 
$Ga=40$ to 
$Bo_o\approx 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\lesssim 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 Part 1. Beyond a second threshold ranging from 
$Bo_c\approx 0.5$ for 
$Ga\leq 40$ to 
$Bo_c\approx 0.7$ for 
$Ga\geq 50$, the two bubbles always collide. We did not attempt to compute the collision process in detail, but analysed 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\gtrsim Bo_o$ and 
$Ga\lesssim 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 CIZ regime. Increasing the Bond number beyond the upper limit of this regime, one first encounters a 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 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 re-enter 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\gtrsim 70$), the CIZ regimes was not observed. Instead, for low Bond numbers (
$Bo\lesssim 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 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 if slight, a non-zero 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\approx 0.5$), i.e. the system stands in the IFS regime for 
$\theta _0=0^\circ$. 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. blackFor this reason, the IFS regime no longer subsists for 
$\theta _0=2^\circ$ within the 
$(Ga,Bo)$ range we explored. The reverse happens when the bubble oblateness is moderate (
$Bo\approx 0.3$), i.e. the system stands in the CIZ regime for 
$\theta _0=0^\circ$. 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 
$\bar {S}_0$ only change the geometry of the final arrangement in the CIZ and NIZ regimes, i.e. for 
$Bo\lesssim 0.4$, reducing (respectively increasing) the inclination of the bubble pair if the initial separation is increased (respectively reduced). Things are more complex in the IFS regime, i.e. for 
$Bo\approx 0.5$ and 
$Ga\gtrsim 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 further direct interaction happens. Increasing 
$\bar {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 
$\bar {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 
$\bar {S}_0$. With 
$\bar {S}_0\lesssim 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 with the wake-induced suction to produce a sufficiently large deviation of the TB avoiding collision. For larger 
$\bar {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),\chi \gtrsim 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 Adoua et al. (Reference Adoua, Legendre and Magnaudet2009). As a consequence, the TB stays within the LB wake even for large 
$\bar {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 non-intuitive phenomenology. Beyond the role of the control parameters 
$Ga,\ Bo,\ \theta _0$ and 
$\bar {S}_0$ identified in Part 1 for 
$Ga\leq 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.
The 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 the 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 
$\theta _0=0^\circ$ is limited by the experimental uncertainly on 
$\theta _0$, which can hardly be reduced to less than 
$1^\circ$. 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 a 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 per cent. 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 hundred 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 time 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 millimetre-size bubbles.
Acknowledgements
The authors acknowledge support 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. Grid convergence in a high Reynolds number case
Several tests were reported in Part 1 to assess the quality of the computational results. In particular, a grid convergence study with minimum grid spacings 
$\varDelta _{min} = R/34,\ R/68$ and 
$R/136$ showed that grid convergence is obtained with 
$\varDelta _{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 approximately 
$470$, so that with 
$\varDelta _{min} = R/68$ (respectively 
$R/136$), approximately three (respectively six) grid cells stand within the boundary layer. If grid convergence is reached, the rise 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 Part 1 for a detailed discussion). The test was carried out on the reference grid with 
$\varDelta _{min} = R/68$ and on a twice as fine grid with 
$\varDelta _{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 \approx 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\times 10^{-4}$.
The results of this test are summarized in figure 23. The inset in panel (a) reveals that the rise 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 rise speed, is approximately 
$3.5\,\%$ at this Reynolds number (Moore Reference Moore1965). Therefore, it can be concluded that this contribution, while small, is properly captured with 
$\varDelta _{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 rise 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\lesssim t\lesssim 10$). However, this shorter 
$\bar {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 rise 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\approx 11$. The difference stops increasing at the end of the TB drift but persists in the long term, making the final 
$\bar {S}$ slightly larger on the finer grid.
Influence of grid resolution on the evolution of a bubble pair with 
$Ga = 90,\ Bo = 0.05$. (a) Rise 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 
$\varDelta _{min} = R/68$ and 
$\varDelta _{min} = R/136$, respectively. The inset in (a) provides a precise view of the residual differences between the rise speeds predicted on the two grids at 
$t\approx 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.
