1. Introduction
Identifying invariant objects (steady states, periodic orbits, invariant tori etc.) of high-dimensional, 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 identifying these objects is to perform a number of initial-value problems (IVPs), 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 (Kerswell, Pringle & Willis Reference Kerswell, Pringle and Willis2014) or parameter continuation (Kuznetsov Reference Kuznetsov1998; Net & Sánchez Reference Net and Sánchez2015). 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 (Kerswell Reference Kerswell2005; Schneider, Eckhardt & Yorke Reference Schneider, Eckhardt and Yorke2007; Schneider et al. Reference Schneider, Gibson, Lagha, De Lillo and Eckhardt2008; Kawahara, Uhlmann & Van Veen Reference Kawahara, Uhlmann and Van Veen2012), 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 (Gaillard et al. Reference Gaillard, Keeler, Thompson, Hazel and Juel2021), 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 was 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 perturbed Hele-Shaw channel. 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.
Experimental time snapshots of the evolution of a bubble in a perturbed Hele-Shaw channel, as viewed from above: (a) 
$Q=0.02$. In this case, post-breakup, the smaller bubble crosses the rail and eventually coalesces with the larger bubble in an asymmetric configuration; (b) 
$Q=0.03$. In this case, post-breakup, the two bubbles separate indefinitely on either side of the rail.
A feature of this system is that the topology of the system changes when a bubble breaks up or when 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,\ldots, n-1$-bubble systems in such a way that a hierarchy of 
$1,2,\ldots, n$-bubble states can be constructed from smaller bubble systems. The broad phenomenon of ‘lower-order’ 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 Drazin & Johnson Reference Drazin and Johnson1989) and nonlinear optics (see, for example Akhmediev & Ankiewicz Reference Akhmediev and Ankiewicz2005), spatially localised states in convection systems (see, for example Mercader et al. Reference Mercader, Batiste, Alonso and Knobloch2010) and oscillons in granular particulate flow (see, for example Umbanhowar, Melo & Swinney Reference Umbanhowar, Melo and Swinney1996). 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 non-trivial and this system represents a rather different example of interacting localised states, from those 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 Taylor & Saffman Reference Taylor and Saffman1959; Tanveer Reference Tanveer1987; Tanveer & Saffman Reference Tanveer and Saffman1987; Khalid, McDonald & Vanden-Broeck Reference Khalid, McDonald and Vanden-Broeck2015; Green, Lustri & McCue Reference Green, Lustri and McCue2017; Lustri, Green & McCue Reference Lustri, Green and McCue2020) and experiments (see, for example Maxworthy Reference Maxworthy1986; Kopf-Sill & Homsy Reference Kopf-Sill and Homsy1988; Wang et al. Reference Wang, Klaasen, Degrève, Blanpain and Verhaeghe2014; Zhang et al. Reference Zhang, Mines, Lee and Jung2016; Madec Reference Madec2021; Sirino, Mancilla & Morales Reference Sirino, Mancilla and Morales2021). A key result from these studies is that there is only one stable single-bubble solution for a channel of constant rectangular cross-section. However, 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. (Reference Franco-Gómez, Thompson, Hazel and Juel2017, Reference Franco-Gómez, Thompson, Hazel and Juel2018); Keeler et al. (Reference Keeler, Thompson, Lemoult, Juel and Hazel2019) and Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021). The solution branches interact in a highly non-trivial 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 (Keeler et al. Reference Keeler, Thompson, Lemoult, Juel and Hazel2019; Gaillard et al. Reference Gaillard, Keeler, Thompson, Hazel and Juel2021).
Although multi-bubble steady propagation modes in a Hele-Shaw channel have been studied in unperturbed channels, see, for example Crowdy (Reference Crowdy2009), Green et al. (Reference Green, Lustri and McCue2017), Vasconcelos (Reference Vasconcelos2015) and Lustri et al. (Reference Lustri, Green and McCue2020), 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, (Pumir & Aref Reference Pumir and Aref1987; Sarig, Starosvetsky & Gat Reference Sarig, Starosvetsky and Gat2016; Green Reference Green2018), 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, are also neutrally stable, whilst two rows of bubbles which are located asymmetrically about the horizontal centreline are unstable, see Pumir & Aref (Reference Pumir and Aref1987). 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 (Maxworthy Reference Maxworthy1986; Rohilla & Das Reference Rohilla and Das2020; Madec Reference Madec2021).
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 is 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 §§ 2 and 3 we present the experimental and numerical methods used to investigate the dynamics of the system. In § 4.1 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 (§ 4.2.1), and ‘offset’ states where the bubbles are staggered on either side of the rail (§ 4.2.2). Next, in § 4.3, we compare the solution structures of the two-bubble system with the single-bubble system before we conclude with a discussion of the implications of our results for 
$n$-bubble systems (§ 5).
2. 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 Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021). 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 \pm 0.1$ mm and height 
$H^* = 1.00 \pm 0.01$ mm, with an aspect ratio 
$\alpha = W^*/H^* = 40$. A centred rail of width 
$w^* = 10.0 \pm 0.1$ mm and thickness 
$h^* = 24 \pm 1\,\mathrm {\mu }$m consisted of a translucent adhesive tape strip bonded to the bottom glass plate, see figure 3(b).
(a) Schematic of the experimental set-up and (b) experimental channel in cross-sectional view.
The channel was filled with silicone oil (Basildon Chemicals Ltd) of dynamic viscosity 
$\mu = 0.019$ Pa s, density 
$\rho = 951\,{\rm kg}\,{\rm m}^{-3}$ and surface tension 
$\sigma = 21\,{\rm mN}\,{\rm m}^{-1}$ at the laboratory temperature of 
$21\pm 1\,^{\circ }{\rm C}$, and connected to oil reservoirs through inlet and outlet ports located at the upstream and downstream ends of the channel, respectively (see figure 3a). 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. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021) 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 3a). 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.
Sketch of the non-dimensional computational domain, which is in a frame of reference centred on the overall centre of mass that moves with speed 
$U_{{b}}$. The horizontal domain is truncated at a value of 
$x = L$ (typically 
$L=6$ in our calculations) and the boundaries of the rectangular rail are marked using dotted horizontal lines. The fluid domain is denoted 
$\varOmega$, and the two air bubble boundaries are denoted 
$\varGamma _1$ and 
$\varGamma _2$. The centroids of the bubbles are denoted 
$(\bar {x}_i,\bar {y}_i)$. The inset diagram shows the smoothed obstacle in the 
$(y,z)$ plane (solid line), and a rectangular step (dashed line), corresponding to the limit 
$s\to \infty$.
Bubbles were propagated from their initial configuration at a constant dimensionless flow rate 
$Q = \mu U_0^* / \sigma$ 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_1 + A_2 = 0.54^2 {\rm \pi}$, which is the size of the single bubbles used in Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021), 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 leading position.
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 up to 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. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021) 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$). These lubricating oil films are a reason why the quantitative prediction of the model (discussed below) does not precisely match the experiment, although the qualitative agreement is excellent, as shown extensively in Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021).
3. 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 McLean & Saffman (Reference McLean and Saffman1980) to account for a non-uniform channel height and has been used extensively in studies of the propagation of a semi-infinite air finger (Thompson, Juel & Hazel Reference Thompson, Juel and Hazel2014; Franco-Gómez et al. Reference Franco-Gómez, Thompson, Hazel and Juel2016), single closed air bubbles (Franco-Gómez et al. Reference Franco-Gómez, Thompson, Hazel and Juel2017, Reference Franco-Gómez, Thompson, Hazel and Juel2018; Keeler et al. Reference Keeler, Thompson, Lemoult, Juel and Hazel2019) and most recently single and multiple air bubbles (Gaillard et al. Reference Gaillard, Keeler, Thompson, Hazel and Juel2021). We use the model to compute steady states of the system, calculate their linear stability and perform numerical time-dependent 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 non-dimensional 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
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} (b^3 \boldsymbol{\nabla} p) = 0\quad (x,y)\in \varOmega, \end{equation}
where the mobility 
$b(y)$ represents the variable depth of the channel, modelled as a smoothed tanh profile
\begin{equation} b(y) = 1 - \tfrac{1}{2}h\left[\tanh(s(y + w)) - \tanh(s(y - w))\right], \end{equation}
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, which in the limit as 
$s\to \infty$, corresponds to a rectangular step; see the inset in figure 4. We use 
$h=0.024$ and 
$w=0.25$ consistent with experiments and we choose 
$s=40$ (Thompson et al. Reference Thompson, Juel and Hazel2014). We impose no-penetration conditions on the upper and lower walls, which yield 
$p_y = 0$ on 
$y=\pm 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 
$\boldsymbol {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 
$\boldsymbol {R}_i$ (where 
$i=1,2,\ldots$ indicates the 
$i$th bubble in decreasing size order) is given by
\begin{equation} \frac{\partial\,\boldsymbol{R}_{i}}{\partial\,t}\boldsymbol{\cdot} \boldsymbol{n}_i + \boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{n}_i + b^2\boldsymbol{\nabla} p\boldsymbol{\cdot} \boldsymbol{n}_i = 0, \end{equation}
where 
$\boldsymbol {n}_i$ is the unit normal vector directed away from the 
$i$th bubble and 
$\boldsymbol {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 is based on a Young–Laplace law (see, for example, De Gennes Reference De Gennes2004), where the pressure jump across the bubble interface is proportional to the mean curvature. In our system this can be written as
\begin{equation} p_i - p = \frac{1}{3\alpha Q}\left(\frac{\kappa}{\alpha} + \frac{1}{b(y)}\right), \end{equation}
where 
$\kappa$ 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 (which assumes that the interface is a semi-circle occupying the entire channel height in the 
$(y,z)$ plane). The aspect ratio is set to 
$\alpha = 40$. We expect the qualitative features of the bifurcation structure to remain similar over a broad range of 
$\alpha$, as found for air fingers (Franco-Gómez et al. Reference Franco-Gómez, Thompson, Hazel and Juel2016). Franco-Gómez et al. (Reference Franco-Gómez, Thompson, Hazel and Juel2016) also found excellent agreement between the model and experimental results when 
$\alpha > 30$. 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$ and area 
$A_i$, are defined via
\begin{equation} V_i = \int_{\varGamma_i} b(y) \, \mathrm{d} x \, \mathrm{d} y= \int_{\partial \varGamma_i} x b(y) \, \mathrm{d} y, \quad A_i= \int_{\varGamma_i} \, \mathrm{d} x \, \mathrm{d} y, \end{equation}
where 
$\varGamma _i$ is the interior of bubble 
$i$ and 
$\partial \varGamma _i$ its bounding curve. We describe two-bubble systems by the ratio of the largest bubble's volume to the total volume, i.e. 
$V_{{r}} = V_1/(V_1 + V_2)$, and the overall vertical centre of mass of the system, defined as 
$\bar {Y} = (V_1\bar {y}_1 + V_2\bar {y}_2)/(V_1 + V_2)$, where 
$(\bar {x}_i,\bar {y}_i)$ are the centroid coordinates of the 
$i{\mathrm {th}}$ bubble. The total bubble volume 
$V_1+V_2$ is set to 
$0.54^2{\rm \pi}$ unless specified otherwise, consistent with the experiments of Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021). 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.4 %. 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 § 2.
We solve the system of equations, (3.1)–(3.5a,b) on the domain shown in figure 4 to determine 
$p$, 
$\boldsymbol {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 using the open-source oomph-lib library (Heil & Hazel Reference Heil and Hazel2006). During time-dependent calculations a second-order backwards differentiation Euler method is used with a typical time step of 
$\Delta t = 0.01$. For more details of this method, applied to a similar problem, see Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021).
When performing time-dependent simulations, the initial shape of each bubble was chosen to be an ellipse with contour coordinates
\begin{equation} \boldsymbol{R}_i(t=0) = (\bar{x}_{i}(0) + \ell_i\cos\theta , \, \bar{y}_{i}(0) + d_i\sin\theta),\quad 0\leq\theta<2{\rm \pi}. \end{equation}
In all the numerical time-dependent simulations presented in this paper, the volume ratio is 
$V_r= 2/3$ 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 
$\ell _i = A_i/{\rm \pi} d_i$ were set to ensure the prescribed volumes.
To account for topology changes that may occur, such as bubble breakup and coalescence events, we use a procedure detailed in Appendix B.1. 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-dependent 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.2. 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.
The one- and two-bubble steady states described in §§ 4.1 and 4.2 are labelled in the form ‘
$nX_m$’, where 
$n$ corresponds to the number of bubbles in the system, 
$X$ is a Latin character which is either ‘A/S’, corresponding to whether the state is asymmetric (
$\bar {Y}\neq 0$) or symmetric (
$\bar {Y} = 0$), respectively. For two-bubble states, this identifier could also be ‘F’, which corresponds to an ‘offset’ state (two bubbles on either side of the occlusion so that 
$\bar {y}_1\bar {y}_2 < 0$). The subscript 
$m$ is a label that corresponds to the order in which the branches appear, as 
$Q$ increases from zero. We also use the subscript to distinguish between branches that, although connected, have different linear stability properties (cf. 1S
$_2$ and 1S
$_3$ in the section below).
4. Results
4.1. 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. For our region of interest, there are three distinct solutions: a stable asymmetric bubble, 1A
$_1$, that exists for all flow rates; an unstable symmetric double-tipped bubble, 1S
$_1$, that exists for small flow rates; and an alternative symmetric bubble, 1S
$_2$/1S
$_3$, that exists for larger flow rates and is unstable/stable, respectively. 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.
Single-bubble solution space. The velocity 
$U_{b}$ is plotted as a function of dimensionless flow rate 
$Q$ for (a) the different solution branches of the theoretical model and (b,c) for the stable steady modes of propagation observed in experiments, where (c) focuses on the lowest flow rates, for different bubble volumes/areas. The experimental data presented in (b) are reproduced from Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021) with permission. The inset panels in (a) correspond to bubble profiles specified by solid circular markers on the branches. On (a), the bubbles have volumes 
$V={\rm \pi} 0.54^2,{\rm \pi} 0.46^2,{\rm \pi} 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 1S
$_2$ state and the solid circular markers in (b,c) correspond to the asymmetric 1A
$_1$ 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, in the experiments, 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 Taylor & Saffman (Reference Taylor and Saffman1959), 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, (3.3). Related results have also been found for buoyant rise of bubbles in Hele-Shaw cells (Maxworthy Reference Maxworthy1986) 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 Maxworthy (Reference Maxworthy1986) 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) 
$\to 0$’ (Maxworthy Reference Maxworthy1986, p. 108). We do not have a specific explanation for this discrepancy, but it is most likely to be a violation of the modelling assumptions used for the boundary conditions at the bubble interface. For example, in our model, the dynamic condition (3.4) only considers normal stresses in the 
$(x,y)$ plane, lateral stresses are ignored, which can cause smaller bubbles (in an unperturbed cell without side walls) to slow down (Nagel Reference Nagel2014). However, it is not clear how lateral stress effects interact with bubble shape deformation, which can be more pronounced at larger 
$Q$, nor with the channel sidewalls or rail. In addition, we assume that the transverse curvature is constant in the transverse 
$(y,z)$ plane, corresponding to a semi-circle perfectly wetting the lower and upper walls, which is not strictly valid when considering small bubbles. 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.
4.2. Two-bubble systems
From the results for single bubbles, § 4.1, 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 (Maxworthy Reference Maxworthy1986; Madec Reference Madec2021).
We first consider aligned states, see § 4.2.1, in which the two bubbles are on the same side of the rail, 
$\bar{y}_1\bar{y}_2 \geq 0$, and then ‘offset’ states, see § 4.2.2, where the two bubbles are on opposite sides of the rail (
$\bar {y}_1 \bar {y}_2 < 0$).
4.2.1. 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 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.
Time sequences of two-bubble evolution for different initial bubble-pair configurations, 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\,{\rm ml}\,{\rm min}^{-1}$), the total bubble area is 
$A_1 + A_2 = 0.54^2 {\rm \pi}$ 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 = 2 U_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).
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 § 4.1, 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 figure 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 (Maxworthy Reference Maxworthy1986; Madec Reference Madec2021). The decay of the perturbations with distance means that if the bubbles are far enough apart, like in figure 6(a,c), 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 figure 6(a–d) 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 figure 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. figure 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-dependent simulations calculated for bubbles of volume ratio 
$V_{{r}} = 2/3$, 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$ (distance between the two centroids), against the offset of the centre of mass 
$\bar {Y}$. Initial conditions with various initial global offsets 
$\bar {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 
$\bar {Y}$ less than approximately 
$0.1$ ultimately lead to one or two symmetric bubbles 
$(\bar {y}_1 = \bar {y}_2 = 0)$ (blue curves) while initial conditions with a global offset 
$\bar {Y}$ larger than approximately 
$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 2S
$_2$ and 2A
$_2$, respectively, and were calculated using the method detailed in Appendix B.2.
Time-dependent calculations when 
$Q = 0.04>Q_{{s}}$ so that larger bubbles are faster, 
$V_{{r}} = 2/3, V_1+V_2 = {\rm \pi}0.54^2$. The lines are trajectories in the projected phase plane 
$(D,\bar {Y})$. Hollow circles indicate initial conditions (IC) and the arrows on the lines indicate increasing values of time, 
$t$. The triangles indicate the unstable steady states, shown as dashed bubble contours in the insets. Insets with solid line contours indicate the shapes of the bubbles at the stated time and correspond to solid markers on the trajectories. The dashed lines indicate a region where time trajectories feature the breakup of at least one of the two bubbles.
For intermediate initial conditions, 
$0.1 \le \bar {Y}(t=0) \le 0.2$, we often observe bubble breakup 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 
$\bar {Y}(t=0) = 0.10$ and 
$0.13$ for the same initial separation distance 
$D(t=0) = 2.4$ in panels (a,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 Keeler et al. (Reference Keeler, Thompson, Lemoult, Juel and Hazel2019) 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 § 4.2.2.
Numerical time snapshots of the evolution of the system at times indicated in each panel, starting from initial conditions in the breakup zone in figure 7 at 
$Q=0.04$. Initial conditions are (a) 
$D=2.4,\bar {y}_1 = \bar {y}_2 = 0.10$ and (b) 
$D=2.4,\bar {y}_1 = \bar {y}_2 = 0.13$. In (b), the two trailing bubbles ultimately propagate at the same velocity, as indicated by the dashed box.
Motivated by the fact that, as explored in Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021), 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 2A
$_2$ (red) and 2S
$_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 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 over-predict 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}}\to 1$ (i.e. 
$V_1/V_2 \to \infty$), see insets 3 and 4, while increasing sharply as 
$V_{{r}}\to 1/2$ (i.e. 
$V_1/V_2 \to 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}}\to 1/2$, indicating that the state with two equal bubbles is neutrally stable, which is consistent with the results of Pumir & Aref (Reference Pumir and Aref1987), and as expected 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.
Solution space as projected in the 
$(V_{{r}},D_{{c}})$ plane for the asymmetric 2A
$_2$ (red) and symmetric 2S
$_2$ (blue) aligned edge states at fixed flow rate 
$Q=0.04$ and for a fixed total bubble volume 
$V_1+V_2=0.54^2{\rm \pi}$. Dashed lines correspond to numerical results while data points with error bars correspond to experimental values for which the bubble sizes are quantified by their area 
$A_i$ instead of their volume 
$V_i$ (see (3.5a,b)). 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.
Figure 10 shows a bifurcation diagram of the different aligned two-bubble states calculated through parameter continuation, where the velocity 
$U_{{b}}$ of each state is plotted against 
$Q$ for a constant volume ratio 
$V_{{r}} = 2/3$. 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 § 4.1, there are two branches discussed in figure 7 featuring symmetric (2S
$_2$) and asymmetric (2A
$_2$) bubbles. In this case the symmetric 2S
$_2$ branch only exists after a finite value of 
$Q\approx 0.03$ and experiences a pitchfork bifurcation, after which the two bubbles are slightly asymmetric, see inset labelled 10. We note that the asymmetric 2A
$_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 (2S
$_1$) and asymmetric (2A
$_1$) 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 § 4.1, 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 § 4.3.
Solution space for aligned bubbles, as projected in the 
$(Q,U_{{b}})$ plane for a fixed bubble volume 
$V_1+V_2=0.54^2{\rm \pi}$ and volume ratio 
$V_{{r}} = 2/3$. The circular markers indicate solutions that are shown in the inset panels. Each inset has its own numerical label. Branches are dashed as the solutions are unstable and the dotted vertical line marks the position of 
$Q_{{s}}\approx 0.0096$.
4.2.2. 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 ‘offset’ configuration. The experimental protocol to prepare these configurations is outlined in Appendix A.2.
In figure 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 § 4.2.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.
(a–d) Experimental time snapshots of two bubbles propagating at flow rate 
$Q=0.04$ from initial offset configurations shown at 
$t=0$. The larger bubble is initially leading in (a–c), with decreasing separation distance from (a) to (c), and trailing in (d). Bubbles have a total area 
$A_{1}+A_{2} = 0.54^2{\rm \pi}$ and a 
$2:1$ area ratio. The time labelling is the same as in figure 6. (e) Corresponding numerical time-dependent simulations presented as trajectories in a phase-plane projection using the coordinated 
$(\bar {x}_1, \bar {y}_1)$ of the larger bubble. The flow rate is 
$Q=0.04$ and bubbles have total volume 
$V_1+V_2 = 0.54^2{\rm \pi}$ with 
$V_{{r}} = 2/3$. Initial conditions are denoted by hollow circles, steady states by solid triangles and a 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.
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 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-dependent simulations are shown in figure 11(e) and presented as trajectories in a 
$(\bar {x}_1,\bar {y}_1)$ projection of the phase space, where 
$(\bar {x}_1,\bar {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 (
$\bar {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 2A
$_2$ unstable steady state discussed in § 4.2.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 (
$\bar {x}_1<0$), all (red) trajectories in figure 11(e) converge towards a stable steadily propagating state labelled 2F
$_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 2F
$_1$ and 2F
$_3$. The red trajectory starting from the initial condition furthest to the right (
$\bar {x}_1(t=0) = -0.1$) is first attracted towards the weakly unstable steady state 2F
$_1$. A transient bubble configuration extracted from such a time-dependent simulation is used as an initial guess for calculating the 2F
$_1$ 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 Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021), along with another unstable state labelled ‘Barrier’ state that is not discussed in the present paper. In contrast, the state 2F
$_3$ appears to have no influence on the dynamics.
We now use parameter continuation to examine the behaviour of the three offset two-bubble steady states 2F
$_1$, 2F
$_2$ and 2F
$_3$ 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 2F
$_2$ stable state are also shown as solid, circular markers.
Solution space for the offset steady states. Solid/dashed lines indicate stable/unstable solutions. The black circular markers indicate solutions that are shown in the inset panels whilst solid blue circular markers correspond to experimental results for the stable 2F
$_2$ state for bubbles of total area and area ratio equivalent to numerical bubbles volumes. (a) The solution space in the 
$(Q,U_{{b}})$ plane when 
$V_1+V_2 = 0.54{\rm \pi} ^2$, 
$V_{{r}} = 2/3$. The vertical dotted line indicates the value of 
$Q_{{s}}$. (b) The solution space in the 
$(V_{{r}},D_{{c}})$ parameter plane for 
$Q=0.04$ and the same total volume. The vertical dotted line indicates when the two bubbles are the same size, i.e. 
$V_{{r}}=1/2$.
Figure 12(a) shows that the stable 2F
$_2$ and unstable 2F
$_3$ solutions are connected through a limit point denoted L1. Unlike 2F
$_1$, the 2F
$_3$ unstable state appears to have no influence on the transient dynamics of the system according to figure 11(e). The limit point L1 occurs at a flow rate that coincides with the transition flow rate 
$Q_{{s}} \approx 0.0096$, discussed in figure 5, where the mathematical model predicts that the smaller bubble propagates faster than the larger one in their respective 1A
$_1$ single-bubble mode of propagation. This is consistent with the fact that a smaller leading bubble would then propagate faster in numerical time-dependent simulations for 
$Q< Q_{{s}}$ so that the two bubbles would separate without reaching a steady state. Surprisingly, although we know from § 4.1 that there is no evidence for 
$Q_{{s}}$ in the experiments, a threshold, consistent with a limit point, of the 2F
$_2$ state is also found experimentally at a critical flow rate 
$Q_{{L}1} = 0.0040 \pm 0.0002$ for a bubble area ratio 
$A_1/(A_1 + A_2) \approx 2/3$. A representative time evolution observed experimentally for 
$Q< Q_{L1}$ 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 2F
$_1$ 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 sidewall of the channel, see associated inset snapshot.
Time sequence snapshots of experiments where 2F
$_2$ is not reached. (a) 
$A_1/(A_1 + A_2) \approx 2/3$ (bubble ratio 
$2:1$) and 
$Q=0.002 < Q_{{L}1}$ and (b) 
$A_1/(A_1 + A_2) = 0.753$ and 
$Q=0.04$.
Figure 12(b) shows that the stable 2F
$_2$ and unstable 2F
$_1$ solutions, which were disconnected in the 
$(Q,U_{{b}})$ projection of figure 12(a), are in fact connected through a limit point denoted L2 under variations in volume ratio. This means that there is a critical volume ratio above which the stable 2F
$_2$ state does not exist, which we find to be 
$V_{{r}} \approx 0.79$ 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_1/(A_1 + A_2) > 0.744 \pm 0.007$ at the same flow rate. A representative experimental time evolution at 
$A_1/(A_1 + A_2) = 0.753$ is shown in figure 13(b) and equivalent numerical simulations are qualitatively similar. The overall dynamics is essentially the same as that at low 
$Q$ below L1 in the experiments: for values of the volume ratio above L2 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 2F
$_2$ and 2F
$_3$ branches. The chosen fixed length of the computational domain means that the 2F
$_2$ branch solution could not be calculated close to that limit. By contrast, the 2F
$_1$ 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$.
4.3. 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 (Keeler et al. Reference Keeler, Thompson, Lemoult, Juel and Hazel2019; Gaillard et al. Reference Gaillard, Keeler, Thompson, Hazel and Juel2021), 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_1+V_2 = {\rm \pi}0.54^2$ and the bubble volumes have a ratio 
$V_{{r}} = 2/3$. The bubble volume in the single-bubble system is the same as the volume of the larger bubble of the two-bubble system, i.e. 
$V=V_1$. A selection of the branches are illustrated by snapshots at given flow rates indicated by a circular marker.
Comparison of the single-bubble and two-bubble solution spaces plotted in the 
$(Q,U_{{b}})$ plane. The single-bubble space is denoted by coloured lines and the two-bubble system by black lines, with solid/dashed lines indicating stable/unstable branches. The 1A
$_1$ solution branch for the single-bubble system is shown for the larger (red) and smaller (blue) bubble. Solid markers indicate specific solutions on the branch shown in the inset panels. The vertical dotted line indicates 
$Q_{{s}}$. The volume of the single bubble is identical to the larger of the bubbles in the two-bubble system.
For 
$Q>Q_{{s}}$ the two-bubble 2A
$_2$ and 2F
$_3$ solution branches overlap and closely match the single-bubble asymmetric state, 1A
$_1$. Furthermore the two-bubble symmetric 2S
$_2$ solution branch is almost indistinguishable from the symmetric one-bubble state, 1S
$_2$/1S
$_3$. For 
$Q< Q_{{s}}$ there is no single-bubble state corresponding to 2A
$_1$, but the symmetric two-bubble 2S
$_1$ solution branch has a similar structure to the single-bubble 1S
$_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 two-bubble 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 2A
$_2$ and 2S
$_2$ and their single-bubble counterparts) but, as can be seen from the figure, the stable two-bubble 2F
$_2$ state is slightly slower than the corresponding 1A
$_1$ state corresponding to leading (smaller) bubble.
Finally, we would expect the position of the two-bubble bifurcation points, as 
$h$ and 
$V_{i}$ are varied, to be similar to the single-bubble equivalent, which have been calculated previously (Keeler et al. Reference Keeler, Thompson, Lemoult, Juel and Hazel2019). For the two-bubble system the calculations are more computationally expensive and we therefore leave a detailed analysis of the loci of the bifurcation points with parameter variations to a future study.
5. 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 (Gaillard et al. Reference Gaillard, Keeler, Thompson, Hazel and Juel2021). 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 long-term 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 (Maxworthy Reference Maxworthy1986; Madec Reference Madec2021), in which the multiple bubbles will always separate or aggregate. There are large regions of overlap between particular single- and two-bubble 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 2F
$_3$ and 2A
$_2$ have the same speed as the stable asymmetric single bubble 1A
$_1$; and the unstable symmetric two-bubble state 2S
$_2$ has the same speed as the stable symmetric single-bubble 1S
$_3$. 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 (Beatus, Tlusty & Bar-Ziv Reference Beatus, Tlusty and Bar-Ziv2006; Beatus, Bar-Ziv & Tlusty Reference Beatus, Bar-Ziv and Tlusty2012) 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.
Experimental evidence of a multi-bubble stable steady states, for 
$Q = 0.011$. (a) A three-bubble stable steady state (dimensionless areas 
$0.33^2{\rm \pi}, 0.40^2{\rm \pi}, 0.47^2{\rm \pi}$). (b) A four-bubble stable steady state (dimensionless areas 
$0.33^2{\rm \pi}, 0.40^2{\rm \pi}, 0.47^2{\rm \pi}, 0.54^2{\rm \pi}$).
Finally, although there is not a direct equivalent of the 2F
$_1$ state in the single-bubble system, the 2F
$_1$ 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 increase indefinitely as 
$h\to 0$, and, as far as we know, no further bifurcations occur, but the 2F
$_1$ state barely changes. The question of whether a stable two-bubble steady state exists when 
$h = 0$ remains open but certainly if one does exist it 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 2F
$_2$ solution (an ‘offset’ configuration with the smaller bubble ahead) and that the dynamics is incredibly slow, indicating neutral stability. This exploration of two-bubble steady states in the experiment will form part of a future investigation.
Funding
We acknowledge funding from the Engineering and Physical Sciences Research Council through and grant numbers EP/P026044/1 EP/T021365.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available in figshare at https://doi.org/m9.figshare.17277401.
Author contributions
J.S.K. performed all numerical simulations; A.G. performed the experiments in figures 6 and 9; J.L. performed the experiments in figures 5(c), 11, 12, 13 and 15. A.L.H, A.J and A.B.T conceived the study. All authors contributed equally to analysing data and reaching conclusions, and in writing the paper.
Appendix A. Experimental protocols
A.1. Experimental protocol: aligned bubbles
Aligned pairs of bubbles were obtained by producing two bubbles sequentially with the trailing bubble generated once the leading bubble had propagated by a prescribed distance from the air port, which was set by the volume of oil injected in the intervening time. A dimensionless flow rate 
$Q_i$ was then imposed to propagate the pair of bubbles to a position a few centimetres downstream of the centring device before interrupting the flow for half a second. For 
$Q_i=0.029$, the two bubbles straddled the rail symmetrically about the channel centreline. However, asymmetric bubbles which overlapped the rail from only one side could also be obtained for sufficiently small values of the flow rate (
$Q_i = 0.0056$) because of the absence of stable one-bubble symmetric modes of propagation, which meant that initially centred bubbles migrated off the rail. This migration was always to the same side of the rail because of unavoidable bias in the levelling of channel; see Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021) for details.
A.2. Experimental protocol: offset bubbles
The simplest method to generate bubbles on opposite sides of the rail is to break up a single bubble into two parts, as illustrated in figure 2(b). This was achieved by initially propagating a single symmetric bubble downstream of the centring device at a flow rate 
$Q_i = 0.03$ before interrupting the flow for a controlled duration to allow the bubble to widen and initiate its sideways migration towards one of the deeper regions of the channel, see Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021) for a description of the bubble relaxation process at 
$Q = 0$. A flow rate of 
$Q = 0.05$ was then imposed which led the bubble to break into two parts of different sizes, overlapping the rail from opposite sides. The difference in bubble sizes after breakup was set by the 
$y$-offset of the single bubble upon imposing the flow. Depending on the direction of the oil flow used to break up the initial bubble (either from the inlet towards the outlet or vice versa), different initial bubble-pair configurations could be prepared, with the larger bubble either leading (closer to the outlet) or trailing (closer to the inlet). The initial separation between the two bubbles was controlled by imposing a flow rate of 
$Q = 0.007$ in the same direction as the flow used to split the single bubble. This flow rate was chosen for two reasons. Firstly, at this flow rate, the larger leading bubble propagated faster than the smaller trailing bubble so that the distance between the two bubbles increased with time. Secondly, this flow rate was also sufficiently small to prevent the smaller bubble from crossing the rail because of its hydrodynamic interaction with the larger bubble. At lower flow rates, the relative strength of the capillary force acting on each bubble increased which in turn meant that they propagated more asymmetrically, ensuring that both bubbles remained on their respective sides of the rail. Finally, when the bubbles reached a prescribed separation, the flow was interrupted for half a second before beginning the two-bubble experiment.
Appendix B. Numerical methods
B.1. Topology changes
Figure 8 shows the numerous topology changes the system of bubbles experiences. To facilitate topology change (i.e. bubble coalescence and bubble breakup) in the numerical code we employ the procedure described in Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021). After each time step, we measure the minimum distance between each of the pairs of bubbles and also check if each individual bubble has self-intersected. If the minimum distance between two bubbles is lower than a pre-defined threshold (in these calculations we choose 
$10^{-2}$), the bubbles are merged into a single bubble. Alternatively, if self-intersection has been identified we split the bubble into two separate bubbles. In each case, the simulations continue after the topology change, with the number of volume constraints in the system deleted/added, as appropriate. For more details we refer the reader to Gaillard et al. (Reference Gaillard, Keeler, Thompson, Hazel and Juel2021).
B.2. Interval-bisection algorithm
The method for calculating the unstable two-bubble aligned states is now described in more detail. The algorithm is initiated by solving two IVPs, one in the case where the two bubbles are initially sufficiently far from one another to separate and one where they are sufficiently close to one another to ultimately coalesce, as illustrated in figure 7. After this initial step, a new IVP is solved where the initial distance between the bubbles is
\begin{equation} D_{{edge}}(t = 0) = (D_{{s}} + D_{{c}})/2, \end{equation}
where 
$D_{{s}}$ and 
$D_{{c}}$ are the values of the initial bubbles distance in the previous simulations that lead to separation and coalescence respectively. Once the final dynamical outcome 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)\to 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 criterion for the interval-bisection procedure is that the bubbles remain within a small distance, 
$\varepsilon$, of each other for a large time, 
$T$, i.e.
\begin{equation} |D(t) - D_{{edge}}(t=0)|<\varepsilon,\quad \forall t< T_{{end}}. \end{equation}
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 
$\varepsilon = 0.1$ and 
$T=20$ is sufficient to ensure that we converged to a steady state.
Open access
