The interaction of multiple bubbles in a Hele-Shaw channel

We study the dynamics of two air bubbles driven by the motion of a suspending viscous fluid in a Hele-Shaw channel with a small elevation along its centreline via physical experiment and numerical simulation of a depth-averaged model. For a single-bubble system we establish that, in general, bubble propagation speed monotonically increases with bubble volume so that two bubbles of different sizes, in the absence of any hydrodynamic interactions will either coalesce or separate in a finite time. However our experiments indicate that the bubbles interact and that an unstable two-bubble state is responsible for the eventual dynamical outcome: coalescence or separation. These results motivate us to develop an edge-tracking routine and calculate these weakly unstable two-bubble steady states from the governing equations. The steady states consist of pairs of `aligned' bubbles that appear on the same side of the centreline with the larger bubble leading. We also discover, through time-simulations and physical experiment, another class of two-bubble states which, surprisingly, consist of stable two-bubble steady states. In contrast to the `aligned' steady states, these bubbles appear either side of the centreline and are `offset' from each other. We calculate the bifurcation structures of both classes of steady states as the flow-rate and bubble volume ratio is varied. We find that they exhibit intriguing similarities to the single-bubble bifurcation structure, which has implications for the existence of $n$-bubble steady states.


I. INTRODUCTION
Identifying invariant objects (steady states, periodic orbits, invariant tori etc.) of highdimensional, nonlinear systems and how they influence the transient dynamics is crucial in understanding how a system evolves towards an eventual dynamical outcome. One approach to identify these objects is to perform a number of initial-value problems (IVP), either experimentally or theoretically, and observe how the system behaves. The inherent disadvantage of this approach is that the outcome is binary; either the system settles to a stable invariant object, or long-term transient behaviour emerges. To capture unstable invariant objects, bespoke techniques are required, for example edge-tracking [16] or parameter continuation [19,25]. Although these invariant objects may be unstable, they still influence the dynamics of the system in a crucial way that would remain hidden in an IVP. In highly complex systems, such as the transition to turbulence in pipe flow [13,15,29,30], these invariant solutions are often of high dimensionality and difficult to compute. We surmise, however, that these ideas are applicable to a large range of nonlinear systems and can be applied to systems which, although nonlinear and high-dimensional, are more amenable to theoretical and experimental analysis.
As a model 'playground' to test these ideas we consider the steady state structure and transient dynamics of two finite air bubbles propagating in a Hele-Shaw channel with a prescribed depth perturbation when the surrounding fluid is extracted at a constant flow rate (see figure 1). In a previous work [9], we showed that a single bubble may break up into two (or more) bubbles depending on its initial spatial configuration and on the flow rate and that, post breakup, the bubbles may either merge back into a single or compound bubble or separate indefinitely (see figure 2). A key result of this study was that the post breakup dynamics were strongly influenced by the existence of weakly unstable steady states that are specific to the two-bubble system. It was hence hypothesised that the complexity of the dynamics may increase with the number of bubbles, owing to the increase in the number of underlying (stable or unstable) steady states of the system.
A feature of this system is that the topology of the system changes when a bubble breaks up or two (or more) bubbles coalesce. Following such topological events, a different family of invariant solutions influence the transient dynamics. For a given system of, say, n-bubbles we might expect the steady states of the system to be related to the steady states of the lower-order 1, 2, · · · n − 1-bubble systems in such a way that a hierarchy of 1, 2, · · · n-bubble states can be constructed from smaller bubble systems. The broad phenomenon of 'lowerorder' states interacting to form new coherent structures has been seen in other physical systems. For example, the interaction of solitons in water waves [see, for example 5] and nonlinear optics [see, for example 1], spatially localised states in convection systems [see, for example 24] and oscillons in granular particulate flow [see, for example 36]. A particular anomaly of our system is that we cannot smoothly move from a n-bubble state to a m-bubble state by continuation or branch-switching methods because the topologies of the systems are different. How the steady states of n and m-bubbles relate to each other is therefore nontrivial and this system represents a rather different example of interacting localised states, from the previously highlighted.
The propagation of finite-air bubbles in a Hele-Shaw channel of uniform depth is a classical problem in fluid dynamics with a long and rich history. Transient behaviour and steady propagation modes have been investigated extensively in the case of a single-bubble using a mixture of analytical and numerical techniques [see, for example 10,17,20,[32][33][34] and experiments [see, for example 18,21,22,31,38,39]. If a depth-perturbation is added to the bottom of the channel as shown in figure 1, the range of existence and stability of steady propagation modes changes dramatically, as mapped out by Franco-Gómez et al. [7,8], Gaillard et al. [9], Keeler et al. [14]. The solution branches interact in a highly nontrivial manner, resulting in a number of different bifurcations and regions of bi-stability in the system; features absent when there is no geometric perturbation in the channel. Recently it has been shown that the transient behaviour of a single-bubble in a perturbed Hele-Shaw channel is heavily influenced by so-called 'edge-states' of the system whose stable and unstable manifolds separate different dynamical outcomes [9,14].
Although multi-bubble steady propagation modes in a Hele-Shaw channel have been studied in unperturbed channels, see, for example [4,10,20,37], these works have focused on steady state solution construction at zero surface tension and their significance to the underlying dynamics, including stability results, have not been investigated. Dipole models have been proposed to understand the dynamics of multiple bubbles in an infinite, unbounded Hele-Shaw cell, [11,26,28], by treating the bubbles as small and circular, and forming a system of ordinary differential equations describing the position of the individual bubbles based on interactions between each of them. The dipole model in an infinite, unbounded Hele-Shaw cell of uniform depth predicts that a single row of identical bubbles is neutrally stable but is prone to instability if 'nudged' out of line. Also, relevant to this study, two rows of identical bubbles, located symmetrically about the horizontal centreline, is also neutrally stable, whilst two rows of bubbles which are located asymmetrically about the horizontal centreline is unstable, see [26]. We remark that no stable multiple-bubble states have been observed in other confined systems and that in general the bubbles will always either separate or coalesce [21,22,27].
In this paper we concentrate on a two-bubble system in a depth-perturbed Hele-Shaw channel and investigate the existence of steady states and their dependence on the flow-rate and bubble volume. We calculate the two-bubble solution structure and find a number of two-bubble steady states, each playing a unique role in the underlying transient dynamics as the system parameters are varied. Surprisingly, we find that a stable steady state exists with a bubble on either side of the centreline and the smaller bubble leading. Furthermore, by comparing the two-bubble and single-bubble bifurcation diagrams, we uncover an underlying solution structure that may have implications for the existence of n-bubble steady states in general. We also make the observation that the dynamics of the two-bubble system are not necessarily dominated by the larger bubble, but rather it is the leading bubble that has the largest influence on the system. The paper is organised as follows. In § II and § III we present the experimental and numerical methods used to investigate the dynamics of the system. In § IV A we summarise the known results of the single-bubble system and extend these to explore the relationship between bubble speed and volume, which is fundamental to understanding the theoretical construction of two-bubble states. We then describe two classes of two-bubble states; aligned states where the bubbles have a similar vertical offset ( § IV B 1), and offset states where the bubbles are staggered on either side of the rail ( § IV B 2). Next, in § IV C, we compare the solution structures of the two-bubble system to the single-bubble system before we conclude with a discussion of the implications of our results for n-bubble systems ( § V).
Fluid is extracted a constant flux, Q * , at one end, so that the bubbles propagate down the channel. The perturbation takes the form of a constant height and width rectangular rail at the bottom of the channel.

II. EXPERIMENTAL METHODS
We performed experiments in which two bubbles propagated through the channel from prescribed initial configurations imposed prior to flow initiation. The experimental Hele-Shaw channel presented in figure 3 has been comprehensively described by [9]. Thus, we only recall the salient details here. The channel consisted of two float glass plates separated by walls (strips of stainless steel shim), which were accurately positioned to make a channel of length L * = 170 cm, width W * = 40 ± 0.1 mm and height H * = 1.00 ± 0.01 mm, with an aspect ratio α = W * /H * = 40. The channel was sealed with clamps and levelled horizontally to within 0.03 • . A centred rail of width w * = 10.0 ± 0.1 mm and thickness h * = 24 ± 1 µm consisted of a translucent adhesive tape strip bonded to the bottom glass plate, see figure 3(b).
The channel was filled with silicone oil (Basildon Chemicals Ltd) of dynamic viscosity µ = 0.019 Pa.s, density ρ = 951 kg/m 3 and surface tension σ = 21 mN/m at the laboratory temperature of 21 ± 1 • C, and connected to oil reservoirs through inlet and outlet ports located at the upstream and downstream ends of the channel, respectively (see figure 3(a)). Flow in the channel was imposed by injecting oil through the inlet port with constant volume flux Q * using a bank of three syringe pumps, and letting oil escape through the outlet port. Air bubbles were generated by injecting prescribed volumes of air in the channel through an air port positioned slightly downstream of the inlet port; see Gaillard et al. [9] for details on the bubble generation protocol. Once formed, the bubbles were propagated through a centring device consisting of a section of channel of reduced width followed by a region of linear expansion, as shown schematically in figure 3(a). Experiments were performed with pairs of bubbles, each of prescribed area as measured from above, which were arranged in reproducible initial configurations in terms of their shapes and relative positions. We distinguish 'aligned' initial bubble configurations from 'offset' configurations. The former correspond to axially aligned bubbles with both bubbles either positioned symmetrically about the channel centreline ('on-rail') or asymmetrically ('off-rail') but on the same side of the rail ( figure 3(a)). In the 'offset' configuration, two off-rail bubbles are positioned on opposite sides of the rail as shown schematically in figure 4. These initial bubble configurations were generated using two different experimental protocols described in appendix A.
Bubbles were propagated from their initial configuration at a constant dimensionless flow rate Q = µU * 0 /σ where U * 0 = Q * /(W * H * ) is the average oil velocity in an equivalent channel without the rail. The bubbles were filmed in top-view using a CMOS camera mounted on a motorised stage, which translated at a constant velocity value chosen to ensure that the bubbles remained within the field of view of the camera for the duration of the experiment. We refer to each initial bubble with a numerical index in order of decreasing size, so i = 1 corresponds to the largest bubble. The projected area A * i (i = 1, 2) and centroid position of each bubble were measured from the bubble contour detected using an edge detection algorithm. The distance between the two bubbles is quantified by D = 2D * /W * where D * is the dimensional distance between the centroids of the two bubbles. Unless otherwise specified, the combined bubble size is A total = A 1 +A 2 = 0.54 2 π which is the size of the single bubbles used in Gaillard et al. [9], where A i = A * i /(W * /2) 2 (i = 1, 2) is the non-dimensional area of each bubble. We investigated initial configurations with either the larger or smaller bubble in the lead position and report results in terms of the initial fractional size of the lead bubble given by the ratio A r = A lead /A total , where A lead may equal A 1 or A 2 depending on the initial order of the bubbles. For convenience, we sometimes also refer to the ratio of the bubble areas in the form A 1 : A 2 . For example, bubbles of volume ratio 2:1 have A r = 2/3 if the larger bubble is (initially) leading and A r = 1/3 otherwise.
In the experiment, we only measure the projected area directly during bubble propagation. For a fixed volume of injected air, the projected area of the bubble can decrease sharply when flow is initiated because of air compression which increases with flow rate as the associated pressure head increases. However, the bubble retains an approximately constant projected area during each experiment, with a small increase of less than 7% at the highest flow rates. Conversely, the presence of lubricating oil films separating the bubble from the top and bottom plates, whose thickness increases with increasing flow rate, tends to increase the projected area of the bubble. These effects are discussed in detail in Gaillard et al. [9] and a suitable calibration of the injected volume of air was performed to obtain propagating bubbles with the required prescribed areas A i (i = 1, 2).

III. MATHEMATICAL MODEL AND NUMERICAL METHODS
The depth-averaged model for the propagation of multiple bubbles in our Hele-Shaw channel has been previously described and we only summarise its key features below. Our approach extends that of [23] to account for a non-uniform channel height and has been used extensively in studies of the propagation of a semi-infinite air finger [6,35], single closed air bubbles [7,8,14] and most recently single and multiple air bubbles [9]. We use the model to compute steady states of the system, calculate their linear stability and perform numerical time-simulations.
We work in a frame moving with the centroid position of the entire collection of bubbles and non-dimensionalise the physical system shown in figure 3 using W * /2 and H * as characteristic length scales in the (x * , y * ) plane and z * direction, respectively, and U * 0 = Q * /(W * H * ) as the velocity scale. The resulting nondimensional computation domain is shown in figure 4.
The two-dimensional depth averaged lubrication model reduces to an equation for the pressure in the fluid domain where the mobility b(y) represents the variable depth of the channel, modelled as a smoothed tanh profile where h = h * /H * and w = w * /W * are the non-dimensional height and width of the rail respectively; and s sets the 'sharpness' of the sides of the rail, as shown in figure 4. We use h = 0.024 and w = 0.25 consistent with experiments and we choose s = 40 [35]. We impose no-penetration conditions on the upper and lower walls, which yield p y = 0 on y = ±1. The pressure is fixed to zero at the inflow, and a non-zero constant at the outflow to ensure the dimensionless volume flux is consistent with the inflow dimensionless volume flux. Equations are solved in the reference frame moving at velocity U and we assume that the air bubble fills the height of the channel so that the kinematic boundary conditions on the contour of each bubble denoted by R i (where i = 1, 2, . . . indicates the i th bubble in decreasing size order) is given by where n i is the unit normal vector directed away from the i th bubble and U = (U b , 0) is the velocity of the centre of mass of the system along x. The centre of mass speed, U b , is an unknown in the problem which is obtained by requiring that the x coordinate of the centre of mass of the system remains at zero. The dynamic boundary condition on each bubble is where κ denotes the curvature of the bubble in the (x, y) plane and the effects of the variable depth on the transverse curvature are accounted for by the 1/b(y) term. The pressure p i in each bubble is not known a priori and is determined by ensuring that the dimensionless bubble volume V i remains constant, where the volume V i is defined via where Γ i is the interior of bubble i and ∂Γ i its bounding curve. The total bubble volume V total = V 1 + V 2 is set to 0.54 2 π unless specified otherwise, consistent with the experiments of [9]. In the model, the dimensionless volume and area of a given bubble are almost identical because we assume that the bubble fills the entire channel height, which differs from 1 by at most 2.5 %. In the experiment, we measure the size of the experimental bubbles by their projected areas A i (i = 1, 2) because the air volume required to yield a bubble of fixed projected area varies with flow rate as discussed in §II. We solve the system of equations, (1)-(5) on the domain shown in figure 4 to determine p, R i , p i and U b and we use the flow rate Q and bubble volumes V i as control parameters. The spatial discretisation of the equations is obtained by using a finite-element method, see appendix §B 1. When performing time-simulations, we use a procedure detailed in appendix §B 2 to account for the topology changes that may occur, such as bubble breakup and coalescence events. Stable and unstable steady solutions of the governing equations are calculated using Newton's method. Convergence of this method requires a good initial guess for the bubble configuration. For stable steady states, an initial guess can be obtained by performing a time-simulation from an initial condition where the system converges towards the stable state. For unstable steady states however, finding a good initial guess requires bespoke methods for each individual state, see appendix B 3. Once a stable or unstable steady state has been identified for a given set of control parameters, we use continuation methods to map the solution space as the control parameters are varied.

A. Two-bubble metrics
We characterise the two-bubble system by the coordinates (x i , y i ) of the centroids of each bubble given by where A i is defined in equation (5). From these, we compute the distance D between the centroids of each bubble, given by as well as the y-coordinate of the centre of mass of the system defined by which we refer to as the offset of the two-bubble system. In time-simulations, these quantities will evolve as functions of time, x i (t), y i (t), D(t) and Y (t), and will be used to characterise the state of the system. In numerical time-simulations, the initial shape of each bubble was chosen to be an ellipse with contour coordinates In all the numerical time-simulations presented in this paper, the volume ratio was V 1 : V 2 = 2 : 1 and we chose initially slender bubbles with d 1 = 0.3 and d 2 = 0.2 so that the bubbles did not break up before they interacted. The values of i = V i /πd i were set to ensure the prescribed volumes.

A. Single-Bubble Systems
Before we discuss two-bubble systems, we present an overview of the steady propagation of single bubbles and examine the influence of bubble volume over our range of interest. Figure  5 presents theoretical, panel (a), and experimental results, panels (b,c), for the dimensionless speeds, U b , of individual bubbles as functions of dimensionless flow rate Q. The different colours correspond to different bubble volumes and we find that the structure of the steadily propagating solutions in the theoretical model is independent of bubble volume within the range investigated. We label the different solution branches as in [14] and [9]. For our region of interest, there are three distinct solutions: a stable asymmetric bubble, AS1, that exists for all flow-rates; an unstable symmetric double-tipped bubble, S1, that exists for small flow-rates; and an alternative symmetric bubble, S2/S3, that exists for larger flow-rates and can be stable or unstable. Inset snapshots in figure 5(a), show the shapes of the bubble at flow rates indicated by dots on the solution branches. For an intermediate range of flow rates, the system is bistable with both symmetric and asymmetric, stable propagation modes available.  The experimental data presented in (b) is reproduced from [9] with permission. The inset panels in (a) correspond to bubble profiles specified by solid circular markers on the branches. The bubbles have volumes V = π0.54 2 , π0.46 2 , π0.35 2 , indicated by the different colours of the branches. In the theoretical results, the flow rate Q s indicates where a 'switchover' occurs in the relative speeds of larger and smaller bubbles. There is no evidence for the existence of Q s in the experiments. The hollow circular markers in (b) correspond to the symmetric S3 state and the solid circular markers in (b) and (c) correspond to the asymmetric AS1 state.
The experimental data shown in figure 5(b) correspond to the symmetric and asymmetric states in the bistable regime and we find the following general trends in both experiments and the theoretical model: (i) the dimensionless speed U b = U * b /U * 0 , of the bubble relative to the average speed of the surrounding fluid, typically increases with flow rate (note that for high flow rates the relative bubble speed saturates for the symmetric state and decreases slightly for the asymmetric state); (ii) for a fixed flow rate, the dimensionless speed increases with bubble volume, for the range of volumes investigated here. At low flow rates in the theoretical model, however, the opposite trend is observed.
The result that larger bubbles travel faster at a fixed flow rate follows from the theoretical analysis of [34], under the assumption of fixed bubble width. The same theory predicts that for a fixed volume, wider bubbles travel faster. The explanation for these results is that either increasing bubble width for fixed volume, or increasing volume for fixed width leads to increased viscous dissipation. Increased viscous dissipation is balanced by an increase in the work done by fluid pressure on the bubble, which results in a higher local fluid pressure gradient. The increased pressure gradient leads to a higher local fluid velocity around the bubble, which leads to faster bubble speeds via the kinematic condition, equation (3). Related results have also been found for buoyant rise of bubbles in Hele-Shaw cells [22] in which the bubble speed also increases with bubble width. In that case, however, the increase in viscous dissipation is balanced by an increase in the work done by the buoyancy force, which itself increases with bubble volume. Similar arguments can explain why in the bistable regime the asymmetric bubble travels faster than the symmetric one, for a fixed flow rate. The symmetric bubble spans the rail which means that it displaces a smaller area of fluid within each cross-section as it propagates, leading to lower dissipation, lower local pressure gradient and hence a lower bubble speed.
At lower flow rates (Q < 0.02), there is a qualitative disagreement between the model and the experimental results. In the experiments, the relative propagation speed of single bubbles is an increasing function of the bubble size at all values of the flow rate investigated, see figure 5(c). In the model, for both the asymmetric and symmetric solutions there is a critical flow rate below which the relative bubble speed decreases as the bubble volume increases. The value of Q where this trend 'switches' over is denoted Q s . This discrepancy has been previously noted by [22] who states "It is also clear also that the theory overestimates the bubble velocities for the smaller widths, having the wrong behaviour as D (the diameter) → 0" [22, pp.108].
We accept, therefore, that for small bubbles and low-rates the model does not reflect the experiments and we confine the majority of our analysis to Q > Q s where the model and experiment predict the same speed-volume relationship.

B. Two-Bubble Systems
From the results for single bubbles, §IV A, we know that two bubbles of different volumes will always travel at different speeds when the surrounding fluid moves at a fixed flow rate Q. Hence, in the absence of any hydrodynamic interactions and irrespective of the initial separation, two different bubbles will either coalescence in finite time or separate indefinitely depending on whether the slower bubble is initially ahead or behind the faster, respectively. In this section, we demonstrate that bubbles interact hydrodynamically, which results in the existence of stable and unstable two-bubble steady modes of propagation. As far as we are aware, such modes have not been observed for bubble of different sizes in related confined systems. In those cases, the bubbles will either separate or coalesce [21,22].
We first consider aligned states, see §IV B 1, in which the two bubbles have similar centroid y−coordinates (y 1 ≈ y 2 ), and then offset states, see §IV B 2, where the two bubbles are on opposite sides of the rail (y 1 y 2 < 0).

Aligned Bubbles
Bubble pairs propagating from aligned initial configurations were prepared experimentally using the protocol outlined in appendix A 1. For both symmetric ('on rail') and asymmetric ('off rail') pairs of bubbles, when the larger bubble was initially placed behind the smaller where the leading bubble is larger than the trailing bubble: (a, b) aligned asymmetric bubbles overlapping the rail from one side, (c, d) aligned bubbles straddling the rail symmetrically about its centreline. The non-dimensional initial distance between the centroids of each bubble (which are indicated by crosses) is D(t = 0) = 2.10 (a), 1.77 (b), 2.39 (c) and 1.98 (d). The flow rate is Q = 0.04 (Q * = 106 ml/min), the total bubble area is A 1 + A 2 = 0.54 2 π and the bubble size ratio is A 1 /(A 1 + A 2 ) = 0.60. Each row of top view images shows the evolution of the system in terms of the non-dimensional time t = 2U * 0 t * /W * elapsed since flow initiation at t = 0. The dimensional time t * is indicated in the last snapshot of each time-sequence. (e) Time evolution of the non-dimensional distance D = 2D * /W * between the centroids of two asymmetric bubbles propagated from different initial separation distances. The two curves with time labels correspond to the time-sequences shown in (a,b). D c is the critical bubble distance delineating aggregation (blue curves) and separation (red curves).
one, the bubbles always coalesced. Hence, any hydrodynamic interactions were not sufficient to prevent the behaviour predicted from the single-bubble results, in which larger bubbles move faster. Figure 6 shows propagation experiments in which the larger bubble is initially leading. For a sufficiently large initial separation, the bubbles separate indefinitely (figures 6(a), asymmetric, and 6(c), symmetric), as predicted from the single-bubble results. For smaller initial separation distances, however, the two bubbles aggregate to form a compound bubble (figures 6(b), asymmetric, and 6(d), symmetric); a process that must be driven by hydrodynamic interactions between the two bubbles, which lead to an increase in the relative speed of the trailing bubble.
The hydrodynamic interactions arise from the changes in the bulk pressure field, which in the absence of the bubbles would decrease linearly along x. Single bubbles always propagate faster than the surrounding oil, as seen in §IV A, with accompanying local increases in pressure gradient. Hence, the fluid pressure at the front and rear of a bubble propagating in the channel is respectively higher and lower than the background pressure. The pressure perturbation decays with distance from the bubble, but increases with bubble volume because larger bubbles are faster. Consequently, when the smaller bubble is placed behind the larger, the net result of the perturbations due to both bubbles is that the trailing bubble experiences a lower local pressure near its tip, but a higher local pressure gradient causing the bubble to extend and narrow, see t = 5.8 in figure 6(b) and t = 8.3 in panel 6(d). The resulting changes in bubble shape cause an increase in speed of the trailing bubble and eventually it catches the bubble in front. The trailing bubble's speed continues to increase as the bubbles approach because the local pressure gradient increases, which further modifies the bubble shape. The interaction just described is generic and has been observed in two-bubble interactions in other confined systems [21,22]. The decay of the perturbations with distance means that if the bubbles are far enough apart the trailing bubble's speed does not increase sufficiently to allow it to catch the leading bubble.
The transition between the separation and aggregation outcomes observed in figures 6(ad) was investigated by performing successive experiments with a variety of initial bubble distances D(t = 0). Results are shown in figure 6(e), which presents the time evolution of the distance D(t) between two initially asymmetric bubbles similar to that of figures 6(a,b) for a variety of different initial D(t = 0). A transition between the two possible outcomes (aggregation in blue and separation in red) appears to occur for a threshold value of D(t = 0). All experiments, regardless of outcome, feature an initial decrease of D(t) for 0 < t < 1.4 which is associated with the rapid change in bubble shape following flow initiation; see e.g. figures 6(a,b) where bubbles are more slender at t = 1.4 than at t = 0. This is followed by a monotonic increase in the case of separation or a steepening decrease in the case of aggregation. The neighbouring red and blue curves which bound the range of initial separations where the transition occurs feature an approximately flat region after their initial decrease, indicating that bubbles initially travel with approximately constant separation. This suggests the existence of an unstable two-bubble steady mode of propagation where the two bubbles would neither separate or aggregate but always remain at the same critical distance D c from one another. We estimate D c to be the average between the values of D(t) for the (blue and red) curves adjacent to the threshold following initial decrease, i.e. D(t = 1.4) in figure 6(e). This unstable state is a so-called edge state that marks the boundary between bubble separation and bubble aggregation.
The evolution of two aligned bubbles in simulations of the theoretical model is very similar to that in the experiments. In figure 7, time-simulations calculated for bubbles of volume ratio 2:1, propagating at flow rate Q = 0.04 from different aligned initial conditions are presented as trajectories in a projection of the phase space plotting the bubble separation, D, against the offset of the centre of mass Y . Initial conditions with various initial global offsets Y (t = 0) and separation distances D(t = 0) are denoted by hollow markers labelled 'IC' and lead to either aggregation and then coalescence or separation of the two bubbles  depending on the value of D(t = 0), as shown in the inset snapshots of the final outcomes with solid-line bubble contours. Initial conditions with a global offset Y less than about 0.1 ultimately lead to one or two symmetric bubbles (y 1 = y 2 = 0) (blue curves) while initial conditions with a global offset Y larger than than about 0.2 ultimately lead to one or two asymmetric bubbles (red curves). Moreover, as suggested by the experimental results, we find that there are unstable steadily propagating states in the model that divide the different dynamical outcomes. The unstable steady states corresponding to the symmetric and asymmetric configurations are labelled S + 2 and AS + 2 , respectively, and were calculated using the method detailed in appendix §B 3.
For intermediate initial conditions, 0.1 ≤ Y (t = 0) ≤ 0.2, we often observe bubble break up leading to three bubbles, as either a transient part of the evolution or a permanent outcome. Figure 8 shows two examples with initial global offsets Y (t = 0) = 0.10 and 0.13 for the same initial separation distance D(t = 0) = 2.4 in panels (a) and (b) respectively. In panel (a), the bubbles oscillate until the smaller trailing bubble breaks up before finally coalescing to form a single bubble which later coalesces with the leading bubble, ultimately generating a single steady symmetric bubble. The initial oscillations are reminiscent of the unstable periodic orbit identified in [14] for single bubbles. In panel (b), the leading bubble is initially more asymmetric and breaks up as it 'hesitates' between an off-rail and on-rail configuration. However, here the bubbles do not coalesce and the final outcome is a three-bubble system with the larger leading bubble propagating faster than the two trailing bubbles which ultimately propagate steadily at the same speed on opposite sides of the rail.
These two examples leading to two radically different final outcomes illustrate the sensitivity of the system to the initial bubble offset. The second example also opens up the possibility of a stable two-bubble state featuring offset bubbles, which will be explored in §IV B 2. Motivated by the fact that, as explored in [9], a single bubble can break up into bubbles of arbitrary volume ratio, we now use parameter continuation to determine the effect of the volume ratio V r on two-bubble steadily propagating solutions. Figure 9 shows the bubble separation distance D c associated with the AS + 2 (red) and S + 2 (blue) edge states against the volume ratio V r for a fixed total bubble volume and a fixed flow rate Q = 0.04. The circular markers with error bars indicate experimental results and which are in reasonable agreement with numerical results. The agreement is generally within the experimental error for the symmetric states at smaller volume ratios, but the theoretical results consistently overpredict the separation distance for the asymmetric states, suggesting that the hydrodynamic interactions between two bubbles located near one edge of the rail are weaker in reality than in the model. The inset snapshots show the bubble configuration of the edge states at values of V r indicated by numbered markers on the solution branches. For both edge states, the separation distance D c decreases with increasing volume ratio and appears to converge to a finite value as V r → 1 (i.e. V 1 /V 2 → ∞), see insets 3 and 4, while increasing sharply as V r → 1/2 (i.e. V 1 /V 2 → 1), see insets 1 and 2. Linear stability results indicate that both branches are unstable with a single positive eigenvalue and that the least unstable eigenvalue approaches the imaginary axis as V r → 1/2, indicating that the state with two equal bubbles is neutrally stable, which is consistent with the results of [26], and as expected . Solution space as projected in the (V r , D c ) plane for the asymmetric AS + 2 (red) and symmetric S + 2 (blue) aligned edge states at fixed flow rate Q = 0.04 and for a fixed total bubble volume V total = 0.54 2 π. Dashed lines correspond to numerical results while data points with error bars correspond to experimental values for which the bubble sizes as quantified by their area A i instead of their volume V i (see equation (5)). Four snapshots of the numerical bubble shapes are shown for values of V r indicated by black circular markers on the numerical lines labelled by digits from 1 to 4.
because identical bubbles will travel at the same speed, assuming negligible hydrodynamic interactions. We note, however, that it is of course impossible in practice to have two bubbles of the exact same size in the experiments. Figure 10 shows a bifurcation diagram of the different aligned two-bubbles states calculated through parameter continuation, where the velocity U b of each state is plotted against Q for a constant volume ratio 2 : 1. Each solution branch is illustrated by at least one snapshot corresponding to a given value of Q indicated by a circular marker on the branch. When Q is larger than the transition flow rate Q s discussed in §IV A, there are two branches discussed in figure 7 featuring symmetric (S + 2 ) and asymmetric (AS + 2 ) bubbles. In this case the symmetric S + 2 branch only exists after a finite value of Q ≈ 0.03 and experiences a pitchfork bifurcation, after which the two bubbles are slightly asymmetric, see inset labelled 10. We note that the asymmetric AS + 2 branch persists for all values of Q > Q s calculated but, as Q approaches Q s from above, the branch terminates as the bubbles become increasingly further apart and the limits of the computational domain are reached. For completeness, we also include unstable symmetric (S − 2 ) and asymmetric (AS − 2 ) solutions calculated for Q < Q s where the leading bubble is now smaller than the trailing one. This is because the model predicts that for single steady bubbles, smaller bubbles propagate faster. However, as discussed in §IV A, numerical results for Q < Q s do not reflect the experiments since there is no such transition flow rate in the experiments. The similarities between the two-bubble and single-bubble bifurcation diagrams presented in figures 10 and 5 will be discussed in section IV C.

Offset bubbles
We now consider offset bubble-pair configurations in which the two bubbles are initially positioned on opposite sides of the rail in an 'Up-Down' (later denoted UD) configuration. The experimental protocol to prepare these configurations is outlined in appendix A 2.
In figures 11(a-d), we present experimental time-sequences for two bubbles of area ratio 2:1 propagating at flow rate Q = 0.04 from different initial conditions. As in the case of aligned bubbles, if the larger bubble leads and the bubbles are initially well separated, figure  11(a), there is no significant hydrodynamic interaction and the bubbles separate indefinitely, remaining on their respective sides of the rail. As the distance between the bubbles decreases, then the hydrodynamic interaction between the bubbles is such that the trailing bubble migrates across the rail as it responds to the locally increased pressure gradient introduced by the leading bubble, see t = 4.0 and t = 3.6 in figures 11(b) and 11(c) respectively. Once the bubbles are on the same side of the rail, the system is in the aligned configuration, see §IV B 1. The two bubbles separate indefinitely in figure 11(b) and aggregate in figure 11(c) owing to the different values of D after bubble migration.
In figure 11(d), we consider the reverse initial configuration where the larger bubble is initially trailing. We observe that at first the trailing bubble propagates faster, as it would in the absence of hydrodynamic interaction, and hence the distance D between the two bubbles decreases with time. As the bubbles approach, the trailing bubble starts to migrate Bubbles have a total area A total = 0.54 2 π and a 2:1 area ratio. The time labelling is the same as in figure 6. (e): Corresponding numerical time-simulations presented as trajectories in a phase-plane projection using the coordinated (x 1 , y 1 ) of the larger bubble. The flow rate is Q = 0.04 and bubbles have total volume V total = 0.54 2 π with a 2:1 volume ratio. Initial conditions are denoted by hollow circles, steady states by solid triangles and an star represents coalescence. Inset panels with label 'IC' show three typical initial bubble configurations, inset panels with dashed bubble contours correspond to four steady states and inset panels with solid bubble contours and coloured outlines correspond to four different final outcomes for selected trajectories.
over the rail, which causes it to slow down owing to the reduced viscous dissipation resulting from a smaller volume of fluid being displaced. In fact, the leading bubble also migrates further over the rail and so both bubbles slow down before ultimately reaching a steadily propagating state with a constant separation distance D c , see figure 11(d) at t = 25.3. This suggests the existence of a stable steady state of the two-bubble system, which is confirmed by numerical simulations.
Corresponding numerical time-simulations are shown in figure 11(e) and presented as trajectories in a (x 1 , y 1 ) projection of the phase space, where (x 1 , y 1 ) are the coordinates of the centroid of the largest bubble. We recall that the x-coordinate of the centre of mass of the system is constrained to be zero throughout calculations. The behaviour is similar to that observed in the experiments: in the case of an initially larger leading bubble (x 1 > 0), the final outcome of the system switches from offset separation (orange trajectories) to aligned separation (blue trajectories) and finally to bubble coalescence (black trajectories) as the initial distance between the two bubbles is decreased. Figure 11(e) shows that when the smaller trailing bubble crosses the rail, the AS + 2 unstable steady state discussed in §IV B 1 acts as an edge state delineating aligned separation from coalescence outcomes, as blue and black trajectories approaching it from different sides are deflected in different directions. In the case of an initially larger trailing bubble (x 1 < 0), all (red) trajectories in figure 11(e) converge towards a stable steadily propagating state labelled UD 1 2 , irrespective of the initial distance between the bubbles and consistent with experimental observations. We identify two further unstable steadily propagating states in figure 11(e), labelled UD 0 2 and UD 2 2 . The red trajectory starting from the initial condition furthest to the right (x 1 (t = 0) = −0.1) is first attracted towards the weakly unstable steady state UD 0 2 . A transient bubble configuration extracted from such a time-simulation is used as an initial guess for calculating the UD 0 2 state. This state plays a role in transient dynamics involving two bubbles that are almost on top of each other, which occurs for example after breakup of a single bubble like in figure 2 and was previously identified in [9], along with another unstable state labelled 'Barrier' state that is not discussed in the present paper. In contrast, the state UD 2 2 appears to have no influence on the dynamics. We now use parameter continuation to examine the behaviour of the three offset twobubble steady states UD 0 2 , UD 1 2 and UD 2 2 as we vary the flow rate Q and the bubble volume ratio V r . The two associated bifurcation diagrams are presented in figure 12 where in panel (a) we present the velocity U b associated with each state against Q for a constant volume ratio V r = 2/3 and where in panel (b) we present the bubble separation distance D c against V r for a constant flow rate Q = 0.04, similar to figures 10 and 9 respectively for the aligned bubble states. Solid and dashed lines indicate stable and unstable solutions respectively. Each branch is illustrated by at least one inset snapshot corresponding to a flow rate indicated by a circular marker on the branch. Experimental data points corresponding to measurements on the UD 1 2 stable state are also shown. Figure 12(a) shows that the stable UD 1 2 and unstable UD 2 2 solutions are connected through a limit point denoted LP1. Unlike UD 0 2 , the UD 2 2 unstable state appears to have no influence on the transient dynamics of the system according to figure 11(e). The limit point LP1 occurs at a flow rate that coincides with the transition flow rate Q s ≈ 0.0096, discussed in figure 5, where the mathematical model predicts that the smaller bubble propagates faster than the larger one in their respective AS1 single-bubble mode of propagation. This is consistent with the fact that a smaller leading bubble would then propagate faster in numerical time-simulations for Q < Q s so that the two bubbles would separate without  reaching a steady state. Surprisingly, although we know from §IV A that there is no evidence for Q s in the experiments, a threshold, consistent with a limit point, of the UD 1 2 state is also found experimentally at a critical flow rate Q LP1 = 0.0040 ± 0.0002 for a bubble area ratio A r ≈ 1/3. A representative time-evolution observed experimentally for Q < Q LP 1 is shown in figure 13(a). In this case, the migration of the trailing bubble over the rail induced by the bubble interaction does not cause a sufficient speed reduction to reach a steadily propagating two-bubble state. Instead, the larger bubble passes over the smaller bubble so that both bubbles ultimately propagate steadily and separate indefinitely. This low flow rate scenario is not captured using our model. We also note that the UD 0 2 branch in figure 12(a) is distinct from the two others and could not be calculated numerically below a flow rate indicated by the label '1' at which the smaller bubble touches the side wall of the channel, see associated inset snapshot. Figure 12(b) shows that the stable UD 1 2 and unstable UD 0 2 solutions, which were disconnected in the (Q, U b ) projection of figure 12(a), are in fact connected through a limit point denoted LP2 under variations in volume ratio. This means that there is a critical volume ratio below which the stable UD 1 2 state does not exist, which we find to be V r ≈ 0.21 at Q = 0.04. The existence of such a critical volume ratio is supported by our experiments, in which no stable state is found for an area ratio A r < 0.256 ± 0.007 at the same flow rate. A representative experimental time-evolution at A r = 0.247 is shown in figure 13(b) and equivalent numerical simulations are qualitatively similar. The overall dynamics are essentially the same as those at low Q below LP1 in the experiments: for values of the volume ratio below LP2 the interaction between the bubbles does not reduce the speed of the trailing bubble sufficiently to establish a steadily propagating two-bubble state. The only qualitative difference to the time evolution shown in figure 13(a) is that the smaller bubble migrates over the rail once the larger bubble has moved ahead. Figure 12(b) also shows that, like in figure 9 for aligned bubbles, the distance D c between the bubbles increases when approaching the limit of two bubbles of equal sizes (V r = 1/2) for the UD 1 2 and UD 2 2 branches. The chosen fixed length of the computational domain means that the UD 1 2 branch solution could not be calculated close to that limit. By contrast, the UD 0 2 solution could be calculated up to the V r = 1/2 limit where it features two identical bubbles on either side of the rail propagating at the same x-position with opposite offsets y 2 = −y 1 .

C. Comparison of single-and two-bubble solution structures
There is a striking similarity between the solution structure for single bubbles presented in figure 5 and discussed in our previous papers [9,14], and the solution structure for two bubbles presented in figures 7 and 12. For ease of comparison, figure 14 shows a direct comparison between the solutions in the (Q, U b ) plane for a single bubble (coloured lines) and a two-bubble system (black lines). The total bubble volume in the two-bubble system is V total = π0.54 2 and the bubble volumes have a ratio 2:1. The bubble volume in the singlebubble system is the same as the volume of the larger bubble of the two-bubble system, i.e. V = 2/3V total . A selection of the branches are illustrated by snapshots at given flow rates indicated by a circular marker.
For Q > Q s the two-bubble AS + 2 and UD 2 2 solution branches overlap and closely match the single-bubble asymmetric state, labelled AS + 1 . Furthermore the two-bubble symmetric S + 2 solution branch is almost indistinguishable from the symmetric one-bubble state, labelled S + 1 . For Q < Q s there is no single-bubble state corresponding to AS − 2 , but the symmetric two-bubble S − 2 solution branch has a similar structure to the single-bubble S − 1 branch, albeit without the excellent quantitative agreement found in Q > Q s case. We examined variations in the single-bubble volume and found that the closest match between the single and twobubble solution branches for Q > Q s occurs when the single bubble volume is chosen to be equal to that of the leading bubble of the two-bubble system, as shown in figure 14, which indicates that the leading bubble (not necessarily the fastest) sets the speed of the bubble pair. This comparison is particularly evident when comparing two-bubble branches with a larger leading bubble (see the AS + 2 and S + 2 and their single-bubble counterparts) but, as can be seen from the figure, the stable two bubble UD 1 2 state is slightly slower than the corresponding AS1 state corresponding to leading (smaller) bubble.

V. DISCUSSION
In this paper we studied the propagation of two bubbles through a geometrically perturbed Hele-Shaw channel under constant flow rate: the simplest configuration that introduces bubble-bubble interactions in the system. The study is of fundamental interest in a variety of applications, but our initial motivation was to conduct a controlled investigation of the post break-up dynamics of the single bubbles that we studied in our previous paper [9]. When a single bubble breaks up, the relative positions and volumes of the resulting multiple bubbles are very sensitive to perturbations in the system and are extremely difficult to control. In this paper, we fixed the sum of the two-bubble volumes, but varied the relative sizes of each bubble to simulate different break-up configurations.
We find that the general behaviour of the two-bubble system falls into four different longterm outcomes: (i) indefinite separation of the two bubbles; (ii) aggregation and coalescence of the two bubbles; (iii) a steadily propagating two-bubble state and (iv) break up to form a larger number of bubbles with potential future break-up and aggregation/coalescence events. We demonstrate that, as in the case of the single bubble, the overall dynamics is orchestrated by both stable and unstable steadily-propagating one-and two-bubble states. The ranges of existence and stability of the states depend on the flow rate and relative sizes of the two bubbles. The general behaviour of the system is qualitatively described by a depth-averaged theory, provided that the flow rate is sufficiently large, Q > Q s .
The existence of steadily-propagating two-bubble states is striking because they have not been in observed in other confined systems [21,22], in which the multiple bubbles will always separate or aggregate. There are large regions of overlap between particular single-and twobubble steadily propagating states in the relationship between bubble speed and flow rate for a fixed volume ratio between the two bubbles, see figure 14. The overlap regions suggest that these single-and two-bubble states are closely related. In particular, over a wide range of flow rates, the unstable asymmetric two-bubble states UD 2 2 and AS + 2 have the same speed as the stable asymmetric single bubble AS 1 ; and the unstable symmetric two-bubble state S 3 has the same speed as the stable symmetric single-bubble S + 1 . In these comparisons, the single bubble always has the same volume as the larger of the two bubbles in the two-bubble state, indicating that the two-bubble state moves at the speed of the leading bubble and that the leading bubble is not significantly affected by the interaction. In other words, the leading bubble is driving the dynamics and the trailing bubble is carried along with it. This is only possible because the presence of the rail allows the trailing bubble to experience different local geometric confinements depending on its lateral position within the channel.
Having established the existence of two-bubble states, we can extend the methods used in this paper to construct a variety of multiple bubble states and the number of possible states increases dramatically with the number of bubbles, in line with the number of permutations of increasing numbers of discrete objects. Experimental confirmation of the existence of what appear to be stable three-and four-bubble steadily propagating states is given in figure 15. The existence of stable and unstable n-bubble steadily propagating states will have a potential influence on the dynamics of bubble trains in confined systems [2,3] in the presence of imperfections in both geometry and bubble volume. The ranges of existence of the multiple-bubble solutions and their sensitivity to perturbations will be pursued in a future investigation.
Finally, although there is not a direct equivalent of the UD 0 2 state in the single-bubble system, the UD 0 2 state is the only two-bubble state that persists if the height of the rail, h, is decreased to 0. All of the other steady states cease to exist because their separation distances increases as h → 0, but the UD 0 2 state barely changes and exists exactly when h = 0. The question of whether a stable two-bubble steady state exists when h = 0 is an open question but certainly if one does exist it is is unlikely to be related to the steady states constructed here. Indeed, numerical IVP calculations confirm that a two-bubble steady state does not exist in the same form as the UD 1 2 solution (smaller bubble head on opposite sides) and that the dynamics are incredibly slow, indicating neutral stability. This exploration of two-bubble steady states in the experiment will form part of a future investigation. is established, for this IC, either D s or D c is updated to the value of D edge (t = 0), as appropriate, and a new value of D edge (t = 0) is chosen from (B1).
This interval-bisection procedure is repeated so that the initial bubble distance D edge converges to a value so that D(t) → D crit corresponding to the unstable steady state. In each simulation, the volume and the initial offset and shape (slender ellipse) of each bubble is kept constant while only varying the initial distance between the two bubbles. The final dynamical outcome is determined when either the minimum distance between the two bubble contours gets smaller than a cutoff value D(t) < D min , which is small enough to ensure that the bubbles will coalesce, or when the centroid-to-centroid distance D(t) gets larger than a cutoff value D(t) > D max which ensures that the bubbles will separate indefinitely. We use values of D min = 0.01, D max = 3 for the results of this paper.
The convergence criteria for the interval-bisection procedure is that the bubble remain within a small distance, ε, of each other for a large time, T , i.e. |D(t) − D edge (t = 0)| < ε, ∀t < T end . (B2) Once this condition has been satisfied we solve the steady governing equations to get the unstable steady state. In the results presented here find that ε = 0.1 and T = 20 is sufficient to ensure that we converged to a steady state.